Toward a Galactic Distribution of Planets. I. Methodology and Planet Sensitivities of the 2015 High-cadence Spitzer Microlens Sample

All events in our sample were found by the Optical Gravitational Lensing Experiment (OGLE) collaboration in real time through its Early Warning System (Udalski et al. 1994; Udalski 2003), based on observations with the 1.4 deg² camera on its 1.3 m Warsaw Telescope at the Las Campanas Observatory in Chile (Udalski 2003; Udalski et al. 2015a). These events received OGLE-IV observations with cadences varying from 3 to 30 per day. The coordinates, OGLE-IV fields, and cadences of individual events are provided in Table 1.

OGLE data were reduced using the photometry software developed by Wozniak (2000) and Udalski (2003), which was based on the Difference Image Analysis (DIA) technique (Alard & Lupton 1998).

KMTNet consists of three 1.6 m telescopes located at CTIO in Chile, SAAO in South Africa, and SSO in Australia. Observations were initiated on February 3 (JD = 2,457,056.9), February 19 (JD = 2,457,072.6), and June 9 (JD = 2,457,182.9) in 2015 from CTIO, SAAO, and SSO, respectively. Each telescope is equipped with a 4 deg² field-of-view camera and observes the ∼16 deg² prime microlensing fields at a ∼10-minute cadence when the bulge is visible.

The KMTNet data were reduced by the DIA photometric pipeline (Alard & Lupton 1998; Albrow et al. 2009).

As detailed in Yee et al. (2015b), the Spitzer program is designed to maximize the sum of the products ${\sum }_{i}{S}_{i}{P}_{i}$, where S_i is the planet sensitivity of event i and P_i is the probability to measure the microlens parallax of this event. As a consequence, the Spitzer team started selecting targets beginning in early 2015 May, although Spitzer did not start taking data until JD' = JD –2,450,000 = 7180.2 (2015 June 6.7). To enforce our blindness to the existence of planets in any events, we select events if (1) they meet certain objective criteria at the time of one of the uploads of targets to Spitzer, in which case they are considered as "objectively chosen," or (2) they do not meet objective criteria but are nevertheless selected on the basis that the Spitzer team believes that by selecting them the quantity ${\sum }_{i}{S}_{i}{P}_{i}$ can be maximized. Events selected in the latter case are known as "subjectively chosen." For objectively chosen events, planets and planet sensitivities from before or after the Spitzer selection dates can be incorporated into the statistical analysis, while for subjectively chosen events, only planets (and planet sensitivities) from after the Spitzer selection dates can be included in the final sample.²⁰ One relevant point is that any event that is originally subjectively chosen but later meets objective criteria will be considered as objectively chosen (provided that its parallax is measurable based on the restricted set of Spitzer data acquired after the date it became objective).

Events, once selected, are given Spitzer cadences according to suggestions in Yee et al. (2015b). The majority of events received Spitzer observations at a 1-day cadence. Higher cadences were assigned to a few events, if the Spitzer team believed that the nominal cadence would lead to failures in parallax measurements. After all targets were scheduled according to their adopted cadences, the remaining time, if any, was applied to events that appeared or would appear with relatively high magnification as seen from the ground. Spitzer observations stopped if the predefined criteria for stopping observations in Yee et al. (2015b) were met, or the event exited the Spitzer Sun-angle window. Our last Spitzer observation was taken on JD' = 7222.28. In Table 1 we provide the information for Spitzer selection and observation of each individual event.

Spitzer data were reduced using the customized software that was developed by Calchi Novati et al. (2015) specifically for this program. This software improved the performance of Spitzer IRAC photometry in crowded fields, although unknown systematics may persist in some cases. We discuss this in Section 5.1.

The characterization of a microlensing event requires a measurement of the color of the source star. This is usually achieved by using the less frequent V-band observations from survey teams, but it does not work for events that are highly extincted in optical bands. For this reason, we also obtained observations of all Spitzer targets using the ANDICAM (DePoy et al. 2003) dichroic camera on the 1.3 m SMARTS telescope at CTIO. These observations were made simultaneously in I and H bands and were for the specific purpose of inferring the $I-[3.6\mu {\rm{m}}]$ color of the source star. These additional color data were reduced using DoPhot (Schechter et al. 1993).

The separation between Earth and the satellite perpendicular to the line of sight to the microlensing event, ${D}_{\perp }$, causes apparent changes in the angular lens–source separation ${\rm{\Delta }}\theta ={\pi }_{\mathrm{rel}}{D}_{\perp }/\mathrm{au}$, and this in turn gives rise to different microlensing light curves. These light curves, as seen from Earth and from the satellite, appear to peak at different times t₀ and with different impact parameters u₀ (normalized to ${\theta }_{{\rm{E}}}$). In the approximation of rectilinear motion of Earth and the satellite (Refsdal 1966; Gould 1994; Graff & Gould 2002)

$$\begin{eqnarray}&&{{\boldsymbol{\pi }}}_{{\rm{E}}}\approx \displaystyle \frac{\mathrm{au}}{{D}_{\perp }}({\rm{\Delta }}\tau ,{\rm{\Delta }}\beta ),\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{\rm{\Delta }}\tau \equiv \displaystyle \frac{{t}_{0,\mathrm{sat}}-{t}_{0,\oplus }}{{t}_{{\rm{E}}}};\quad {\rm{\Delta }}\beta \equiv {u}_{0,\mathrm{sat}}-{u}_{0,\oplus }.\end{eqnarray} \tag{ 7 }$$

Unfortunately, u₀ is a signed quantity (depending on whether the lens passes the source on its right or left; for sign definition see Figure 4 of Gould 2004), while only $| {u}_{0}|$ is directly measurable from the light curve. Therefore, satellite parallax measurements are subject to a fourfold degeneracy²¹

