A Likely Detection of a Two-planet System in a Low-magnification Microlensing Event

Multiple-planet systems are not uncommon in the planetary systems found by Kepler (Lissauer et al. 2011a; Fabrycky et al. 2014). Of the planetary systems discovered by the transit method, 22.4% are multiple-planet systems.²¹ Transiting multiple systems provide much more information than transiting single systems, i.e., mutual inclination, orbital period ratio, and average density with the mass measurement by transit timing variation. Such rich information allows us to constrain the planet formation models. Tightly packed coplanar multiple systems, such as Kepler-11 (Lissauer et al. 2011b) and TRAPPIST-1 (Gillon et al. 2017), are thought to form in outer regions rather than in the current orbit and migrate inward through a disk. Transiting single-planet systems could be dynamically hotter and less orderly compared to the transiting multiple-planet systems (Morton & Winn 2014). These studies have rapidly developed our understanding of the architecture of close-in, hot multiple-planet systems, but there is a huge parameter space waiting for our probes: cold, wide-orbit multiple-planet systems.

The radial velocity (RV) method is sensitive to planets with wider orbits than transit, and it has detected a lot of multiple-planet systems. The ratio of the RV multiple-planet systems to the number of planetary systems in the RV sample is 22.2%¹, which is, interestingly, almost same as the transit sample. However, the inclination of each planet is not able to be derived with the RV alone. Also, finding wide-orbit planets with a mass of below Saturn mass is difficult with current facilities. Direct imaging (DI) can find multiple-planet systems if the systems have young, bright, massive objects in very wide (mostly ≥10 au) orbits (Marois et al. 2008), but again, detecting planets with sub-Jupiter mass or below is very difficult with current instruments. The ratio of the number of multiple systems to the number of whole systems detected by DI is 4.2%¹, which is probably reflecting the intrinsic low occurrence rate of multiple super-Jupiters in wide orbit.

Microlensing does have a sensitivity down to Earth-mass planets (Bennett & Rhie 1996) orbiting beyond the snow line, at which ices, including water ice, condense to increase the surface density of solid materials by a factor of five in the protoplanetary disks (Ida & Lin 2005; Kennedy & Kenyon 2008). Massive planets with very wide orbits (Bennett et al. 2012; Poleski et al. 2014), even for unbound ones (Sumi et al. 2011; Mróz et al. 2017b), are detectable by microlensing. Thus, in principle, microlensing is the best tool to study the cold multiple-planet systems, especially for small planets, which are not able to be probed by the other methods (including the astrometry by the Gaia mission). The host star–to–planet mass ratio and normalized projected separation can always be directly measured from observed microlensing light curves. The planet mass and physical projected separation are also measurable for about half of planetary systems with combinations of measurements of the finite-source effect, Earth orbital parallax effect (Muraki et al. 2011), satellite parallax (Udalski et al. 2015c; Street et al. 2016), and high angular resolution images for the lens flux measurement (Fukui et al. 2015; Koshimoto et al. 2017b) and resolving the lens and source (Batista et al. 2015; Bennett et al. 2015). However, measuring the orbital parameters could be difficult for most of the planetary microlensing events. Nevertheless, the inclination and eccentricity of the first multiple planets found by microlensing, OGLE-2006-BLG-109Lbc, were constrained with the assumption of a coplanar orbit (Bennett et al. 2010).

So far, only two multiple-planet systems have been found by microlensing against 49 planetary systems, so the observed multiple ratio is 4.1%¹, which is almost same as the DI sample and much smaller than the transit and RV results. There are two reasons for the low observed multiple ratio in the microlensing sample. First, the rate of very high-magnification events, which is the major microlensing channel to detect multiple-planet systems, is low, although the detection efficiency through this channel is high (Gaudi et al. 1998). Second, the detection efficiency for multiple Neptune-mass (or below) planet systems is not high, and the occurrence rate of such systems would also not be very high. (We do not yet know about the occurrence rate for multiple Neptune-mass planets.) The two multiple-planet systems were found through the high-magnification channel. OGLE-2006-BLG-109L is a system of Jupiter and Saturn analogs with orbits approximately half of the Jupiter and Saturn orbits, respectively, around the host with half of the solar mass (Gaudi et al. 2008; Bennett et al. 2010). The second multiple-planet system is OGLE-2012-BLG-0026L (Han et al. 2013), which is a sub-Jupiter and sub-Saturn system with a solar-mass host (Beaulieu et al. 2016). For both systems, the masses were measured with the robust parallax-effect measurement and further constraints from the high angular resolution images. By contrast, low-magnification events are expected to have a detection efficiency of ${ \mathcal O }(1 \% )$ or less²² to detect multiple (two-planet) systems (Gaudi 2012).

In this paper, we report the analysis of OGLE-2014-BLG-1722, a low-magnification microlensing event showing two planetary anomalies that were caused by a lens system having two Saturn-mass planets. The event timescale of this event is ∼24 days, which is a typical value of microlensing events observed toward the Galactic bulge. The mass of the lens system, M, and distance to the lens system, D_L, are related to the event timescale by the relation

$$\begin{eqnarray}&&{t}_{{\rm{E}}}=\displaystyle \frac{{\theta }_{{\rm{E}}}}{{\mu }_{\mathrm{rel}}},\,\,\,\,\,{\theta }_{{\rm{E}}}=\sqrt{\displaystyle \frac{4\mathrm{GM}}{{c}^{2}}\left(\displaystyle \frac{1}{{D}_{{\rm{L}}}}-\displaystyle \frac{1}{{D}_{{\rm{S}}}}\right)},\end{eqnarray} \tag{ 1 }$$

where μ_rel is the relative lens-source proper motion and D_S is the distance to the source star. For most planetary events, we can measure the finite-source effect (or source crossing time), leading to the measurement of the angular Einstein radius θ_E. Furthermore, if we measure the microlensing parallax ${\pi }_{{\rm{E}}}\equiv {\pi }_{\mathrm{rel}}/{\theta }_{{\rm{E}}}$, ${\pi }_{\mathrm{rel}}=\mathrm{au}({D}_{{\rm{L}}}^{-1}-{D}_{{\rm{S}}}^{-1})$, we can robustly measure the mass of and distance to the lens system, assuming the source distance (∼8 kpc). Unfortunately, neither of these effects were measured in this planetary microlensing event; thus, we estimated the lens property with the Bayesian approach. We describe the observation and data reduction in Section 2. In Section 3, we show our light-curve modeling. In Section 4, we estimate the physical parameters: the distance to the system, mass of the host star and two planets, and projected separation between the host and planets based on the Bayesian analysis. We briefly discuss the validity of the second anomaly signal in Section 3.2. In Section 5, we estimate the occurrence rate of multiple cold gas giant planet systems. In Section 6, we summarize the results and conclude.

The microlensing event OGLE-2014-BLG-1722 was alerted by the Optical Gravitational Lens Experiment (OGLE; Udalski 2003) collaboration using their real-time event detection system, the Early Warning System, at HJD' ≡ HJD–2450000 = 6889.1²³ as a new microlensing event at (R.A., decl.)₂₀₀₀ = (17^h55^m0057, −31°28'086). The event occurred in the OGLE-IV field (Udalski et al. 2015b) BLG507, which is observed one to three times per night using the 1.4 deg² camera on the 1.3 m Warsaw Telescope at the Las Campanas Observatory in Chile. This microlensing event is also in the footprint of the MOA (Bond et al. 2001; Sumi et al. 2003) survey, which uses the MOA-II 1.8 m telescope with a 2.2 deg² field-of-view CCD camera, MOA-cam3 (Sako et al. 2008), at Mt. John University Observatory in New Zealand. Following the OGLE alert, the MOA team discovered this event independently at HJD' = 6890.0 and labeled it as MOA-2014-BLG-490. Since this event is found in one of the highest-cadence MOA fields, the images were taken more than 20 times per night. OGLE data were taken with the I band, whereas the MOA data were taken with the custom wide-band MOA-Red filter, which is roughly the sum of the Cousins R and I bands. The observed light curves of OGLE and MOA are plotted with blue and red circles in Figure 1, where the MOA data points are binned into 1.2 hr.

Zoom In Zoom Out Reset image size

No one noticed in real time that the first anomaly occurred a few days before the discovery alerts, as well as the second, weak anomaly at HJD' = 6899. The first anomaly was found about 20 days after the magnification peak. Even at the peak, the event was not bright enough for observations using small telescopes. Also, when the anomaly was found, it was clear that it had ended, and the entire light curve was well covered by the OGLE and MOA teams. Thus, no observations were conducted by the follow-up teams, and the OGLE and MOA survey observations were conducted with the normal cadence throughout the event. As a result, this event can be classified as a planetary event discovered by pure surveys (Koshimoto et al. 2014; Shvartzvald et al. 2014; Suzuki et al. 2014; Yee et al. 2012; Rattenbury et al. 2015; Mróz et al. 2017a), which is a main stream of the current microlensing planet detection opposed to the survey+follow-up discoveries in the early time of microlensing observation.

The OGLE and MOA data were first reduced by using their difference image analysis (DIA) photometry pipelines (Bond et al. 2001; Udalski 2003). After OGLE-2014-BLG-1722 was found to be a planetary event, the OGLE and MOA data were carefully rereduced to precisely characterize the anomalies. The optimized centroid of the event position was used for the rereduced OGLE light curve. For the MOA images, we use a numerical kernel (Bramich 2008) with a modification to allow for a spatial variation of the kernel across the subimages centered on the event. Also, we correct seeing, airmass, and differential refraction correlations in the photometry, which can often be seen for the MOA images. OGLE photometry is calibrated to the OGLE-III photometry map (Szymański et al. 2011). For the MOA data, instrumental magnitude is used for the light-curve modeling below.

There could be possible systematics in the long baseline in the OGLE and MOA data, which can sometimes affect the microlensing parallax measurements and other parameters. To avoid this systematics, we use the photometry data in only 2014 and 2015 for the light-curve modeling. In general, the error bars given by the photometry codes for the crowded region are underestimated, which avoids us finding the correct light-curve model and estimating the accurate uncertainties for each fitting parameter. Thus, we normalize the error bars by using ${\sigma }_{i}^{{\prime} }=k\sqrt{{\sigma }_{i}^{2}+{e}_{\min }^{2}}$, where σ_i is the original error bar of the ith data point in magnitudes, and k and e_min are the renormalizing parameters (Yee et al. 2012). We expect that the cumulative of χ² in each data point sorted with the magnification of the best-fit model is a straight line if the data are normal distribution. The e_min parameter is important for high-magnification events where the error bars of bright points are especially underestimated. In such a case, the cumulative of χ² would have a break, and the e_min factor corrects it. We choose e_min to make the cumulative χ² distribution a straight line, and k is chosen so that each data set gives ${\chi }^{2}/\mathrm{dof}=1$. We use k = 1.305371 and e_min = 0.00 for the OGLE data and k = 1.139291 and e_min = 0.00 for the MOA data for the final result.

To double-check the second anomaly feature, a different point-spread function (PSF) is also used for the MOA data reduction. For this purpose, we reduced the MOA images by using DIA where we construct the PSF with the DoPHOT package (Schechter et al. 1993; Bennett et al. 2012). The generated light curves with the DoPHOT PSF are usually used to measure the source-star color, with which the angular Einstein radius θ_E can be estimated if we measure the source crossing time t_*. In this event, however, the source crossing time was not able to be measured, as the source trajectory did not cross the caustics or approach the cusps. So, we use the DoPHOT PSF light curves only for light-curve modeling as written in Section 3.2.

3.1. One-planet Model

The light curve of OGLE-2014-BLG-1722 is almost symmetric except for the short-term feature at HJD' ∼ 6887, a dip lasting for a few days (see Figure 1). Considering the duration of the dip, the feature at a wing of the light curve is expected to be caused by a small mass ratio companion to the lens located inside of the Einstein radius. Nevertheless, we investigate the vast parameter space to find the best-fit model for the observed light curve.

A microlens light curve produced by a single mass is described by three parameters: the time of the peak magnification, t₀; the event timescale, t_E, which is the time required for the Einstein radius crossing time; and the impact parameter divided by the angular Einstein radius, u₀. If the lens object has a companion, three additional parameters are required: the mass ratio of the planet to the host star, $q(\equiv {m}_{{\rm{b}}}/{M}_{\mathrm{host}})$; the projected separation divided by the angular Einstein radius, s; and the angle between the source trajectory and planet–star axis, α. Most planetary events are characterized by the strong features of caustic crossings, at which the magnification of a point source diverges. But the actual source star has a finite angular size, so the magnification is limited such that a sharp magnification change is smoothed out. Thus, modeling this rounded feature in a light curve allows us to measure the relative source-star radius with an additional parameter, the angular source-star radius, relative to the angular Einstein radius, ρ, or the source radius crossing time, ${t}_{* }\equiv \rho {t}_{{\rm{E}}}$.

In addition to the above seven parameters, two linear parameters for each photometry datum are required to describe the light curve: the flux from the source star, F_S, and from the blended stars, F_B, which is the sum of the flux from the lens and adjacent stars on the sky that are unrelated to the microlensing event. Thus, the light-curve model is given by $F(t)=A(t){F}_{{\rm{S}}}+{F}_{{\rm{B}}}$, where F(t) is the flux at a given time t and $A(t)=A(t;{t}_{0},{t}_{{\rm{E}}},{u}_{0},q,s,\alpha ,\rho )$ is the magnification. The data set observed with different filters and telescopes can be aligned by the linear fitting, yielding F(t), F_S, and F_B for each different data set.