$$\begin{eqnarray}&&{{\boldsymbol{\pi }}}_{{\rm{E}}}\approx \displaystyle \frac{\mathrm{au}}{{D}_{\perp }}({\rm{\Delta }}\tau ,{\rm{\Delta }}{\beta }_{\pm ,\pm }),\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}\left\{\begin{array}{ll}{\rm{\Delta }}{\beta }_{+,+}\equiv +{u}_{0,\mathrm{sat}}-{u}_{0,\oplus }, & \,(+,+)\mathrm{solution}\\ {\rm{\Delta }}{\beta }_{+,-}\equiv -{u}_{0,\mathrm{sat}}-{u}_{0,\oplus }, & \,(+,-)\mathrm{solution}\\ {\rm{\Delta }}{\beta }_{-,-}\equiv -{u}_{0,\mathrm{sat}}+{u}_{0,\oplus }, & \,(-,-)\mathrm{solution}\\ {\rm{\Delta }}{\beta }_{-,+}\equiv +{u}_{0,\mathrm{sat}}-{u}_{0,\oplus }, & \,(-,+)\mathrm{solution}\end{array}\right..\end{eqnarray} \tag{ 9 }$$

In principle, higher-order effects in the light curve itself can break this degeneracy. At first order (in the polynomial expansion of Smith et al. 2003), it can be broken from the different Einstein timescales ${t}_{{\rm{E}}}$ measured from Earth and the satellite owing to their relative motion (even within the approximation of rectilinear motion) (Gould 1995). At third and fourth order, it can be broken owing to parallax effects from the accelerated motion of Earth (Gould 1992). In practice, however, these effects are usually quite weak. First, with current experiments, ${t}_{{\rm{E}}}$ is normally not independently measured from the satellite simply because the observational investment for this would be extremely high (Gaudi & Gould 1997a), and these resources are better applied to observing more events. Ground-based parallaxes are rarely measured because the Einstein timescales are typically small, ${t}_{{\rm{E}}}\lt \mathrm{yr}/2\pi$, so that third- and particularly fourth-order effects are very subtle. This indeed is the reason for going to space. Nevertheless, although these higher-order effects are small, they can contribute to breaking the degeneracy between well-determined but otherwise indistinguishable parallax solutions.

We search for and characterize the four solutions using the code developed in Zhu et al. (2016a). We first find a simple three-parameter $({t}_{0},{u}_{0},{t}_{{\rm{E}}})$ solution based on OGLE data. Next, we include Spitzer data, introduce two parameters ${\pi }_{{\rm{E}}}n$ and ${\pi }_{{\rm{E}}}e$, which are the two components of vector ${{\boldsymbol{\pi }}}_{{\rm{E}}}$ along the north and east directions, respectively, and easily find one of the four parallax solutions by allowing ${\chi }^{2}$ to go downhill. As per the usual convention, these parameters (${\pi }_{{\rm{E}}}n$, ${\pi }_{{\rm{E}}}e$, ${t}_{{\rm{E}}}$) are defined in the geocentric frame (Gould 2004).²² The location of the Spitzer satellite is extracted from the JPL Horizons website,²³ enabling a self-consistent quantification of the microlens parallax effect and the event timescale. In addition, there are two flux parameters for each data set, ${F}_{{\rm{s}}}$ and ${F}_{{\rm{b}}}$. The former is the flux from the source, and the latter is the flux that is blended within the aperture and does not participate in the event. The model for the total flux observed at epoch t_i for data set j is then given by

$$\begin{eqnarray}&&{F}^{j}({t}_{i})={F}_{{\rm{s}}}^{j}\cdot {A}^{j}({t}_{i};{t}_{0},{u}_{0},{t}_{{\rm{E}}},\rho ,{{\boldsymbol{\pi }}}_{{\rm{E}}})+{F}_{{\rm{b}}}^{j}.\end{eqnarray} \tag{ 10 }$$

Once a solution is found, we estimate the uncertainties of parameters via a Markov Chain Monte Carlo (MCMC) analysis, using the emcee ensemble sampler (Foreman-Mackey et al. 2013). The remaining three solutions are also easily found by seeding solutions at the locations expected based on Equation (8). In some cases, typically events with long timescales or events peaking near the beginning of the season, there is no local minimum at ${\chi }^{2}$ surface for one or more solutions owing to strong parallax information from the ground. Within the mathematical formalism that follows, these other solutions can be thought of as "existing" but having very high ${\rm{\Delta }}{\chi }^{2}$ relative to the best solution.

For each of the four solutions we then derive ${{\boldsymbol{\pi }}}_{{\rm{E}}}$, ${\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}$, the uncertainty of the latter quantity, and ${\rm{\Delta }}{\chi }^{2}$ relative to the best solution. Here ${\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}$ is the transverse velocity between the source and the lens projected onto the observer plane, after the correction from geocentric to heliocentric frames,

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}={\tilde{{\boldsymbol{v}}}}_{\mathrm{geo}}+{{\boldsymbol{v}}}_{\oplus ,\perp };\quad {\tilde{{\boldsymbol{v}}}}_{\mathrm{geo}}=\displaystyle \frac{\mathrm{au}}{{t}_{{\rm{E}}}}\displaystyle \frac{{{\boldsymbol{\pi }}}_{{\rm{E}}}}{{\pi }_{{\rm{E}}}^{2}},\end{eqnarray} \tag{ 11 }$$

where ${{\boldsymbol{v}}}_{\oplus ,\perp }$ is the velocity of Earth at the event peak and projected perpendicular to the directory of the event. To facilitate further discussions, we also define here the event timescale in the heliocentric frame

$$\begin{eqnarray}&&{t}_{{\rm{E}}}^{{\prime} }\equiv \displaystyle \frac{{\tilde{r}}_{{\rm{E}}}}{{\tilde{v}}_{\mathrm{hel}}};\quad {\tilde{r}}_{{\rm{E}}}\equiv \displaystyle \frac{\mathrm{au}}{{\pi }_{{\rm{E}}}}.\end{eqnarray} \tag{ 12 }$$

In deriving event solutions, we are able to incorporate "color constraints" (either ${VI}[3.6\,\mu {\rm{m}}]$ or ${IH}[3.6\,\mu {\rm{m}}]$) into the fit. This is either very important or essential for the great majority of cases, as anticipated by Yee et al. (2015b). The naive idea of space-based parallaxes, as outlined by Refsdal (1966) and Gould (1994) and as captured by Equation (6), is that t₀ and u₀ will be measured independently from the satellite and Earth. However, such independent measurements are essentially impossible if the event is not observed over (or at least close to) peak. Hence, in the 2014 pilot program, exceptional efforts were made to observe over peak, which greatly restricted the number of events that could be targeted, given the short (∼38-day) observing window set by Spitzer Sun-angle restrictions and given the delays of 6 ± 3 days in observing targets (Figure 1 of Udalski et al. 2015b). However, based on experience in optical bands (Yee et al. 2012), Yee et al. (2015b) argued that, even if the peak were not observed from the satellite, it would be possible to recover ${({t}_{0},{t}_{{\rm{E}}})}_{\mathrm{sat}}$ provided that the Spitzer source flux could be determined from a combination of (1) the measured source flux of the ground-based light curve, (2) the measured source color in ground-based bands (V − I or I − H), and (3) a color–color relation (e.g., ${VI}[3.6\,\mu {\rm{m}}]$) derived from field stars. In practice, we derive the $(I-[3.6\,\mu {\rm{m}}])$ from the measured color and color–color relation and then impose the 2σ limits of this measurement as hard constraints in the fit.

4.2.1. Stellar Density Profile

The Galactic center has equatorial coordinates $({\alpha }_{\mathrm{GC}},{\delta }_{\mathrm{GC}})\,=({17}^{{\rm{h}}}{45}^{{\rm{m}}}37\buildrel{\rm{s}}\over{.} 224,-28^\circ 56^{\prime} 10\buildrel{\prime\prime}\over{.} 23)$ (Reid & Brunthaler 2004) and heliocentric distance ${R}_{\mathrm{GC}}=8.3\,\mathrm{kpc}$ (Gillessen et al. 2009). The Sun is above the Galactic midplane $(z=0)$ by $27\,$ pc (Chen et al. 2001), which corresponds to a tilt angle $\beta =0\buildrel{\circ}\over{.} 19$.

The total stellar number density ${n}_{\star }$ at given galactocentric coordinates $(x,y,z)$ is the sum of contributions from the bulge and disk components

$$\begin{eqnarray}&&{n}_{\star }(x,y,z)={n}_{{\rm{B}}}({x}^{{\prime} },{y}^{{\prime} },{z}^{{\prime} })+{n}_{{\rm{D}}}(R,z).\end{eqnarray} \tag{ 13 }$$

We assume a triaxial G2 model for the bulge component (Kent et al. 1991; Dwek et al. 1995):

$$\begin{eqnarray}&&{n}_{{\rm{B}}}={n}_{{\rm{B}},0}{e}^{-{r}_{s}^{2}/2};\quad {r}_{s}\equiv {\left\{{\left[{\left(\displaystyle \frac{x^{\prime} }{{x}_{0}}\right)}^{2}+{\left(\displaystyle \frac{y^{\prime} }{{y}_{0}}\right)}^{2}\right]}^{2}+{\left(\displaystyle \frac{z^{\prime} }{{z}_{0}}\right)}^{4}\right\}}^{1/4},\end{eqnarray} \tag{ 14 }$$

where ${n}_{{\rm{B}},0}=13.7$ pc⁻³, ${x}_{0}=1.59\,$ kpc, ${y}_{0}=424\,$ pc, and ${z}_{0}=424\,$ pc. These values are adopted from Robin et al. (2003). The coordinates $(x^{\prime} ,y^{\prime} ,z^{\prime} )$ are derived by rotating the galactocentric coordinates $(x,y,z)$ around the z-axis by ${\alpha }_{\mathrm{bar}}=30^\circ$ (e.g., Cao et al. 2013; Wegg & Gerhard 2013). The disk component in Equation (13) has the form (Bahcall 1986)

$$\begin{eqnarray}&&{n}_{{\rm{D}}}={n}_{{\rm{D}},0}\exp \left[-\left(\displaystyle \frac{R-{R}_{\mathrm{GC}}}{{R}_{0}}+\displaystyle \frac{| z| }{{z}_{{\rm{D}},0}}\right)\right].\end{eqnarray} \tag{ 15 }$$

Here $R\equiv \sqrt{{x}^{2}+{y}^{2}}$, the local stellar number density ${n}_{{\rm{D}},0}=0.14$ pc⁻³, the scale length of the disk ${R}_{0}=3.5\,\mathrm{kpc}$, and the scale height of the disk ${z}_{{\rm{D}},0}=325\,$ pc (Han & Gould 1995). We show in Figure 1 the stellar number density profile toward the Baade's window, which is approximately the center of microlensing fields.

Zoom In Zoom Out Reset image size

4.2.2. Stellar Velocity Distribution

The mean stellar velocity at galactocentric coordinates $(x,y,z)$ has the form

$$\begin{eqnarray}&&{{\boldsymbol{\mu }}}_{{\boldsymbol{v}}}(x,y,z)=\displaystyle \frac{{n}_{{\rm{B}}}}{{n}_{\star }}{{\boldsymbol{\mu }}}_{v,{\rm{B}}}+\displaystyle \frac{{n}_{{\rm{D}}}}{{n}_{\star }}{{\boldsymbol{\mu }}}_{v,{\rm{D}}},\end{eqnarray} \tag{ 16 }$$

and the velocity dispersion is given by

$$\begin{eqnarray}&&{\sigma }_{v,i}^{2}(x,y,z)={\left(\displaystyle \frac{{n}_{{\rm{B}}}}{{n}_{\star }}\right)}^{2}{\sigma }_{v,i,{\rm{B}}}^{2}+{\left(\displaystyle \frac{{n}_{{\rm{D}}}}{{n}_{\star }}\right)}^{2}{\sigma }_{v,i,{\rm{D}}}^{2}\quad (i=x,y,z).\end{eqnarray} \tag{ 17 }$$

We assume that the bulge stars have zero mean velocity and $120\,\mathrm{km}\,{{\rm{s}}}^{-1}$ velocity dispersion along each direction (${\sigma }_{v,i,{\rm{B}}}=120\,\mathrm{km}\,{{\rm{s}}}^{-1}$). The latter is derived from the proper-motion dispersion of bulge stars ${\sigma }_{\mu }=3\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$ (Poleski et al. 2013). Disk stars partake of the flat rotation curve with $240\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (i.e., ${\mu }_{v,z,{\rm{D}}}=0\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and ${\mu }_{v,y,{\rm{D}}}=240\,\mathrm{km}\,{{\rm{s}}}^{-1}$; Reid et al. 2014), and their velocity dispersions are 18 and $33\,\mathrm{km}\,{{\rm{s}}}^{-1}$ in the vertical $(z)$ and rotation $(y)$ directions, respectively. The Sun partakes of the same rotation curve and has a peculiar motion (${V}_{\odot }=12\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and ${W}_{\odot }=7\,\mathrm{km}\,{{\rm{s}}}^{-1}$; Schönrich et al. 2010) relative to the local standard of rest.

4.2.3. Stellar Mass Function

We choose two forms of the lens mass function (MF): (1) a flat MF with $d\xi ({M}_{{\rm{L}}})/d\,\mathrm{log}\,{M}_{{\rm{L}}}\propto 1$ and (2) a Kroupa MF (Kroupa 2001)

$$\begin{eqnarray}\displaystyle \frac{d\xi ({M}_{{\rm{L}}})}{d\,\mathrm{log}\,{M}_{{\rm{L}}}}\propto \left\{\begin{array}{ll}{M}_{{\rm{L}}}^{0.7}, & \,0.013\lt {M}_{{\rm{L}}}/{M}_{\odot }\lt 0.08\\ {M}_{{\rm{L}}}^{-0.3}, & \,0.08\lt {M}_{{\rm{L}}}/{M}_{\odot }\lt 0.5\\ {M}_{{\rm{L}}}^{-1.3}, & \,0.5\lt {M}_{{\rm{L}}}/{M}_{\odot }\lt 1.3\end{array}\right..\end{eqnarray} \tag{ 18 }$$

In both cases, no planetary lenses are included, and the upper end of the MF is truncated at $1.3\,{M}_{\odot }$. As has been demonstrated in Calchi Novati et al. (2015) and will also be shown later, the choice of a different MF has essentially no effect on the result.

Following Calchi Novati et al. (2015), we define a lens "distance" parameter ${D}_{8.3}$ that is a monotonic function of ${\pi }_{\mathrm{rel}}$,

$$\begin{eqnarray}&&{D}_{8.3}\equiv \displaystyle \frac{\mathrm{kpc}}{({\pi }_{\mathrm{rel}}/\mathrm{mas})+1/8.3}.\end{eqnarray} \tag{ 19 }$$

This has the advantage that ${\pi }_{\mathrm{rel}}$ is much better constrained than the lens distance ${D}_{{\rm{L}}}$ (see Equation (5)), and it also informs us more of the Galactic population from which the lens is drawn. That is,

$$\begin{eqnarray}&&{D}_{8.3}\to {D}_{{\rm{L}}}\quad ({D}_{{\rm{L}}}\ll {D}_{{\rm{S}}}),\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&(8.3\,\mathrm{kpc}-{D}_{8.3})\to {D}_{\mathrm{LS}}\quad ({D}_{\mathrm{LS}}\ll {D}_{{\rm{S}}}),\end{eqnarray} \tag{ 21 }$$

where ${D}_{\mathrm{LS}}\equiv {D}_{{\rm{S}}}-{D}_{{\rm{L}}}$. A determination that ${D}_{\mathrm{LS}}\ll {D}_{{\rm{S}}}$ is a much better indicator that the lens is in the bulge than the value of ${D}_{{\rm{L}}}$ (which in any case is less precisely known).

As discussed in Section 1, the lens distance parameter ${D}_{8.3}$ cannot be uniquely determined for the majority of events because of the lack of ${\theta }_{{\rm{E}}}$ measurement. We therefore derive the Bayesian distribution of ${D}_{8.3}$ by imposing a Galactic model. As first pointed out by Han & Gould (1995), such a distribution of ${D}_{8.3}$ is fairly compact if ${{\boldsymbol{\pi }}}_{{\rm{E}}}$ rather than ${\theta }_{{\rm{E}}}$ (which gives ${\mu }_{\mathrm{rel}}$) can be measured. One can understand this by first approximating the Galactic disk lenses as moving exactly on a flat rotation curve and bulge sources as not moving. Then (also approximating the Sun as being at the local standard of rest),

$$\begin{eqnarray}&&{\pi }_{\mathrm{rel}}\to \displaystyle \frac{{\tilde{v}}_{\mathrm{hel}}}{{\mu }_{\mathrm{sgrA}\star }},\end{eqnarray} \tag{ 22 }$$

where ${\mu }_{\mathrm{sgrA}\star }$ is the observed proper motion of the Galactic center and ${\tilde{v}}_{\mathrm{hel}}$ is the magnitude of ${\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}$ from Equation (11). In fact, the velocities of both the sources and lenses are dispersed relative to this naive model. However, since these dispersions (projected on the observer plane) are typically small compared to the projection of the flat rotation curve, the probability distribution of ${\pi }_{\mathrm{rel}}$ (and therefore ${D}_{8.3}$) is typically compact. Then, since ${M}_{{\rm{L}}}={\pi }_{\mathrm{rel}}/(\kappa {\pi }_{{\rm{E}}}^{2})$, ${M}_{{\rm{L}}}$ is also quite well measured.

To further illustrate this point under our adopted Galactic model, we show in Figure 2 the probability distributions of ${\tilde{v}}_{\mathrm{hel}}$ and ${\mu }_{\mathrm{rel}}$ for several different lens distances. Here ${\tilde{v}}_{\mathrm{hel}}$ and ${\mu }_{\mathrm{rel}}$ are the amplitudes of the two vectors ${\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}$ and ${{\boldsymbol{\mu }}}_{\mathrm{rel}}$, respectively, and these vectors are related to lens and source properties by

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}=\displaystyle \frac{{D}_{{\rm{S}}}}{{D}_{\mathrm{LS}}}{{\boldsymbol{v}}}_{{\rm{L}},\mathrm{gc}}-\displaystyle \frac{{D}_{{\rm{L}}}}{{D}_{\mathrm{LS}}}{{\boldsymbol{v}}}_{{\rm{S}},\mathrm{gc}}-{{\boldsymbol{v}}}_{\odot ,\mathrm{gc}},\end{eqnarray} \tag{ 23 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{\mu }}}_{\mathrm{rel}}=\displaystyle \frac{{{\boldsymbol{v}}}_{{\rm{L}},\mathrm{gc}}-{{\boldsymbol{v}}}_{\odot ,\mathrm{gc}}}{{D}_{{\rm{L}}}}-\displaystyle \frac{{{\boldsymbol{v}}}_{{\rm{S}},\mathrm{gc}}-{{\boldsymbol{v}}}_{\odot ,\mathrm{gc}}}{{D}_{{\rm{S}}}}.\end{eqnarray} \tag{ 24 }$$