We use a variation of the Markov chain Monte Carlo (MCMC) algorithm (Verde et al. 2003) to search for the best-fit model efficiently (Bennett 2010) with the rereduced photometries. We start to explore the parameter space with a grid for the parameters q, s, α that are sensitive to anomaly signals. We search for the other parameters by a downhill approach. The range of the grid for the mass ratio, q, is 10⁻⁴–1, which well covers the planetary and stellar binary mass ratio regime. As for the separation, we investigate $-0.5\lt \mathrm{log}\,s\lt +0.5$ because microlensing is sensitive to the companions in this region. We use 22, 11, and 40 grids for $\mathrm{log}s,\mathrm{log}q$, and α, respectively. After the grid search, we selected the best 100 models as initial models to run the fitting code by letting each parameter free.

The best-fit one-planet model we found²⁴ has a mass ratio of q ∼ 4 × 10⁻⁴, which corresponds to a sub-Saturn planet if we assume that the mass of the host star is 0.5 M_⊙. The separation is s ∼ 0.76, which agrees with what we expect from the dip feature at the light-curve wing. The χ² improvement for this planetary model over the single-lens model is Δχ² = 1029, so the planetary signal was significantly detected by the two survey teams. But, unfortunately, the relative source-star radius, ρ, cannot be measured because the source did not cross the caustics. The Einstein radius crossing time of t_E = 24 days is a typical value for single-lens events (Suzuki et al. 2016, Figure 1), and it is usually too short to measure the parallax effect that is often measured in events with t_E ≳ 60 days. Nevertheless, our model includes the parallax-effect parameters ${\pi }_{{\rm{E}}}=({\pi }_{{\rm{E}},{\rm{E}}},{\pi }_{{\rm{E}},{\rm{N}}})$. The parameters of the best-fit one-planet model are shown in Table 1.

Table 1.  Best-fit Parameters

Parameters	Best Close	MCMC Close	Best Wide	MCMC Wide	Best One-planet	MCMC One-planet
χ2	4191.0	⋯	4192.2	⋯	4361.7	⋯
${t}_{0}\ (\mathrm{HJD}^{\prime} )$	6900.224	${6900.228}_{-0.008}^{+0.009}$	6900.215	${6900.219}_{-0.011}^{+0.012}$	6900.191	${6900.191}_{-0.010}^{+0.010}$
${t}_{{\rm{E}}}\ (\mathrm{days})$	23.819	${23.797}_{-0.380}^{+0.418}$	23.867	${23.909}_{-0.582}^{+0.801}$	24.211	${24.234}_{-0.512}^{+0.513}$
u0	−0.131	$-{0.132}_{-0.003}^{+0.003}$	−0.131	$-{0.131}_{-0.003}^{+0.004}$	0.127	${0.127}_{-0.003}^{+0.003}$
$q\ ({10}^{-4})$	4.468	${4.451}_{-0.557}^{+0.693}$	4.146	${4.263}_{-0.607}^{+0.678}$	4.335	${4.534}_{-0.583}^{+0.705}$
s	0.754	${0.754}_{-0.004}^{+0.004}$	0.754	${0.754}_{-0.005}^{+0.006}$	0.757	${0.757}_{-0.004}^{+0.004}$
$\alpha \ (\mathrm{rad})$	0.228	${0.230}_{-0.024}^{+0.020}$	0.239	${0.242}_{-0.032}^{+0.035}$	0.209	${0.205}_{-0.025}^{+0.027}$
ρ	0.005	${0.007}_{-0.005}^{+0.006}$	0.002	${0.006}_{-0.004}^{+0.005}$	0.007	${0.008}_{-0.005}^{+0.006}$
${\pi }_{{\rm{E}},{\rm{N}}}$	0.199	${0.189}_{-0.553}^{+0.570}$	−0.077	$-{0.153}_{-0.924}^{+0.843}$	−0.585	$-{0.665}_{-0.595}^{+0.655}$
${\pi }_{{\rm{E}},{\rm{E}}}$	0.092	${0.094}_{-0.083}^{+0.101}$	0.103	${0.104}_{-0.101}^{+0.110}$	0.002	$-{0.010}_{-0.083}^{+0.086}$
${q}_{2}\ ({10}^{-4})$	6.388	${7.071}_{-1.832}^{+2.371}$	6.245	${7.098}_{-1.855}^{+1.836}$	⋯	⋯
s2	0.851	${0.843}_{-0.030}^{+0.034}$	1.358	${1.368}_{-0.044}^{+0.044}$	⋯	⋯
$\psi \ (\mathrm{rad})$	2.196	${2.190}_{-0.033}^{+0.028}$	2.207	${2.209}_{-0.036}^{+0.036}$	⋯	⋯
${I}_{{\rm{S}}}$ (mag)	18.88	18.87 ± 0.03	18.87	18.88 ± 0.04	18.92	18.92 ± 0.04
${f}_{{\rm{S}}}$	1.16	1.16 ± 0.03	1.16	1.16 ± 0.04	1.11	1.11 ± 0.04

Note. Here $\mathrm{HJD}^{\prime} \equiv \mathrm{HJD}-2450000$, and ${f}_{{\rm{S}}}={F}_{{\rm{S}}}/({F}_{{\rm{S}}}+{F}_{{\rm{B}}})$ is the blending parameter. The Best Close values are plotted as plus signs in Figure 3.

Download table as:  ASCIITypeset image

Note that the best-fit one-planet model shows negative blending (${f}_{{\rm{S}}}={F}_{{\rm{S}}}/({F}_{{\rm{S}}}+{F}_{{\rm{B}}})\gt 1$ in Table 1). As discussed in many papers (Park et al. 2004; Sumi et al. 2006; Smith et al. 2007), negative blending does happen for combinations of the following reasons: an incorrect light-curve model, systematic errors in the photometry, and/or, indeed, negative blending flux within the PSF. The last case occurs when the target happens to be located at a "hole" in the uneven background toward the crowded Galactic bulge fields. The uneven background consists of the mottled distribution of upper-main-sequence and turnoff stars in the bulge that are just below the detection limit of the ground-based microlens survey. If the background close to the event is lower than the average background, then the photometry pipelines will oversubtract the background, yielding a negative blending flux. In this event, the better light-curve models discussed later also show negative blending. Furthermore, the amount of blending, f_S ∼ 1.1–1.2, is quite common (Sumi et al. 2006; Smith et al. 2007), and the magnitude of the hole, which is given as ${I}_{\mathrm{hole}}={M}_{\mathrm{zero}\mathrm{point}}-2.5\,\mathrm{log}(| {F}_{{\rm{B}}}| )$, is consistent with unresolved turnoff stars in the Galactic bulge written in later. Thus, it is unlikely that the incorrect light-curve model yields the negative blending.