Here ${{\boldsymbol{v}}}_{{\rm{L}},\mathrm{gc}}$, ${{\boldsymbol{v}}}_{{\rm{S}},\mathrm{gc}}$, and ${{\boldsymbol{v}}}_{\odot ,\mathrm{gc}}$ are the galactocentric velocities of the lens, the source, and the Sun, respectively. Figure 2 demonstrates again that any knowledge of ${\tilde{v}}_{\mathrm{hel}}$ provides much more information of the lens distance than ${\mu }_{\mathrm{rel}}$ could.

Zoom In Zoom Out Reset image size

The distribution of ${D}_{8.3}$ is derived following a variant of the method in Calchi Novati et al. (2015). Here we provide the mathematical form of this derivation. For a fixed source distance ${D}_{{\rm{S}}}$, the differential event rate of Galactic microlensing is given by

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{4}{\rm{\Gamma }}}{{{dD}}_{{\rm{L}}}d\,\mathrm{log}\,{M}_{{\rm{L}}}{d}^{2}{{\boldsymbol{\mu }}}_{\mathrm{rel}}}={n}_{{\rm{L}},\star }{D}_{{\rm{L}}}^{2}(2{\theta }_{{\rm{E}}}){\mu }_{\mathrm{rel}}{f}_{\mu }({{\boldsymbol{\mu }}}_{\mathrm{rel}})\displaystyle \frac{d\xi ({M}_{{\rm{L}}})}{d\,\mathrm{log}\,{M}_{{\rm{L}}}}.\end{eqnarray} \tag{ 25 }$$

Here ${n}_{{\rm{L}},\star }$ is the local stellar density at position $(\alpha ,\delta ,{D}_{{\rm{L}}})$, ${f}_{\mu }({{\boldsymbol{\mu }}}_{\mathrm{rel}})$ is the two-dimensional probability distribution function of the lens–source relative proper motion ${{\boldsymbol{\mu }}}_{\mathrm{rel}}$, and $d\xi ({M}_{{\rm{L}}})/d\,\mathrm{log}\,{M}_{{\rm{L}}}$ is the stellar MF in logarithmic scale. Equation (25) can be rewritten in terms of microlensing observables $({D}_{8.3},{t}_{{\rm{E}}}^{{\prime} },{\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}})$,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{d}^{4}{\rm{\Gamma }}}{{{dD}}_{8.3}{{dt}}_{{\rm{E}}}^{{\prime} }{d}^{2}{\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}} & = & \displaystyle \frac{{d}^{4}{\rm{\Gamma }}}{{{dD}}_{{\rm{L}}}d\,\mathrm{log}\,{M}_{{\rm{L}}}{d}^{2}{{\boldsymbol{\mu }}}_{\mathrm{rel}}}\\ & & \times \left|\displaystyle \frac{\partial ({D}_{{\rm{L}}},\mathrm{log}{M}_{{\rm{L}}},{{\boldsymbol{\mu }}}_{\mathrm{rel}})}{\partial ({D}_{8.3},{t}_{{\rm{E}}}^{{\prime} },{\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}})}\right|\\ & = & 4{n}_{{\rm{L}},\star }\displaystyle \frac{{D}_{{\rm{L}}}^{4}}{{D}_{8.3}^{2}}{\mu }_{\mathrm{rel}}^{2}{f}_{\tilde{v}}({\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}})\displaystyle \frac{d\xi ({M}_{{\rm{L}}})}{d\,\mathrm{log}\,{M}_{{\rm{L}}}}\end{array}.\end{eqnarray} \tag{ 26 }$$

Here ${f}_{\tilde{v}}({\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}})$ is the two-dimensional probability function of ${\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}$, which can be derived from Equation (23) under a given Galactic model. In the latter evaluation of Equation (26), we have substituted Equation (25) and the following Jacobian determinant:

$$\begin{eqnarray}\left|\displaystyle \frac{\partial ({D}_{{\rm{L}}},\mathrm{log}{M}_{{\rm{L}}},{{\boldsymbol{\mu }}}_{\mathrm{rel}})}{\partial ({D}_{8.3},{t}_{{\rm{E}}}^{{\prime} },{\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}})}\right| & = & {\left(\displaystyle \frac{{D}_{{\rm{L}}}}{{D}_{8.3}}\right)}^{2}\displaystyle \frac{{\mu }_{\mathrm{rel}}}{{\tilde{v}}_{\mathrm{hel}}}\displaystyle \frac{1}{{M}_{{\rm{L}}}}\left|\displaystyle \frac{\partial ({M}_{{\rm{L}}},{\mu }_{\mathrm{rel}})}{\partial ({t}_{{\rm{E}}}^{{\prime} },{\tilde{v}}_{\mathrm{hel}})}\right|\\ & = & {\left(\displaystyle \frac{{D}_{{\rm{L}}}}{{D}_{8.3}}\right)}^{2}\displaystyle \frac{2{\pi }_{\mathrm{rel}}^{2}}{{\mathrm{au}}^{2}{t}_{{\rm{E}}}^{{\prime} }}.\end{eqnarray} \tag{ 27 }$$

For a given set of $({t}_{{\rm{E}}}^{{\prime} },{\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}})$, Equation (26) thus determines the relative (prior) probability distribution of ${D}_{8.3}$ at fixed ${D}_{{\rm{S}}}$. This is then integrated over the posterior distributions of ${t}_{{\rm{E}}}^{{\prime} }$ and ${\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}$ from the light-curve modeling to yield the relative probability distribution of ${D}_{8.3}$ for a fixed ${D}_{{\rm{S}}}$. To account for the variation in ${D}_{{\rm{S}}}$, we average over all possible values of ${D}_{{\rm{S}}}$ (from ${D}_{\min }=6\,\mathrm{kpc}$ to ${D}_{\max }=10\,\mathrm{kpc}$, assuming bulge sources), with each weighted by the number of available sources at that distance, ${n}_{{\rm{S}},\star }{D}_{{\rm{S}}}^{2-\gamma }{{dD}}_{{\rm{S}}}$. Here ${n}_{{\rm{S}},\star }$ is the local stellar density at $(\alpha ,\delta ,{D}_{{\rm{S}}})$, ${D}_{{\rm{S}}}^{2}{{dD}}_{{\rm{S}}}$ is the volume between ${D}_{{\rm{S}}}$ and ${D}_{{\rm{S}}}+{{dD}}_{{\rm{S}}}$, and ${D}_{{\rm{S}}}^{-\gamma }$ is approximately the fraction of stars that have the measured apparent magnitude (Kiraga & Paczynski 1994). We choose $\gamma =2.85$ for our sample, for reasons that are given in Appendix A. Then the non-normalized ("raw") probability distribution of ${D}_{8.3}$ for the given solution is

$$\begin{eqnarray}&&{P}_{\mathrm{raw}}({D}_{8.3})=\displaystyle \frac{{\int }_{{D}_{{\rm{S}},\min }}^{{D}_{{\rm{S}},\max }}{ \mathcal P }({D}_{8.3}| {D}_{{\rm{S}}}){n}_{{\rm{S}},\star }{D}_{{\rm{S}}}^{2-\gamma }{{dD}}_{{\rm{S}}}}{{\int }_{{D}_{{\rm{S}},\min }}^{{D}_{{\rm{S}},\max }}{n}_{{\rm{S}},\star }{D}_{{\rm{S}}}^{2-\gamma }{{dD}}_{{\rm{S}}}},\end{eqnarray} \tag{ 28 }$$

where

$$\begin{eqnarray}{ \mathcal P }({D}_{8.3}| {D}_{{\rm{S}}}) & \equiv & \displaystyle \int \displaystyle \frac{{d}^{4}{\rm{\Gamma }}}{{{dD}}_{8.3}{{dt}}_{{\rm{E}}}^{{\prime} }{d}^{2}{\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}}\\ & & \times P({t}_{{\rm{E}}}^{{\prime} }| \mathrm{Data})P({\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}| \mathrm{Data}){{dt}}_{{\rm{E}}}^{{\prime} }{d}^{2}{\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}.\end{eqnarray} \tag{ 29 }$$

In practice, we assume that the posterior distribution of ${t}_{{\rm{E}}}^{{\prime} }$, $P({t}_{{\rm{E}}}^{{\prime} }| \mathrm{Data})$, is a Dirac δ function and that the posterior distribution of ${\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}$, $P({\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}| \mathrm{Data})$, is a bivariate Gaussian function whose covariance matrix is determined in Section 4.1. The former assumption is reasonable because ${t}_{{\rm{E}}}^{{\prime} }$ (essentially ${t}_{{\rm{E}}}$) is well measured in almost all events and especially because it is much better constrained than ${\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}$. The second assumption is adopted so that the above integration can be computed analytically (see Appendix B). We have nevertheless tested the validity of this assumption with some examples, by comparing the analytic result with numerical integration of the (non-Gaussian) true posterior distribution from MCMC.

To derive the distance distribution of one event, we must weight all degenerate solutions correctly. The weight contains two factors: (1) $\exp (-{\rm{\Delta }}{\chi }^{2}/2)$, which is from the light-curve modeling, and (2) ${\pi }_{{\rm{E}}}^{-2}$, which is based on the so-called "Rich" argument. The Rich argument was originally pointed out by J. Rich (ca. 1997, private communication). It argues that, qualitatively, if ${\pi }_{{\rm{E}},+-}$ (and so ${\pi }_{{\rm{E}},-+}$) are much larger than ${\pi }_{{\rm{E}},++}$ (and ${\pi }_{{\rm{E}},--}$), then the former are likely spurious solutions. This is because the true solutions for events with small ${\pi }_{{\rm{E}}}$ are much more likely to be ${\pi }_{{\rm{E}},\pm \pm }$ solutions and can almost always generate spurious counterpart solutions ${\pi }_{{\rm{E}},\pm \mp }$ that are much larger. However, large ${\pi }_{{\rm{E}}}$ solutions can only rarely generate spurious small ${\pi }_{{\rm{E}}}$ solutions. Calchi Novati et al. (2015) quantified this argument and showed that solutions should be weighted by ${\pi }_{{\rm{E}}}^{-2}$, although they nevertheless only applied this weighting when the ratio of ${\pi }_{{\rm{E}}}$ between solutions was relatively large. Here, we carry the analysis of Calchi Novati et al. (2015) to its logical conclusion and apply this weighting uniformly to all events. The final normalized distribution of ${D}_{8.3}$ for each individual solution i is given by

$$\begin{eqnarray}&&{P}_{i}({D}_{8.3})=\displaystyle \frac{{e}^{-{\rm{\Delta }}{\chi }_{i}^{2}/2}{\pi }_{{\rm{E}},i}^{-2}{P}_{i,\mathrm{raw}}({D}_{8.3})}{{ \mathcal N }},\end{eqnarray} \tag{ 30 }$$

where

$$\begin{eqnarray}&&{ \mathcal N }\equiv \displaystyle \sum _{i=1}^{4}\displaystyle \frac{{e}^{-{\rm{\Delta }}{\chi }_{i}^{2}/2}}{{\pi }_{{\rm{E}},i}^{2}}\int {P}_{i,\mathrm{raw}}({D}_{8.3}){{dD}}_{8.3}.\end{eqnarray} \tag{ 31 }$$

The lens mass distribution is derived in a similar way. One key equation involved is the following:

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{4}{\rm{\Gamma }}}{d\,\mathrm{log}\,{M}_{{\rm{L}}}{{dt}}_{{\rm{E}}}^{{\prime} }{d}^{2}{\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}}=4{n}_{{\rm{L}},\star }{D}_{{\rm{L}}}^{4}{f}_{\tilde{v}}({\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}})\displaystyle \frac{d\xi ({M}_{{\rm{L}}})}{d\,\mathrm{log}\,{M}_{{\rm{L}}}}\displaystyle \frac{{\mu }_{\mathrm{rel}}^{3}}{{\tilde{v}}_{\mathrm{hel}}}.\end{eqnarray} \tag{ 32 }$$

The rest are almost identical to Equations (28) and (30), except for replacing ${D}_{8.3}$ with ${M}_{{\rm{L}}}$ (or $\mathrm{log}{M}_{{\rm{L}}}$).

We apply the planet sensitivity code developed by Zhu et al. (2015b). The method was first proposed by Rhie et al. (2000)²⁴ and further developed by Yee et al. (2015b) and Zhu et al. (2015b) to incorporate space-based observations. Below we provide brief descriptions of the methods and the code, and interested readers can find more details in Yee et al. (2015b) and Zhu et al. (2015b).

The calculation of planet sensitivity requires a certain value for ρ, which is the angular source size ${\theta }_{\star }$ normalized to ${\theta }_{{\rm{E}}}$, $\rho \equiv {\theta }_{\star }/{\theta }_{{\rm{E}}}$. The angular source size ${\theta }_{\star }$ is estimated following the standard procedure, i.e., by comparing the positions of the source star and the red clump centroid on the color–magnitude diagram (Yoo et al. 2004). The determination of ${\theta }_{{\rm{E}}}$ follows the prescription given by Yee et al. (2015b): for a given solution, we derive the transverse velocity ${\tilde{v}}_{\mathrm{hel}}$ using Equation (11) and choose ${\mu }_{\mathrm{rel}}=7\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$ if ${\tilde{v}}_{\mathrm{hel}}$ favors a disk lens and ${\mu }_{\mathrm{rel}}=4\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$ if ${\tilde{v}}_{\mathrm{hel}}$ favors a bulge lens; then ${\theta }_{{\rm{E}}}={\mu }_{\mathrm{rel}}{t}_{{\rm{E}}}$.

We first compute the planet sensitivity S as a function of planet-to-star mass ratio q and the planet/star separation s normalized to the angular Einstein radius ${\theta }_{{\rm{E}}}$. Twenty q values are chosen uniformly in logarithmic scale between 10⁻⁵ and 0.04, which correspond to a mass range from $1\,{M}_{\oplus }$ to $13\,{M}_{{\rm{J}}}$ for a $0.3\,{M}_{\odot }$ host. Twenty s values are chosen also uniformly in logarithmic scale between 0.3 and 3. Our choice of the "lensing zone" covers the region where microlensing is sensitive for nearly all the events. For each set of (q, s), we generate 100 planetary light curves that have other parameters the same except for α, which is the angle between the source trajectory and the lens binary axis. For each simulated light curve, we then find the best-fit single-lens model using the downhill simplex algorithm, the goodness of which is quantified by ${\chi }_{\mathrm{SL}}^{2}$. For events that were subjectively chosen and never met the objective criteria, we additionally find the deviation from the single-lens model in the ground-based data that were released²⁵ before the subjective chosen date ${t}_{\mathrm{sub}}$. If this deviation is significant (${\chi }^{2}\gt 10$; Yee et al. 2015b), we consider the injected planet as having been noticeable and thus reject this α, regardless of how significant ${\chi }_{\mathrm{SL}}^{2}$ is. Otherwise, for these events and events that met objective criteria, we pass the simulated events to the anomaly detection filter. The sensitivity $S(q,s)$ is the fraction of α values for which the injected planets are detectable.