With a careful visual inspection of the residual plot from the best-fit one-planet model, a weak bump in the MOA data just before the peak magnification was found. Although this bump could be explained by systematic scatter, red noise, or other artificial effects, they are unlikely, as discussed in Section 3.2. We conduct further light-curve analysis in the following sections to examine whether the bump is explained by an additional companion to the lens or source star.

3.2. The Second Anomaly Has an Astrophysical Origin

3.3. The Planetary Model Ignoring the First Anomaly

3.4. Two-planet Model

A triple-lens model requires three additional fitting parameters: the mass ratio of the second planet to the host star, ${q}_{2}\,(\equiv {m}_{{\rm{c}}}/{M}_{\mathrm{host}})$; the projected separation between the second planet and host star normalized to the angular Einstein radius, s₂; and the angle between the line connecting the host star and first planet and the line connecting the host star and second planet, ψ. Thus, the magnification at a given time t is described as $A(t)=A(t;{t}_{0},{t}_{{\rm{E}}},{u}_{0},q,s,\alpha ,\rho ,{q}_{2},{s}_{2},\psi )$. As we did for the one-planet model, MCMC is used to localize the χ² minimum for the triple-lens model. We use the values of the one-planet model parameters for the initial parameters for (${t}_{0},{t}_{{\rm{E}}},{u}_{0},\alpha ,\rho$). Considering the lens system and source trajectory configuration from the one-planet model and planetary model ignoring the first anomaly, we estimate that ψ is roughly 2.1 (120^◦) and use this as an initial value. For the parameters of ($q,s,{q}_{2},{s}_{2}$), we use 162 combinations of the initial values, which consist of three, three, three, and six different values for each parameter. Note that the s₂ covers 0.7 < s₂ < 1.6 to investigate the possible close–wide degeneracy.

As we see in the one-planet model, a significant detection of the parallax effect is unlikely in this event, but again, we include the parameters, since including the parallax parameters could affect the uncertainties of the other parameters (Skowron et al. 2016), and the upper limit could be used as a constraint to estimate the mass of and distance to the lens system.

The best-fit two-planet model is plotted in Figure 1, and its parameters are shown in Table 1. The geometry of the lens system and source trajectory is plotted in Figure 2. As we expected in the previous subsection, we found close (s₂ < 1) and wide (s₂ > 1) models for the second planet. The parameters for the two models are consistent with each other, with the exception of s₂. The parameters of the primary lens and first planet in this model are also consistent with those of the best one-planet model. The χ² improvement of the best-fit two-planet model (close model) from the best one-planet model is 171, which is significant if we apply the detection threshold of ${\chi }_{\mathrm{thres}}^{2}=100$. Here 80% of the χ² improvement originates from the 24 data points around HJD' = 6899, 15% comes from the data at HJD' = 6900, and the rest is from the nights around the peak magnification. Only the MOA data contribute to the detection of the second planet, OGLE-2014-BLG-1722Lc. Note that DPB conducted a grid search for the third lens object with all parameters for the primary star and first planet fixed and found the best model, which is consistent with the above best-fit two-planet model.

Zoom In Zoom Out Reset image size

The amount of negative blending in the two-planet model is almost the same level as that in the one-planet model, as written in Table 1. The magnitude of the "hole" is ${I}_{0,\mathrm{hole}}=19.27\,\pm 0.23$ with the extinction correction (see Section 4), which is consistent with the unresolved turnoff stars in the Galactic bulge.

Figure 3 shows parameter correlations for the MCMC chain for the best-fit close model. The mass ratio of planet b, q, is measured with 16% uncertainty, but the mass ratio for planet c, q₂, has a large uncertainty of ∼33%. This is simply because the anomaly of planet c in the light curve is relatively weak. Nevertheless, the projected separations of both planets are well measured, as the separation is very sensitive to the size of the caustics compared to the mass ratio (Chung et al. 2005; Han 2006). Unlike the previous two-planet microlensing events (Gaudi et al. 2008; Han et al. 2013), the relative source-star size and Earth's orbital parallax effect are not well constrained, because the source trajectory did not hit the caustics and/or cusp, and the event timescale is not long enough to measure the parallax. Thus, we can estimate the mass and distance of the lens system only from a Bayesian analysis based on a standard Galactic model and the assumption that the planet hosting probably does not depend on the mass and distance of the lens star (Beaulieu et al. 2006; Suzuki et al. 2014).

Zoom In Zoom Out Reset image size

3.5. Alternative Models

3.5.1. Stellar Companion to the Lens

It is unlikely that a stellar companion to the lens (i.e., a circumbinary planetary system or a planet orbiting one member of a binary system) causes the second anomaly, because the stellar binary solutions were ruled out by Δχ² > 220 when we conducted the binary-lens modeling excluding the data points around the first anomaly.

3.5.2. Binary-source Model

A binary-source model is also a possible alternative explanation of the second anomaly (Gaudi 1998; Jung et al. 2017b). In this event, a hypothetical second source star is expected to be very faint and highly magnified, since the apparent duration of the anomaly is very short (∼1 day) compared to that of the primary event (∼1 month), and we do not see any excess flux when the primary source is demagnified (compared to the single-lens–single-source model) at the first anomaly. To test this binary-source scenario, we fit the light curves with a binary-source model (more specifically, a binary-lens–plus–binary-source model), which is also used for the modeling of MOA-2010-BLG-117, the first²⁶ planetary microlensing event with a binary source (Bennett et al. 2018), and OGLE-2016-BLG-1003, the first resolved caustic-crossing binary-source event (Jung et al. 2017a). Compared to a binary-lens model with a single source star, the binary-source model requires an additional three parameters: the time of the peak magnification for the second source, t_0,2; the impact parameter of the second source divided by the Einstein radius, u_0,2; and the flux ratio of the second source to the primary source, f_ratio. We found that the best-fit binary-source model is also able to explain the small bump. The best-fit parameters are given in Table 2 with the apparent brightness of the source stars and blending parameter. Note that there is a twofold degeneracy in u_0,2: the positive and negative values in u_0,2 show the same magnification. These degenerate models yield different values of projected separation between the two sources. The amount of the blending is also very similar to that of one-planet and two-planet models. The best-fit binary-source model is slightly disfavored compared to the best two-planet model by Δχ² of 4.2 (2.0σ). Note that we use the static model (no microlensing parallax parameters) for the χ² comparison, since the two-planet static model and the one-planet binary-source model have the same number of modeling parameters. Thus, we cannot simply rule out the binary-source scenario from the light-curve fitting.

Table 2.  Parameters for the Binary-source Model