We adopt the following detection thresholds, which are more realistic than that used in Zhu et al. (2015b) and have been used in Poleski et al. (2016): ${\chi }_{\mathrm{SL}}^{2}\gt 300$ and at least three consecutive data points from the same observatory show $\gt 3\sigma$ deviations (C1), or ${\chi }_{\mathrm{SL}}^{2}\gt 500$ (C2). C1 aims for capturing sharp planetary anomalies, and C2 is supplementary to C1 for recognizing the long-term weak distortions.

In principle, the planet sensitivities could be substantially different for the ${\pi }_{{\rm{E}},\pm \pm }$ solutions compared to the ${\pi }_{{\rm{E}},\pm \mp }$ solutions, because source trajectories as seen by Earth and Spitzer pass by the lens on the projected plane from the same side for the former, but on opposite sides for the latter (Zhu et al. 2015b). However, for the data sets under consideration in the present paper, which typically have several dozen observations per day from the ground and only one or a few per day from space, almost all the sensitivity comes from the ground observations. Hence, the sensitivities of the four degenerate solutions are almost identical. See Figure 6 of Poleski et al. (2016) for an example. The small differences between four solutions arise from the different values of ρ used in the computation, because $\rho ={\theta }_{\star }/({\mu }_{\mathrm{rel}}{t}_{{\rm{E}}})$ and the choice of ${\mu }_{\mathrm{rel}}$ rely on the magnitude of ${\pi }_{{\rm{E}}}$.

Current experiments are very far from having the ability to separately measure distance distributions for the individual (s,q). Hence, we also define the sensitivity to a given planet-to-star mass ratio q as

$$\begin{eqnarray}&&S(q)=\int S(q,s)d\,\mathrm{log}\,s.\end{eqnarray} \tag{ 33 }$$

This bears the assumption that the distribution of s is flat in logarithmic scale, which is reasonable according to recent studies (e.g., Dong & Zhu 2013; Fressin et al. 2013; Petigura et al. 2013; Burke et al. 2015; Clanton & Gaudi 2016).

We present the ground-based and space-based light curves of each event in our sample in Figures 3–5. All data sets except OGLE have been rescaled to the OGLE I magnitude based on the best-fit model

$$\begin{eqnarray}&&{\tilde{F}}^{j}=\displaystyle \frac{{F}_{{\rm{s}}}^{\mathrm{OGLE}}}{{F}_{{\rm{s}}}^{j}}({F}^{j}-{F}_{{\rm{b}}}^{j})+{F}_{{\rm{b}}}^{\mathrm{OGLE}}.\end{eqnarray} \tag{ 34 }$$

We suppress KMTNet data sets in these plots and only show OGLE data for clarity. The reader can find an example event that demonstrates the much denser coverage of KMTNet in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The ground-based data of all 50 events in our sample can be well fitted by a single-lens model. However, the Spitzer data of several of them show deviations from this simple description. Some of these deviations are prominent, such as in OGLE-2015-BLG-0081, OGLE-2015-BLG-0461, OGLE-2015-BLG-0703, OGLE-2015-BLG-0961, and OGLE-2015-BLG-1189. However, we believe that these are due to unknown systematics in the Spitzer data rather than indications of companions to the lens. Below we provide two examples to demonstrate this point. Poleski et al. (2016) noticed a strong deviation in the Spitzer data of OGLE-2015-BLG-0448. Although they found that a lens companion with $q=1.7\times {10}^{-4}$ could improve the single-lens model by ${\rm{\Delta }}{\chi }^{2}=128$, they showed that even the best-fit binary-lens model could not remove all the deviations in the Spitzer data. Therefore, the trend in Spitzer data was likely caused by systematics rather than a physical signal from an additional lens object. This is especially true for OGLE-2015-BLG-0961. As shown in Figure 6, the ground-based data can be well fitted by a single-lens model with extremely high magnification (${u}_{0,\oplus }\leqslant 0.005$ at 1σ level), which excludes any lens companions with $q\gtrsim 3\times {10}^{-4}$ if close to the Einstein ring (see Figure 13 below). The Spitzer data show a long-term deviation centered at the time when the event peaked from the ground. This long-term deviation, if attributed to a companion to the lens, could only be caused by the planetary caustic. With ${u}_{0,\mathrm{spitz}}=0.1$ and the position of planetary caustic at $| s-1/s|$, the separation between the hypothetical lens companion and the primary lens should be $\mathrm{log}s=\pm 0.02$. Combining the duration of the deviation (10 days out of ${t}_{{\rm{E}}}=60$ days) and the width of the planetary caustic (Han 2006), we can put a limit on the companion mass ratio $q\gtrsim 2\times {10}^{-3}$. There do not exist any q values that could explain the nondetection in the ground-based data and the significant trend in Spitzer data. Therefore, the trend in Spitzer data is likely due to systematics in the Spitzer data reduction.²⁶

The systematics in Spitzer data can potentially affect the parallax measurements. However, it has been demonstrated that the influence is small in several published events. For example, Poleski et al. (2016) showed that the parallax parameters with and without the systematic trend were almost identical. The agreement between orbital parallax and satellite parallax also indicates that the effect of systematics is less likely an issue (e.g., Udalski et al. 2015b; Han et al. 2017).

We provide in Table 2 the best-fit parameters and associated uncertainties for solutions with ${\rm{\Delta }}{\chi }^{2}\leqslant 100$ of all 50 events. Here ${\rm{\Delta }}{\chi }^{2}$ is the difference between a given solution and the best solution for that event. With these, and following the method in Section 4.3, we derive the lens distance parameter ${D}_{8.3}$ and lens mass ${M}_{{\rm{L}}}$ distributions for every event in our raw sample.

Table 2.  Best-fit Parameters and Associated Uncertainties for the 50 Events in the Raw Sample

OGLE #	Solution	${\rm{\Delta }}{\chi }^{2}$	t0	u0	${t}_{{\rm{E}}}$	${\pi }_{{\rm{E}},{\rm{N}}}$	${\pi }_{{\rm{E}},{\rm{E}}}$	${I}_{\mathrm{base}}$ a	Blending	$I-[3.6\mu {\rm{m}}]$
0011	(+, +)	0.0	7125.99(12)	0.8919(5)	80.96(26)	0.262(6)	0.104(4)	15.948	⋯	2.494(26)
	(−, −)	31.3	7125.53(10)	−0.8919(5)	74.0(3)	−0.403(9)	0.089(4)	15.948	⋯	2.529(24)
0029	(+, +)	0.0	7233.03(7)	0.7603(3)	101.77(11)	−0.1199(17)	−0.1192(17)	15.132	⋯	2.005(12)
0034	(+, +)	0.0	7140.37(3)	0.282(4)	64.2(4)	0.055(29)	−0.057(5)	16.559	0.009(18)	3.06(4)
	(−, −)	1.4	7140.37(3)	−0.283(3)	63.3(5)	−0.13(5)	−0.061(4)	16.559	0.001(15)	3.06(4)
0081	(+, +)	30.3	7135.4(4)	0.7014(14)	85.5(14)	−0.345(25)	−0.090(18)	17.892	⋯	2.758(14)
	(−, −)	0.0	7135.7(5)	−0.7037(13)	97.7(24)	0.475(28)	−0.130(21)	17.892	⋯	2.74(3)
0350	(+, +)	5.1	7214.86(12)	0.445(19)	70.7(20)	−0.049(3)	−0.068(7)	17.295	0.12(5)	2.626(21)
	(−, −)	0.0	7214.99(12)	−0.479(19)	66.7(18)	0.050(4)	−0.090(8)	17.295	0.01(5)	2.617(23)
	(−, +)	19.4	7214.50(11)	−0.448(17)	65.0(17)	0.567(20)	−0.205(9)	17.295	0.11(4)	2.595(23)
0379	(+, +)	2.8	7189.62(6)	0.448(15)	64.4(14)	0.003(11)	−0.098(7)	18.087	0.18(3)	2.61(3)
	(−, −)	0.0	7189.63(6)	−0.461(15)	62.7(14)	−0.001(13)	−0.119(7)	18.087	0.13(4)	2.59(5)
0388	(+, +)	0.0	7161.55(7)	0.508(21)	40.3(11)	0.11(3)	0.048(29)	17.425	0.25(4)	3.67(6)
	(−, −)	2.2	7161.55(7)	−0.518(23)	39.8(12)	−0.15(5)	0.018(29)	17.425	0.25(4)	3.65(6)
	(−, +)	2.1	7161.66(9)	−0.478(18)	46.5(11)	1.40(7)	−0.27(5)	17.425	0.32(3)	3.63(5)
0461	(+, +)	1.8	7161.27(19)	0.87(5)	41.6(17)	−0.20(5)	0.05(4)	17.414	−0.01(10)	1.76(7)
	(−, −)	0.0	7161.31(20)	−0.86(6)	42.9(23)	0.31(8)	0.05(4)	17.414	0.04(10)	1.71(9)

Note.  For each event, we only include solutions that have ${\rm{\Delta }}{\chi }^{2}\lt 100$ from the lowest value. We present the total baseline magnitude in OGLE-IV I band, ${I}_{\mathrm{base}}$, the blending fraction in I band, and the source $I-[3.6\ \mu {\rm{m}}]$ color, rather than the flux parameters $({F}_{{\rm{s}}},{F}_{{\rm{b}}})$ of individual data sets. We assume no blending for events OGLE-2015-BLG-0011, OGLE-2015-BLG-0029, OGLE-2015-BLG-0081, OGLE-2015-BLG-0772, OGLE-2015-BLG-0798, OGLE-2015-BLG-1096, OGLE-2015-BLG-1188, OGLE-2015-BLG-1289, OGLE-2015-BLG-1297, OGLE-2015-BLG-1341, and OGLE-2015-BLG-1470, because free blending would lead to severely negative blending.

^aThe uncertainty of ${I}_{\mathrm{base}}$ is ∼1 mmag, primarily arising from OGLE-IV's data-recording format. In addition, the calibration precision of OGLE-IV I band to the standard system, ∼10 mmag, is not included here, on the basis that it does not affect the determination of microlensing parameters.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Based on all event parameters and the subsequent lens distributions, we can now select events for our final statistical sample. The guideline is that only events with "detected parallax" can be included for the study of the Galactic distribution of planets, as Yee et al. (2015b) pointed out. At first sight, the above guideline seems to suggest a criterion on the measurement uncertainty of ${{\boldsymbol{\pi }}}_{{\rm{E}}}$. However, such an approach would be problematic, in particular because the uncertainty of ${{\boldsymbol{\pi }}}_{{\rm{E}}}$ is determined for individual solutions, but decisions have to be made for individual events, which generally have more than one solution. As shown in Table 2, the $(\pm ,\,\mp )$ solutions are in general better constrained than the $(\pm ,\,\pm )$ solutions, so even though they are statistically disfavored by the Rich argument, they are more likely to survive if a cut on the detection significance of ${\pi }_{{\rm{E}}}$ is applied. Although it is possible to design a criterion for choosing events that balances the two opposite factors, a better approach is to choose events based on the distance parameter ${D}_{8.3}$ and its associated uncertainty $\sigma ({D}_{8.3})$. This is because only events with well-determined distances contribute to the measurement of the Galactic distribution of planets.

We show in Figure 7 the median value and the 1σ uncertainty of the lens distance parameter ${D}_{8.3}$ derived for each event in our raw sample. Here the 1σ uncertainty is the half-width of a 68% confidence interval centered on the median ${D}_{8.3}$. By visually inspecting the ${D}_{8.3}$ distributions of all 50 events, which are shown in Figure 8, we decide to use $\sigma ({D}_{8.3})\leqslant 1.4\,\mathrm{kpc}$ as the criterion for claiming a parallax detection and thus for any event to be included in the final sample. We end up with 41 events in the final sample. The nine events that are excluded all have broad distributions of ${D}_{8.3}$, even though some of them have very good measurements of ${\pi }_{{\rm{E}}}$ (e.g., OGLE-2015-BLG-0029, OGLE-2015-BLG-0843, OGLE-2015-BLG-1167). The broad distribution of ${D}_{8.3}$ arises from the atypical magnitude and direction of ${\tilde{{\boldsymbol{v}}}}_{\mathrm{hel}}$. When combined with the Galactic model, the former favors near- to mid-disk lenses, while the latter favors bulge lenses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The derived lens masses and the fractional uncertainties are shown in Figure 9. As expected, events that do not show compact ${D}_{8.3}$ distributions do not have well-constrained mass estimates either. For events in our final sample, the typical uncertainty of the lens mass estimate is 20%, regardless of whether it is substellar or not. In particular, we note that the lens mass estimate of OGLE-2015-BLG-1482 agrees reasonably well ($\leqslant 2\sigma$) with the direct mass measurement from the finite-source effect (S.-J. Chung et al. 2017) as a demonstration that the mass estimate method employed here is valid. The derived lens mass distributions of all 50 events are presented in Figure 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We show in Figure 11 the cumulative distributions of event timescales ${t}_{{\rm{E}}}$ and impact parameters ${u}_{0,\oplus }$ as seen from the ground, and we compare them with those in the OGLE-III microlensing event catalog (Wyrzykowski et al. 2015), which can be considered to be complete and uniform for our purpose. Because different solutions have only slightly different ${t}_{{\rm{E}}}$ and ${u}_{0,\oplus }$, we simply take values of the solution with lowest ${\chi }^{2}$. For the timescale ${t}_{{\rm{E}}}$ distribution, we notice that events in our final sample are more concentrated within 10–100 days than events in the OGLE-III catalog. The lower limit comes into play because there is a 3- to 9-day lag between events being selected and events getting observed by Spitzer (see Figure 1 of Udalski et al. 2015b). The lack of extremely long timescale (${t}_{{\rm{E}}}\gtrsim 100$ days) events comes as a consequence of our event selection criteria, because a substantial brightness change ($\gtrsim 0.3$ mag) during the ∼40-day Spitzer bulge window is required in order to detect the parallax effect (Yee et al. 2015b). Although the events in our sample (and subsequent larger samples) have a biased ${t}_{{\rm{E}}}$ distribution, this bias applies to both events with and without planet detections in the same way. Therefore, it will not affect the statistical studies of the Galactic distribution of planets. The Spitzer sample u₀ distribution shows similar overall morphology to that of the OGLE-III catalog but is more uniform, which indicates that it shows less magnification bias. This again reflects the fact that OGLE-III detections are possible based on a few days of relatively magnified sources, whereas Spitzer selections are delayed by 3–9 days.

Zoom In Zoom Out Reset image size

We show in Figure 12 the distributions of lens distance parameter ${D}_{8.3}$ and lens mass ${M}_{{\rm{L}}}$, which are averaged over all 41 events in the final sample. We consider the influences of the Rich argument and a different choice of stellar MF (e.g., Kroupa MF). As expected (see Section 4.1; see also Calchi Novati et al. 2015), the lens distance distribution is biased toward more nearby and therefore lower-mass lenses, if the Rich argument is not taken into account. The different choices of the stellar MF have a marginal effect, especially on the lens distance distribution. Our result demonstrates, for the first time, that the peak of the microlens mass distribution is at $0.5\,{M}_{\odot }$ and that the majority microlensing events are caused by M dwarfs.

Zoom In Zoom Out Reset image size