Parameters	Values
${\chi }^{2}$	4196.7
${t}_{0}\ (\mathrm{HJD}^{\prime} )$	${6900.244}_{-0.001}^{+0.001}$
${t}_{{\rm{E}}}\,(\mathrm{days})$	${23.893}_{-0.227}^{+0.202}$
u0	0.133 ± 0.002
$q\ ({10}^{-4})$	${4.218}_{-0.267}^{+0.285}$
s	0.753 ± 0.002
$\alpha \ (\mathrm{rad})$	0.238 ± 0.002
$\rho \ ({10}^{-3})$	${0.709}_{-0.272}^{+0.555}$
${t}_{\mathrm{0,2}}\ (\mathrm{HJD}^{\prime} )$	${6898.960}_{-0.003}^{+0.006}$
${u}_{\mathrm{0,2}}({10}^{-3})$	${7.000}_{-0.235}^{+0.216}$
${f}_{\mathrm{ratio}}({10}^{-3})$	${3.575}_{-1.160}^{+0.320}$
${I}_{{\rm{S}}1}(\mathrm{mag})$	18.88 ± 0.01
${I}_{{\rm{S}}2}(\mathrm{mag})$	25.00 ± 0.23
${f}_{{\rm{S}}12}$	1.16 ± 0.02

Note. Here $\mathrm{HJD}^{\prime} \equiv \mathrm{HJD}-2450000$, I_S1 and I_S2 are the apparent magnitudes of the primary and companion source stars, and ${f}_{{\rm{S}}12}={F}_{{\rm{S}}12}/({F}_{{\rm{S}}12}+{F}_{{\rm{B}}})$, where F_S12 is the total flux from the source system, is the blending parameter.

Download table as:  ASCIITypeset image

However, the occurrence rate of the stellar companion to explain the best binary-source model is not very high. Assuming the projected Einstein radius on the source plane, ${\hat{r}}_{{\rm{E}}}=2.93\,\mathrm{au}$, that is estimated using the Galactic model (Han & Gould 1995) with a constraint of t_E, the projected physical separation of the binary source is ${0.40}_{-0.15}^{+0.20}$ and ${0.44}_{-0.17}^{+0.21}\,\mathrm{au}$ for the positive and negative u_0,2, respectively. Although we do not have color information for the primary source star, it is safe to assume that the primary is a solar-type star. With an assumption that the period distribution of nearby solar-type stars (Raghavan et al. 2010) is applicable to the stars in the Galactic bulge, the estimated projected physical separation between the two sources corresponds to the tail of the period distribution of companion stars in Raghavan et al. (2010). To fairly compare the relative probability of the two-planet and binary-source scenario, we also consider the χ² difference, the prior probability from the mass and separation distribution, and the detection efficiency.

For the binary-source scenario, we use the companion distribution of Raghavan et al. (2010) to estimate the prior probability. For the two-planet scenario, the mass ratio function of Suzuki et al. (2016) is used. The relative prior probability we obtained is 0.041 and 0.043 for the positive and negative u_0,2 binary-source models and 0.152 and 0.192 for the close and wide s₂ two-planet models. The detection efficiency we obtained is 0.098 and 0.052 for the positive and negative u_0,2 binary source and 0.46 and 0.12 for the close and wide s₂ two-planet models, respectively. With the above values, we found a relative probability of 0.991 and 0.009 for the two-planet and binary-source scenarios, respectively. The binary-source scenario is marginally ruled out by 3.1σ, if we consider that the relative probability is equivalent with the weight of ${e}^{-{\rm{\Delta }}{\chi }^{2}/2}$.

The microlensing method allows us to measure the mass of and distance to the lens system when we measure the microlensing parallax and finite-source effects (Gould 2000). It is also possible to measure the mass and distance of the lens if we can determine the flux of the lens star in addition to a measurement of the microlensing parallax or finite-source effect (Bennett et al. 2010, 2015; Batista et al. 2015). But it is important to verify that the excess flux at the position of the source is really due to the lens (Bhattacharya et al. 2017; Koshimoto et al. 2017a). Unfortunately, OGLE-2014-BLG-1722 showed neither the parallax-effect nor the finite-source effect in the observed light curve. Thus, we need to use the Bayesian approach to estimate the mass of the lens system. As a prior, we use the Galactic model of Han & Gould (1995) and the mass function of Sumi et al. (2011). The constraints we use here to obtain the posterior distribution of the lens mass and distance are the observed t_E and upper limit of π_E.

The estimated physical values of the lens system are given in Table 3. As we expect, the physical parameters of the close and wide solutions are similar to each other except for the projected/deprojected separation between the host and planet c (i.e., ${r}_{\perp ,{\rm{c}}}$ and ${a}_{3{\rm{D}},{\rm{c}}}$). Thus, we combine the posterior distribution of close and wide models with a weight of ${e}^{-{\rm{\Delta }}{\chi }^{2}/2}$, where Δχ² is the χ² difference of the two models, Δχ² = 1.2, and report the combined values of each parameter except for ${r}_{\perp ,c}$ and ${a}_{3{\rm{D}},{\rm{c}}}$. The posterior distribution of the distance to the lens system and mass of the host star are shown in Figure 4. We found that the host-star mass is ${M}_{\mathrm{host}}={0.40}_{-0.24}^{+0.36}\,{M}_{\odot }$ and the lens system is located at ${D}_{{\rm{L}}}={6.4}_{-1.8}^{+1.3}\,\mathrm{kpc}$ from our solar system. Assuming that the planet frequency of a system like OGLE-2014-BLG-1722Lb,c is uniform in the host-star mass and distance, we use the distribution of the host-star mass and distance to estimate the physical values of the planets. Figure 5 shows the posterior distribution of the mass and (projected and deprojected) separation of the two planets. We found that the masses of planets b and c are ${m}_{{\rm{b}}}={56}_{-33}^{+51}\,$ and ${m}_{{\rm{c}}}={85}_{-51}^{+86}\,{M}_{\oplus }$. The projected separation between the host and planet b is ${r}_{\perp ,{\rm{b}}}=1.5\pm 0.6\,\mathrm{au}$. For planet c, the projected separation from the host is ${r}_{\perp ,{\rm{c}}}={1.7}_{-0.6}^{+0.7}\,$ and ${r}_{\perp ,{\rm{c}}}={2.7}_{-1.0}^{+1.1}\,\mathrm{au}$ for the close and wide models, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Physical Parameters for Lens System

Parameters	Units	Values
		Close	Wide
DL	kpc	${6.40}_{-1.87}^{+1.28}$	${6.38}_{-1.93}^{+1.29}$
${M}_{\mathrm{host}}$	M⊙	${0.40}_{-0.24}^{+0.36}$	${0.39}_{-0.24}^{+0.36}$
${m}_{{\rm{b}}}$	M⊕	${55.3}_{-33.1}^{+51.3}$	${55.0}_{-33.2}^{+51.4}$
${m}_{{\rm{c}}}$	${M}_{\oplus }$	${83.7}_{-51.2}^{+86.4}$	${83.3}_{-51.3}^{+86.6}$
${r}_{\perp ,{\rm{b}}}$	au	1.5 ± 0.6	1.5 ± 0.6
${r}_{\perp ,{\rm{c}}}$	au	${1.7}_{-0.6}^{+0.7}$	${2.7}_{-1.0}^{+1.1}$
${a}_{3{\rm{D}},{\rm{b}}}$	au	${1.9}_{-0.8}^{+1.3}$	${1.9}_{-0.8}^{+1.3}$
${a}_{3{\rm{D}},{\rm{c}}}$	au	${2.5}_{-1.1}^{+2.2}$	${3.5}_{-1.4}^{+2.4}$

Download table as:  ASCIITypeset image

Revealing the semimajor axis, inclination, and eccentricity of the two cold planets can provide stronger constraints on the formation and evolution of this system. The Jupiter–Saturn analog system, OGLE-2006-BLG-109, provided us its orbital information, since the effects of the orbital motion in the lens system were observed (Bennett et al. 2010). For OGLE-2014-BLG-1722, however, such a constraint cannot be obtained because we did not detect the orbital motion effect of the planets due to the short planetary perturbations and no caustic crossings. We assume coplanar and circular orbits to estimate the semimajor axis of the planets. Also, we apply the Hill stability constraint on planet c (Gladman 1993) with a given ${a}_{3{\rm{D}},{\rm{b}}}$. With the random inclination and orbital phase, we estimated that the semimajor axis of planet b is ${a}_{3{\rm{D}},{\rm{b}}}\,={1.9}_{-0.8}^{+1.3}\,\mathrm{au}$. The estimated semimajor axis of planet c is ${a}_{3{\rm{D}},{\rm{c}}}\,={2.5}_{-1.1}^{+2.2}\,$ and ${a}_{3{\rm{D}},{\rm{c}}}={3.5}_{-1.4}^{+2.4}\,\mathrm{au}$ for the close and wide models, respectively. Assuming a snow line of ${a}_{\mathrm{snow}}=2.7{M}_{\mathrm{host}}/{M}_{\odot }\,\mathrm{au}$ (Kennedy & Kenyon 2008), both planets are likely located at 1.4–2.4 times the snow line, which is typical for microlens planets.

Figure 6 shows the expected apparent lens brightness in the I, H, and K bands, as well as the expected relative proper motion, μ_rel, in the geocentric frame. These distributions are derived based on the Bayesian analysis described above with a constraint of the observed t_E and upper limit on π_E. The extinction in the I band to the source star, ${A}_{{\rm{I}},{\rm{S}}}=1.75\pm 0.09$, is estimated by comparing the apparent position of a red clump giant (RCG) on the OGLE-III Color–Magnitude Diagram (CMD), ${(V-I,I)}_{\mathrm{RCG}}=(2.650,16.258)\pm (0.014,0.003)$, with the intrinsic value ${(V-I,I)}_{0}=(1.06,14.51)\pm (0.06,0.09)$ from Bensby et al. (2011) and Nataf et al. (2013). The extinction in the H and K bands to the source was estimated by using the BEAM calculator (Gonzalez et al. 2012), where we take an average for the extinction law of Cardelli et al. (1989) and Nishiyama et al. (2009). Then, following Bennett et al. (2015), the extinction to the foreground lens is given by

$$\begin{eqnarray}&&{A}_{{\rm{\Lambda }},{\rm{L}}}=\displaystyle \frac{1-{e}^{-{D}_{{\rm{L}}}/{h}_{\mathrm{dust}}}}{1-{e}^{-{D}_{{\rm{S}}}/{h}_{\mathrm{dust}}}}{A}_{{\rm{\Lambda }},{\rm{S}}},\end{eqnarray} \tag{ 2 }$$

where the index Λ refers to the passband I, H, or K, and ${h}_{\mathrm{dust}}=(0.10\pm 0.02\,\mathrm{kpc})/\sin \,b$. Unfortunately, the expected lens brightness is much fainter than the source star, even for the K band, and thus it will be difficult to detect the lens flux on the superposition of the relatively bright source. Instead, the color-dependent centroid shift (Bennett et al. 2006, 2007, 2015; Bhattacharya et al. 2017) in the I and K bands, ${\rm{\Delta }}{x}_{I-K}$, would be detectable in 5 yr, as shown in the bottom panel of Figure 6. Bennett et al. (2006) detected ${\rm{\Delta }}{x}_{V-I}$ of 0.6 mas for OGLE-2003-BLG-235 using the Hubble Space Telescope (HST). We expect that the James Webb Space Telescope (JWST) will measure this effect with even better accuracy and put much tighter constraints on the physical parameters of the lens system.

Zoom In Zoom Out Reset image size

The light curve of this event was well observed by the survey teams. Especially, MOA had dense coverage each night while the source star was magnified, which resulted in higher detection efficiencies to planets than the average efficiency of the MOA survey in 2007–2012 (Suzuki et al. 2016). Therefore, an extended study of the statistical analysis of the MOA data will automatically include this event. Here we can roughly (but robustly) estimate the occurrence rate of multiple cold giant planet systems by using the μFUN (Gould et al. 2010) and extended MOA survey data and the two detections of such systems: OGLE-2006-BLG-109Lbc and OGLE-2014-BLG-1722Lbc. We can constrain the occurrence rate of the multiple systems only for cold giant planets from the current microlensing data for the following reasons. First, the most sensitive region of the microlensing method is located at a few times of the snow line. Second, the detection efficiency of two-planet systems with masses of Neptune mass (or below) is too low to put a strong constraint on its occurrence rate (because the sensitivity is approximately the product of detection efficiency for each planet). Third, the detected multiple-planet systems included in this analysis are a Jupiter–Saturn analog and two-Saturn system. A more specific definition of multiple cold giant planet systems will be described later in this section.

To estimate the occurrence rate, we need to consider the detection efficiency (Gould et al. 2010; Cassan et al. 2012; Shvartzvald et al. 2016; Suzuki et al. 2016) of both planets in each observed event. First, we assume that the detection efficiency of two-planet systems can be approximated as ${\epsilon }_{12}\equiv \epsilon ({q}_{1})\times \epsilon ({q}_{2})$. This is safe, especially for the sample dominated by low-magnification events (Gaudi et al. 1998; Han & Park 2002) like the MOA survey data. Second, we assume that the MOA's survey sensitivity (the total detection efficiency) in 2013–2016 is similar to that of 2007–2012. We randomly pick up 4 yr from 2007 to 2012 and use the detection efficiency in each event as the hypothetical MOA data of 2013–2016. (The full analysis of the 2013–2016 data is beyond the scope of this paper.) Then, one event randomly selected in the hypothetical 2014 data is substituted by the real detection efficiency of OGLE-2014-BLG-1722 (whose MOA name is MOA-2014-BLG-490). We repeated this process and confirmed that the occurrence rate does not depend on the data used for the hypothetical survey data. For the μFUN sample, we use the detection efficiencies estimated by Gould et al. (2010). For the MOA sample, we use the detection efficiencies in Suzuki et al. (2016), which are determined by using the method of Rhie et al. (2000) with an assumption of logarithmically uniform planet distribution in both mass ratio and separation. Thus, the total sample we use here consists of 13 events from μFUN (2005–2008) and about 2400 events from MOA (2007–2016). The overlapping events are removed from the MOA sample, and we use the μFUN detection efficiencies for them.