We present in Figures 13 and 14 the planet sensitivity plots of individual events in our final sample. Events are divided according to their final status of Spitzer selections, with objectively selected events shown in Figure 13 and subjectively selected events shown in Figure 14. For all objective events and most subjective events, the sensitivity curves are smooth and triangle-like, with either a single horn (for relatively high magnification events; see also Gould et al. 2010) or double horns (for relatively low magnification events; see also Gaudi et al. 2002). In the remaining subjective events, however, the sensitivity curves show discontinuity especially at large q values. This was caused by the way that the planet sensitivity of a subjectively chosen event was computed. As described in detail in Yee et al. (2015b) and Zhu et al. (2015b) and summarized in Section 4.4, for events that were chosen subjectively and never met the objective selection criteria, all (hypothetical) planet detections must be censored from the statistical sample if they would have betrayed their existence in the data that were released before the subjective selection date ${t}_{\mathrm{sub}}$. This has only a marginal effect if ${t}_{\mathrm{sub}}$ is well before the event peak ${t}_{0,\oplus }$, because the bulk of planet sensitivities come from the region near the peak ($| t-{t}_{0,\oplus }| \lesssim {u}_{0,\oplus }{t}_{{\rm{E}}}$). If ${t}_{\mathrm{sub}}$ is close to ${t}_{0,\oplus }$, then the above procedure could affect the final sensitivity curves significantly. In particular, planets that are more massive and closer to the Einstein ring are more easily excluded in the sensitivity computation; for given combinations of q and s, some choices of α are more easily discarded as well. As an example, we show in Figure 15 the ${\chi }^{2}$ maps for three different q values for two events, OGLE-2015-BLG-0987 and OGLE-2015-BLG-1189, which have similar impact parameters ${u}_{0,\oplus }$ but show very different sensitivity curves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We provide constraints on the planet distribution function, based on the null detection in our sample. We adopt as the planet distribution function the form

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\,\mathrm{log}q}={ \mathcal A }{\left(\displaystyle \frac{q}{{q}_{\mathrm{ref}}}\right)}^{\alpha },\end{eqnarray} \tag{ 35 }$$

and we choose ${q}_{\mathrm{ref}}=5\times {10}^{-4}$, which is the typical q value of microlensing planets (e.g., Gould et al. 2010). We first show in the left panel of Figure 16 the sensitivity curves averaged over the 41 events in the final sample. Assuming Poisson-like noise and that "planets" should have $q\leqslant {10}^{-2}$ (to be consistent with previous studies; e.g., Gould et al. 2010), we are able to derive the constraints on the slope of the planet MF α and the normalization factor ${ \mathcal A }$ based on the null detection in our sample. The results are shown in the right panel of Figure 16. Our constraints are consistent, at the 2σ level, with previous statistical studies based on samples of microlensing planets (Gould et al. 2010; Shvartzvald et al. 2016; Suzuki et al. 2016). In particular, we find ${ \mathcal A }\lt 0.49$ at the 95% confidence level for a flat ($\alpha =0$) planet MF, which is consistent with the result (${ \mathcal A }=0.36\pm 0.15$) from Gould et al. (2010).

Zoom In Zoom Out Reset image size

We derive the cumulative distribution of planet sensitivities of our sample based on the lens distribution $P({D}_{8.3})$ and the planet sensitivity S(q). The results are presented here in terms of both the planet-to-star mass ratio q,

$$\begin{eqnarray}{{ \mathcal C }}_{q}({D}_{8.3}) & = & \displaystyle \frac{1}{{{ \mathcal C }}_{q}({R}_{\mathrm{GC}})}\displaystyle \sum _{i,j}{\int }_{0}^{{D}_{8.3}}{P}_{i}^{j}(D^{\prime} ){dD}^{\prime} \\ & & \times {\int }_{{q}_{\min }}^{{q}_{\max }}{S}_{i}^{j}(q)d\,\mathrm{log}\,q,\end{eqnarray} \tag{ 36 }$$

and the planet mass ${m}_{{\rm{p}}}$,

$$\begin{eqnarray}{{ \mathcal C }}_{m}({D}_{8.3}) & = & \displaystyle \frac{1}{{{ \mathcal C }}_{m}({R}_{\mathrm{GC}})}\\ & & \times \displaystyle \sum _{i,j}{\displaystyle \int }_{0}^{{D}_{8.3}}{P}_{i}^{j}(D^{\prime} ){dD}^{\prime} {\displaystyle \int }_{{q}_{\min }(D^{\prime} )}^{{q}_{\max }(D^{\prime} )}{S}_{i}^{j}(q)d\,\mathrm{log}\,q.\end{eqnarray} \tag{ 37 }$$

Here ${P}_{i}^{j}({D}_{8.3})$ and ${S}_{i}^{j}(q)$ are the lens distance distribution and the planet sensitivity for solution i of event j, respectively. In Equation (36), we choose ${q}_{\min }={10}^{-5}$ and ${q}_{\max }={10}^{-2}$. In Equation (37), we solve for the boundaries on q for individual ${D}_{8.3}$ values that lead to the planet mass ranging from $1\,{M}_{\oplus }$ to $3\,{M}_{{\rm{J}}}$. These two distributions are normalized so that ${ \mathcal C }({R}_{\mathrm{GC}})=1$. The results are shown as the solid black curves in Figure 17.

Zoom In Zoom Out Reset image size

In terms of the Galactic distribution of planets, we derive the ratio of planets in the bulge to planets in the disk (to which our survey is sensitive),

$$\begin{eqnarray}&&{\eta }_{{\rm{b}}2{\rm{d}}}\equiv \displaystyle \frac{{\int }_{0}^{{R}_{\mathrm{GC}}}{ \mathcal C }^{\prime} ({D}_{8.3}){f}_{B}({D}_{8.3})d\,{D}_{8.3}}{{\int }_{0}^{{R}_{\mathrm{GC}}}{ \mathcal C }^{\prime} ({D}_{8.3})[1-{f}_{B}({D}_{8.3})]d\,{D}_{8.3}}.\end{eqnarray} \tag{ 38 }$$

Here ${ \mathcal C }^{\prime} ({D}_{8.3})$ is the derivative of ${ \mathcal C }({D}_{8.3})$, and ${f}_{B}({D}_{8.3})$ is the contribution of bulge events (lens in the bulge) to all events at given ${D}_{8.3}$, which is given by

$$\begin{eqnarray}&&{f}_{B}({D}_{8.3})=\displaystyle \frac{{\int }_{{D}_{{\rm{S}},\min }}^{{D}_{{\rm{S}},\max }}{n}_{{\rm{B}}}({D}_{{\rm{L}}}){/n}_{\star }({D}_{{\rm{L}}}){n}_{\star }({D}_{{\rm{S}}}){D}_{{\rm{S}}}^{2-\gamma }{D}_{{\rm{L}}}^{2}d\,{D}_{{\rm{S}}}}{{\int }_{{D}_{{\rm{S}},{\min }}}^{{D}_{{\rm{S}},\max }}{n}_{\star }({D}_{{\rm{S}}}){D}_{{\rm{S}}}^{2-\gamma }{D}_{{\rm{L}}}^{2}d\,{D}_{{\rm{S}}}}.\end{eqnarray} \tag{ 39 }$$

Here ${D}_{{\rm{L}}}$ is derived for given ${D}_{8.3}$ and ${D}_{{\rm{S}}}$. As illustrated in Figure 17, for the current sample we find that ${\eta }_{{\rm{b}}2{\rm{d}}}=28 \%$ if "planet" is defined by mass ratio q in the range 10⁻⁵ to 10⁻², and that ${\eta }_{{\rm{b}}2{\rm{d}}}=35 \%$ if "planet" is defined by mass in the range $1\,{M}_{\oplus }$ to $3\,{M}_{{\rm{J}}}$. Assuming that the planet formation is no different between the bulge and the disk, these suggest that ∼1/3 of all planet detections in our experiment should come from bulge events. In other words, any deviation from the above value would indicate that the bulge planet population is different from the disk planet population.

We also investigate the influence of impact parameters on the two cumulative distributions ${{ \mathcal C }}_{q}({D}_{8.3})$ and ${{ \mathcal C }}_{m}({D}_{8.3})$, for the purpose of better planning future similar experiments. We find that the sample of events with maximum magnifications ${A}_{\max }\gt 8$, which account for 24% of all events in our sample, contributes 40%–45% of all planet sensitivities.