The occurrence rate of multiple cold giant planet systems, η_multi, can be estimated by maximizing the likelihood,

$$\begin{eqnarray}&&L({\eta }_{\mathrm{multi}})\equiv P({\eta }_{\mathrm{multi}}| d)\propto P({\eta }_{\mathrm{multi}})\displaystyle \prod _{i}^{N}{f}_{i}\displaystyle \prod _{j}^{M-N}(1-{f}_{j}),\end{eqnarray} \tag{ 3 }$$

where $P({\eta }_{\mathrm{multi}}| d)$ shows the probability of η_multi conditioned on the data d; N and M are the number of detected multiple systems and all the events included, respectively; and f is the product of the detection efficiency and the occurrence rate, $f={\epsilon }_{12}\,\times {\eta }_{\mathrm{multi}}$. We assume that $P({\eta }_{\mathrm{multi}})$ is a uniform distribution; thus, the likelihood function is given by the binomial distribution. For each event, we integrate the detection efficiency within $-0.5\lt \mathrm{log}\,s\lt 0.5$ where the sensitivity is highest and use the mass ratio range of $5\times {10}^{-5}\lt q\lt 0.03$ to compute ₁₂. As indicated in Figure 7, we use a box that has 0.5 dex in each side on the q₁–q₂ plane and calculate the average detection efficiency in this box. The occurrence rate at the center of the box is calculated by using the averaged detection efficiency. Then, by moving the box on this plane, we calculate the occurrence rate of the system with any combination of q₁ and q₂, where q₂ > q₁. Thus, η_multi is the occurrence rate for the planetary systems whose host stars are orbited by at least two planets, both of which are located in the separation of $-0.5\lt \mathrm{log}\,s\lt 0.5$ and in a 0.5 dex bin centered on $\mathrm{log}\,{q}_{j}$, where j is 1 or 2.

Zoom In Zoom Out Reset image size

We computed the upper limit on the occurrence rate at the 90% confidence level, which is indicated in Figure 7. As we expect from the detection efficiency, the upper right corner of the figure has a lower occurrence rate limit, and the bottom left corner has a higher occurrence rate limit. We averaged over the upper limit in the area where the box can include the two detections and found η_multi < 0.25 with a 90% confidence level. Furthermore, we estimate the mode, median, and 1σ uncertainty of the occurrence rate in the same parameter space. As indicated in Figure 8, we found ${\eta }_{\mathrm{multi}}={0.12}_{-0.06}^{+0.09}$, with a mode of 0.09. This result is not sensitive to the size of the sliding box. Also note that the mass ratios of Jupiter and Saturn in our solar system happen to be included in the area where we calculate the η_multi.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The most likely value of the occurrence rate, η_multi = 0.09, can be compared with the frequency of the solar-like systems estimated by Gould et al. (2010), 0.17, which is consistent with η_multi due to the large uncertainty but a factor of 1.9 larger than our mode value. They considered a scaled solar system including Uranus and Neptune²⁷ analogs, so the definition of the solar-like system is different from that of the multiple cold gas giant planet system. In their simulation, however, most of the detected multiple systems are Jupiter–Saturn analogs. Also, their estimated number (0.17) would decrease by about a factor of two if we consider an extended analysis of Gould et al. (2010). They found six planets with one multiple system in their 4 yr of μFUN data, but we expect that only three planets (Bachelet et al. 2012; Yee et al. 2012; Fukui et al. 2015) and no multiple system would be included in the next 6 yr. Note that OGLE-2012-BLG-0026 (Han et al. 2013) would not be included in the simple extended study, because its peak magnification does not satisfy the u₀ < 0.005 criterion in Gould et al. (2010). Therefore, with the assumption of a linear extrapolation of the 4 yr sensitivity to 10 yr, the frequency of solar-like systems would be roughly 0.07, which is closer to our mode value.

Considering that the average number of cold planets with a mass ratio of q = 4.5 × 10⁻⁴ per star is 0.14, which is estimated by integrating the best-fit mass ratio function of Suzuki et al. (2016) within the separation of $-0.5\lt \mathrm{log}\,s\lt +0.5$ and mass ratio of $-3.6\lt \mathrm{log}\,q\lt -3.1$, the estimated mode of η_multi could be overestimated. If we assume that all planetary systems with cold gas giants are two-planet systems, the occurrence rate of two-planet systems will be 0.07 at maximum. Thus, we can use the average number of planets for a prior probability of the upper limit on η_multi. We use the uncertainty in the mass ratio function in Suzuki et al. (2016) to estimate the upper limit distribution on η_multi, whose median and 1σ uncertainty are 0.058 and 0.008, respectively. Here we assume that the upper limit on η_multi of the systems with ${q}_{1}={q}_{2}=4.5\times {10}^{-4}$ is more conservative than that of the systems with ${q}_{1}=4.5\,\times {10}^{-4}$ and q₂ = 3.0 × 10⁻³. Adopting a 3σ uncertainty for the prior distribution, we estimate η_multi = 0.06 ± 0.02, with a mode of 0.07, as shown in Figure 9. The estimated occurrence rate is strongly dependent on the assumed prior distribution of η_multi. Nevertheless, the constrained η_multi is more consistent with the hypothetically extended study of Gould et al. (2010) in the previous paragraph. Although there are still large uncertainties in η_multi, future studies of the occurrence rate of multiple cold planet systems are important to understand planet formation beyond the snow line.

The microlensing event OGLE-2014-BLG-1722 showed two relatively short deviations from the expected single-lens light curve. The triple-lens model consisting of a primary lens star and two planetary mass ratio objects is the best-fit model to explain the two anomalies. Thus, if the two-planet model is true, this is the first detection of a two-planet system in low-magnification microlensing events and the third detection of a two-planet system found by the microlensing method. Also, this event is the first two-planet system detected in pure survey data. We found that the first anomaly was caused by planet b with a mass ratio of $q=({4.5}_{-0.6}^{+0.7})\times {10}^{-4}$ and separation of s = 0.753 ± 0.004. The second anomaly revealed planet c by the χ² improvement of Δχ² ∼ 170 compared to the one-planet model, but its separation has degenerate close–wide solutions. For the close solution, the mass ratio and separation for planet c are ${q}_{2}=({7.0}_{-1.7}^{+2.3})\times {10}^{-4}$ and s₂ = 0.84 ± 0.03, respectively. For the wide solution, they are ${q}_{2}=({7.2}_{-1.7}^{+1.8})\times {10}^{-4}$ and s₂ = 1.37 ± 0.04.

The physical properties of the lens system could not be measured because neither finite-source nor microlensing parallax effects were measured. Instead, we estimated the mass of and distance to the lens system with a Bayesian analysis using a standard Galactic model assuming that the planet-hosting probability does not depend on host-star mass and distance. The estimated host-star mass and distance are ${M}_{\mathrm{host}}={0.40}_{-0.24}^{+0.36}\,{M}_{\odot }$ and ${D}_{{\rm{L}}}={6.4}_{-1.8}^{+1.3}\,\mathrm{kpc}$, respectively. We estimate that the masses of planets b and c are ${m}_{{\rm{b}}}={56}_{-33}^{+51}\,$ and ${m}_{{\rm{c}}}={85}_{-51}^{+86}\,{M}_{\oplus }$. The projected distance between the host star and planet b is ${r}_{\perp ,{\rm{b}}}=1.5\pm 0.6\,\mathrm{au}$. For the close and wide models of planet c, the projected separations are ${r}_{\perp ,{\rm{c}}}={1.7}_{-0.6}^{+0.7}\,$ and ${r}_{\perp ,{\rm{c}}}={2.7}_{-1.0}^{+1.1}\,\mathrm{au}$, respectively. We also estimated the semimajor axis of the planets assuming circular and coplanar orbits. The semimajor axis of planet b is ${a}_{3{\rm{D}},{\rm{b}}}={1.9}_{-0.8}^{+1.3}\,\mathrm{au}$, and that of planet c is ${a}_{3{\rm{D}},{\rm{c}}}={2.5}_{-1.1}^{+2.2}\,$ and ${a}_{3{\rm{D}},{\rm{c}}}={3.5}_{-1.4}^{+2.4}\,\mathrm{au}$ for the close and wide models, respectively. Therefore, OGLE-2014-BLG-1722L is a late-type star orbited by two Saturn-mass planets that are located at a separation of a few times the snow line.

Although the two-planet scenario is preferred by 3.1σ, one cannot rule out the interpretation where the source is a binary. For OGLE-2005-BLG-390 (Beaulieu et al. 2006) and OGLE-2016-BLG-1195 (Bond et al. 2017), a binary-source model was ruled out by Δχ² = 46 and 120, respectively. Nevertheless, the marginal detection of planet c is important for future statistical studies. The statistical study of 6 yr MOA data (Suzuki et al. 2016) dealt with an ambiguous event where the planetary and stellar binary models comparably well explain the observed light curves.

Generally, triple-lens events are challenging to find the best model because of not only the vast parameter space we need to explore but also possible degenerate models. The lens system of a circumbinary planet event, OGLE-2007-BLG-349 (Bennett et al. 2016), was initially thought to be a two-planet system. Also, the one-planet event with a binary source, MOA-2010-BLG-117 (Bennett et al. 2018), initially mimicked a circumbinary planet event. Sometimes, the orbital motion of a lens system could produce an anomaly in the light curve, which apparently seems to be a signal from the third body (Bennett et al. 1999; Albrow et al. 2000; Udalski et al. 2015a; Han et al. 2016). In these cases, the possible degenerate solutions were finally ruled out, but there is one triple-lens event waiting for modeling.

Although the detection efficiency of multiple planets is very high in high-magnification events (Gaudi et al. 1998), the separations could suffer from the close–wide degeneracy. (OGLE-2006-BLG-109Lbc (Gaudi et al. 2008) was detected through the high-magnification channel, but the degenerate models were ruled out nonetheless.) On the other hand, low-magnification events are generally expected to have less-degenerate models and unique solutions. However, we need to overcome the low detection efficiency by continuous light-curve coverage.

The two two-planet systems previously detected by microlensing, OGLE-2006-BLG-109L and OGLE-2012-BLG-0026L, were well characterized because both finite-source and parallax effects were measured, and the follow-up high-resolution images put further constraints on the lens physical properties. For most of the planetary events, including two-planet systems, however, we expect low magnification with a moderate event timescale like that of OGLE-2014-BLG-1722, where the mass measurement is difficult from ground-based telescopes alone. Future space-based microlensing surveys, such as WFIRST (Spergel et al. 2015) and hopefully Euclid (Penny et al. 2013) as well, are expected to measure the physical parameters of the lens systems for most of the observed microlensing events, including two-planet events, which will eventually allow us to study the occurrence rate of multiple cold planet systems as a function of host mass and galactocentric distance.

The occurrence rates of multiple cold giant systems η_multi were estimated in Section 5 with the two detections of such systems (OGLE-2006-BLG-109Lbc and OGLE-2014-BLG-1722Lbc) in the sample of the 4 yr μFUN data (Gould et al. 2010) and 10 yr MOA survey data, which consist of the 6 yr full analysis of detection efficiency (Suzuki et al. 2016) and the simulated hypothetical 4 yr data. We found ${\eta }_{\mathrm{multi}}\,={0.12}_{-0.06}^{+0.09}$ with a mode value of 0.09. Thus, roughly 10% of stars, which specifically could become lens stars of microlensing events toward the Galactic bulge, are orbited by two cold giant planets. To confirm this result, we need larger statistics of multiple-planet system detections. The occurrence rate of the 90% confidence level upper limit (Figure 7) shows lower and higher occurrence rates for massive and less massive planet pairs. This can be inferred because massive planets easily affect the orbits of other planets, and they might let the other planets migrate in/outward. Also, we know that there are many tightly packed coplanar systems consisting of super-Earth or Earth-mass planets found by Kepler. Thus, we would speculate that multiple cold Neptunes and/or super-Earth systems could be common. Such systems can be found only by microlensing surveys. The current ground-based microlensing surveys do have sensitivities to such systems, but the sensitivity is just very low. The 24 hr continuous light-curve coverage by OGLE, MOA, and KMTNet (Kim et al. 2016) will increase the number of detections of multiple cold planet systems. The ultimate light-curve coverage by WFIRST in terms of continuous coverage and high precision photometry is expected to dramatically increase the detection of multiple systems (Bennett & Rhie 2002), as its sensitivity is much higher for both the mass ratio and separation parameter space. The precisely measured occurrence rate of multiple systems can be directly compared with the population synthesis simulations (e.g., Ida & Lin 2004; Mordasini et al. 2009). Statistical studies with the larger sample will allow us to put constraints on the planet formation models from the point of view of multiple cold planets.

DS and DPB acknowledge support from NASA grants NNX13AF64G, NNX14AG49G, and NNX15AJ76G. TS acknowledges financial support from JSPS23103002, JSPS24253004, and JSPS26247023. AY acknowledges financial support from JSPS25870893. The MOA project is supported by grants JSPS25103508, JSPS23340064, and JP17H02871. The OGLE Team thanks Profs. M Kubiak and G Pietrzyński for their contribution to the collection of the OGLE photometric data. The OGLE project has received funding from the National Science Centre, Poland, grant MAESTRO 2014/14/A/ST9/00121 to AU. Work by C. Han was supported by a grant (2017R1A4A1015178) from the National Research Foundation of Korea.