This site uses cookies. By continuing to use this site you agree to our use of cookies. To find out more, see our Privacy and Cookies policy.

THE EXOPLANET MASS-RATIO FUNCTION FROM THE MOA-II SURVEY: DISCOVERY OF A BREAK AND LIKELY PEAK AT A NEPTUNE MASS

, , , , , , , , , , , , , , , , , , , , , , , , , , , , , and

Published 2016 December 13 © 2016. The American Astronomical Society. All rights reserved.
, , Citation D. Suzuki et al 2016 ApJ 833 145 DOI 10.3847/1538-4357/833/2/145

0004-637X/833/2/145

ABSTRACT

We report the results of the statistical analysis of planetary signals discovered in MOA-II microlensing survey alert system events from 2007 to 2012. We determine the survey sensitivity as a function of planet–star mass ratio, q, and projected planet–star separation, s, in Einstein radius units. We find that the mass-ratio function is not a single power law, but has a change in slope at q ∼ 10−4, corresponding to ∼20 M for the median host-star mass of ∼0.6 ${M}_{\odot }$. We find significant planetary signals in 23 of the 1474 alert events that are well-characterized by the MOA-II survey data alone. Data from other groups are used only to characterize planetary signals that have been identified in the MOA data alone. The distribution of mass ratios and separations of the planets found in our sample are well fit by a broken power-law model of the form ${{dN}}_{\mathrm{pl}}/{(d\mathrm{log}qd\mathrm{log}s)=A(q/{q}_{\mathrm{br}})}^{n}{s}^{m}\,{\mathrm{dex}}^{-2}$ for q > qbr and ${{dN}}_{\mathrm{pl}}/{(d\mathrm{log}qd\mathrm{log}s)=A(q/{q}_{\mathrm{br}})}^{p}{s}^{m}\,{\mathrm{dex}}^{-2}$ for q < qbr, where qbr is the mass ratio of the break. We also combine this analysis with the previous analyses of Gould et al. and Cassan et al., bringing the total sample to 30 planets. This combined analysis yields $A={0.61}_{-0.16}^{+0.21}$, n = −0.93 ± 0.13, $m={0.49}_{-0.49}^{+0.47}$, and $p={0.6}_{-0.4}^{+0.5}$ for qbr ≡ 1.7 × 10−4. The unbroken power-law model is disfavored with a p-value of 0.0022, which corresponds to a Bayes factor of 27 favoring the broken power-law model. These results imply that cold Neptunes are likely to be the most common type of planets beyond the snow line.

Export citation and abstract BibTeX RIS

1. INTRODUCTION

More than 50 exoplanets have been discovered by the gravitational microlensing method (Mao & Paczyński 1991) since the first discovery more than a decade ago (Bond et al. 2004). This is a small fraction of the number of planets found by the transit (Mullally et al. 2015) and radial velocity (RV; Butler et al. 2006) methods, but the microlensing method can detect a class of planets that are invisible to other methods. Microlensing is uniquely sensitive to planets down to an Earth mass (Bennett & Rhie 1996) beyond the snow line, where the other detection methods are ineffective (Bennett 2008; Gaudi 2012). According to the core accretion theory (Lissauer 1993; Ida & Lin 2004a), the snow line is the radius in the protoplanetary disk beyond which the temperature is cold enough for ices to condense. The presence of ices increases the density of solids in the protoplanetary disk by a factor of 3–4 (Kennedy & Kenyon 2008), and this is thought to allow the ice-rock cores of giant planets to form more quickly than solid planets can condense in the inner parts of protoplanetary disks. The core accretion theory predicts (Ida & Lin 2004a) that gas-giant planet cores can grow to ∼10 Earth masses of solid material beyond the snow line, and this is enough to trigger the rapid accretion of hydrogen and helium gas needed to form a gas giant. However, there is a danger that the hydrogen and helium might dissipate before the core grows large enough to trigger rapid accretion. In this case, we are likely to be left with a failed gas-giant-core planet. It is thought that these failed gas-giant-core planets may be especially common around low-mass host stars (Laughlin et al. 2004). The early results from microlensing support this idea. Microlensing has found numerous examples of planets that are consistent with the expectations for these failed gas-giant-core planets (Beaulieu et al. 2006; Gould et al. 2006, 2014; Muraki et al. 2011; Furusawa et al. 2013; Batista et al. 2015; Bennett et al. 2015). Other methods have not been able to detect such planets, so the microlensing method is uniquely sensitive to this population of failed gas-giant-core planets.

The microlensing method is also uniquely sensitive to old, free-floating planets that may have been ejected from the planetary system of their formation. While it is possible to detect young, unbound planetary mass objects directly in the infrared (Delorme et al. 2012; Beichman et al. 2013), the infrared surveys do not yet have the sensitivity to detect planets of ∼1 Jupiter mass, where the signal of a large free-floating, or rogue, planet population has been found by microlensing (Sumi et al. 2011). Some have suggested that this population of apparently rogue planets might be bound in orbits so wide that microlensing cannot detect the host stars (Quanz et al. 2012), but the microlensing results indicate that the median projected distance from rogue planets to potential host stars is likely to be >30 au (Bennett et al. 2012).

There have been several previous statistical studies of exoplanets found by microlensing. Statistical studies (Gaudi et al. 2002; Snodgrass et al. 2004) from early microlensing survey data with no planet detections were able to put upper limits on only the cold gas-giant planet frequency. Using 10 microlensing planet detections, Sumi et al. (2010) looked at the distribution of planets as a function of their mass ratio, q, and found that the exoplanet distribution beyond the snow line could be explained by a power-law mass-ratio function of the form ${dN}/d\,\mathrm{log}q={q}^{-0.68\pm 0.20}$, but this calculation was done without a calculation of the absolute exoplanet detection efficiency. (Throughout this paper, we use log to refer to the base-10 logarithm.) They used the detection efficiency for each planetary event, but did not calculate the efficiency for events without planetary signals. This allowed the determination of the relative efficiency for different types of planets, but it was insufficient to determine the absolute detection efficiency. So, Sumi et al. (2010) could only calculate the power-law index of the mass-ratio function, but not its normalization. Gould et al. (2010) took advantage of the very high detection efficiency provided by high-magnification microlensing events (Griest & Safizadeh 1998; Rhie et al. 2000) to do the first calculation of the absolute detection efficiency with a sample of 13 high-magnification microlensing events including six planets. While the very high planet detection efficiency is provided by high-magnification events, there is a well-known degeneracy between close and wide planet solutions (Dominik 1999; Chung et al. 2005) that implies that, for most high-magnification events, there is a significant ambiguity in the star–planet separation, unless it is close to the Einstein radius. As a result, it is difficult to measure the separation dependence of the planet population with a sample of high-magnification events, and Gould et al. (2010) did not attempt to do so. They found a planet frequency of ${d}^{2}N/d\,\mathrm{log}q\,d\,\mathrm{log}s=0.36\pm 0.15$ at q = 5 × 10−4 beyond the snow line.

The PLANET collaboration (Cassan et al. 2012) also published a statistical analysis of planetary events found in their 2002–2007 data. They had a much larger sample of microlensing events, but the number of planets discovered was smaller, three, compared to six planets found in the Gould et al. (2010) sample. Because one planet is common to both samples, the combined sample includes eight planets, and Cassan et al. (2012) use this combined sample and the slope of the power-law mass-ratio function from Sumi et al. (2010) to find the cold-planet mass function ${d}^{2}N/d\,\mathrm{log}M\,d\,\mathrm{log}a={0.24}_{-0.10}^{+0.16}{(M/{M}_{\mathrm{Sat}})}^{-0.73\pm 0.17}$, where MSat = 95.2 M. This result is quoted in terms of planetary mass instead of mass ratio because Cassan et al. (2012) used a Bayesian method (Dominik 2006) to convert from the mass ratios, q, measured in the microlensing light curves to masses. This method implicitly assumes that the probability for a star to host a planet of a given mass ratio is independent of the stellar mass and location of the system through the Galaxy. This assumption is probably safe with such a small sample, but it could be problematic with samples large enough to measure the dependence of the exoplanet mass function on the mass of and distance to the host star.

A more recent paper is based on three-site survey data using combined data from the Optical Gravitational Lensing Experiment (OGLE), MOA, and Wise wide-field microlensing surveys (Shvartzvald et al. 2016). This has some overlap with our own sample, as it includes data from 2011 to 2014, but it only includes events from the 8 deg2 covered by all three surveys. They used a coarse grid in mass ratio and separation for the light-curve modeling of some of the binary and planetary events found in their sample. This appears to have resulted in important errors in the parameters for some events, as we discuss in Section 6.1. Nevertheless, they estimate a slope of ${dN}/d\,\mathrm{log}q={q}^{-0.50\pm 0.17}$. They also estimated the total exoplanet frequency for 10−4.9 < q < 10−1.4.

Recent studies (Clanton & Gaudi 2014a, 2014b; Montet et al. 2014) have compared the exoplanet distribution found by microlensing with the results of RV observations of M-dwarfs and found that the results from both methods are consistent, although the RV is only sensitive to planets of Jupiter mass or greater beyond the snow line.

The Kepler telescope (Borucki et al. 2010; Mullally et al. 2015) has revealed the planet population in close orbits, mostly within 1 au around solar-type stars. Most of detected planets or candidates are located within the orbit of Mercury and many multiple planets systems consisting of super Earths have been found (Howard et al. 2012). Microlensing also has detected a number of super Earths (Beaulieu et al. 2006; Gould et al. 2006, 2014; Bennett et al. 2008; Muraki et al. 2011), but the planets found by microlensing are in much longer period orbits around somewhat lower mass stars than the small planets found by Kepler. Also, the masses for small Kepler planets with longer orbital periods are rarely determined, and when they are determined, they show surprisingly low densities (Lissauer et al. 2013). Attempts to measure a mass–radius relation find a scatter much larger than the measurement errors (Weiss & Marcy 2014; Wolfgang et al. 2016), indicating a range of compositions for exoplanets at a fixed mass or radius. The average exoplanet composition could change dramatically between the well-measured Kepler planets at ≪1 au and microlens planets beyond the snow line. Thus, it is difficult to compare statistical analyses of exoplanet frequency from Kepler data with those from current, ground-based microlensing survey data. Future space-based exoplanet microlensing surveys (Bennett & Rhie 2002; Penny et al. 2013), such as the WFIRST survey (Spergel et al. 2015), will explore a much wider range of orbital periods and will therefore have overlap with Kepler.

The one observable that relates to interesting lens properties in single-mass microlensing light curves is the event timescale, the time required for a source star to cross the angular Einstein radius,

Equation (1)

where μrel is the relative lens-source proper motion, M is mass of the lens system, and DL and DS are the distance to the lens and source star, respectively. In planetary microlensing events, where signals of both the planet and its host star are seen, we can also measure the planet–star mass ratio and angular separation relative to the angular Einstein radius. In most planetary events, the sharp planetary light-curve features resolve the finite angular size of the source star, allowing the source radius crossing time, t*, to be measured. Since the angular radius of the source star, θ*, can be determined from its brightness and color (Boyajian et al. 2014), the measurement of θ* also provides a measurement of the angular Einstein radius, θE = θ*tE/t*. Unfortunately, this is not enough to determine the mass or distance of the lens system. Instead, we are left with the following mass–distance relation:

Equation (2)

where x = DL/DS. In some cases, the actual mass can be determined if a microlensing parallax signal (Gould 1992; Alcock et al. 1995) is measured (Bennett et al. 2008; Gaudi et al. 2008; Muraki et al. 2011; Furusawa et al. 2013; Bennett et al. 2016) or when the host star is detected with high angular resolution adaptive optics (Batista et al. 2014, 2015) or Hubble Space Telescope (Bennett et al. 2006, 2007, 2015) observations, but in most cases, the mass and distance to the lens are not directly determined. Often, microlens planet discovery papers estimate the mass and distance with a Bayesian analysis, under the assumption that exoplanet properties do not depend on the host mass and distance. We would eventually like to measure this dependence on the host mass and distance, so we cannot use such mass estimates that assume the answer to a question that we are trying to investigate.

In this paper, we report on our analysis of the exoplanet mass-ratio function as a function of separation (in the angular Einstein radius units) using microlensing events detected by the Microlensing Observations in Astrophysics (MOA) collaboration alert system during the years 2007–2012. We use planet-to-host-star mass ratio instead of planet mass, because the former is measurable without any assumption regarding the estimate of host-star masses. This is the first statistical analysis of the prevalence and properties of exoplanets based on microlensing survey data, which includes full light-curve analysis of binary and planetary events to precisely determine the mass ratio and separation for events with significant planetary or binary signals. Previous analyses based on microlensing survey data either identified no definitive planetary events (Tsapras et al. 2003, 2016; Snodgrass et al. 2004), or included an incomplete light-curve analysis that did not definitively identify planetary events (Shvartzvald et al. 2016). Previous complete statistical analyses (Gould et al. 2010; Cassan et al. 2012) have been based on planets detected in the data from microlensing follow-up groups who focus on events thought to have relatively high sensitivities to exoplanets. In contrast, microlensing survey groups use large field-of-view (FOV) telescopes to observe a much larger number of events. Many of these events have low planet detection efficiency, but overall, the sensitivity of the surveys is higher because of the much larger number of events observed (Tsapras et al. 2016). This is why new microlensing projects just coming online now (Kim et al. 2016) and planned for the future (Spergel et al. 2015) are all high-cadence surveys that will find exoplanet signals in high-cadence, wide-field, survey data.

This paper is organized as follows. We describe the MOA survey observations, as well as observations by other telescopes that are used to help classify and characterize anomalous microlensing events in Section 2. In Section 3, we explain how microlensing events from the full 2007–2012 sample are selected for use to search for planetary signals and determine our planet detection efficiency. We describe our detection efficiency to planets in Section 4. The exoplanet mass-ratio function is presented in Section 5, and our results are discussed in Section 6. Finally, we summarize our results in Section 7.

2. OBSERVATIONS

The microlensing method requires the monitoring of tens of millions of stars to find gravitational microlensing events. It also requires high-cadence photometric monitoring to detect planetary signals, which have short timescales, ranging from a few hours to a few days for Earth- and Jupiter-mass planets, respectively. To find such rare, unpredictable planetary signals, the survey and follow-up teams have collaborated to comprise a large observational network.

Since 2006, the MOA collaboration has surveyed the Galactic bulge using the dedicated 1.8 m MOA-II telescope at Mt. John University Observatory in New Zealand. This telescope employs the MOA-Cam3 (Sako et al. 2008) CCD camera consisting of ten 2 k × 4 k detectors with a plate scale of 0.58 arcsec pixel−1. It is mounted at the prime focus and has a 2.2 deg2 FOV. This wide FOV allows us to monitor ∼44 deg2 of the Galactic bulge with cadences ranging from 15 minutes for the fields with the highest microlensing rate, to 45 minutes for fields with a medium microlensing rate, to 95 minutes for the fields with the lowest microlensing rate. These high-cadence observations allow the MOA survey to detect planetary signals in all the observed microlensing events, rather than just the fraction of events that are monitored for a fraction of their duration in the survey plus follow-up strategy suggested by Gould & Loeb (1992) and employed in the previous statistical analysis of Gould et al. (2010)16 and Cassan et al. (2012).17 The larger number of events that can be surveyed for planets leads to a higher planet detection rate. We find nearly four planetary events per year compared to an average of about one planetary event per year from the alert and follow-up strategy.

Once MOA identifies a candidate microlensing event in progress, we issue a microlensing event alert (Bond et al. 2001), often within 1 hr of the observation that triggered the alert. It is these events, found by the MOA alert system that we use for our sample of microlensing events. Our sample includes events announced by the MOA alert system in 2007–2012, and they consist of 488, 477, 563, 607, 485, and 680 events, respectively in the years 2007–2012. The MOA alert system was only partially functional in 2006, so the 2006 data are not included in this analysis.

The other major survey group that operated during the 2007–2012 seasons is the OGLE collaboration (Udalski et al. 2015b), which uses the 1.3 m Warsaw Telescope at the Las Campanas Observatory in Chile. The years 2007, 2008, and very early 2009 were the final years of the OGLE-III survey using a CCD camera with a FOV of 0.25 deg2, which was too small an FOV for a high-cadence survey, but in 2010, OGLE upgraded its CCD camera to the 1.4 deg2 OGLE-IV camera and began the high-cadence OGLE-IV survey.

A third high-cadence survey was started in 2011 at the Wise Observatory in Israel in order to fill the longitude gap between the MOA-II and OGLE-IV (Shvartzvald & Maoz 2012). The Wise survey uses a 1.0 deg2 CCD camera mounted on the Wise 1.0 m telescope. Because of the Wise Observatory's latitude of 30fdg60 N, the Wise survey can observe the bulge for no more than 5 hr per night. As a result, they have limited their high-cadence survey observations to 8 deg2 of the Galactic bulge.

In addition to these survey teams, most of the planetary and high-magnification events are observed by the microlensing follow-up groups, such as the Probing Lensing Anomalies NETwork (PLANET), the μlensing Follow-up Network (μFUN), RoboNet, and the Microlensing Network for the Detection of Small Terrestrial Exoplanets (MiNDSTEp). Originally, before the high-cadence surveys came online, these follow-up programs were needed to get high-cadence observations of the events thought to be most promising for planet detection, following the alert plus follow-up strategy laid out by Gould & Loeb (1992). This strategy led to the discovery of number of high-magnification planetary events (e.g., Gould et al. 2006; Dong et al. 2009a; Janczak et al. 2010; Miyake et al. 2011; Bachelet et al. 2012b; Yee et al. 2012), as well as a couple of planetary signals in low-magnification events (Beaulieu et al. 2006; Batista et al. 2011).

Follow-up observations were also used to help characterize planetary signals that were detected in progress in the MOA-II survey data (Sumi et al. 2010; Muraki et al. 2011; Furusawa et al. 2013; Skowron et al. 2015). It is important to treat such follow-up observations carefully in a statistical analysis. If these observations are triggered by the detection of a planetary signal, they cannot be used to increase the significance of the planetary signal above the planet detection threshold, because they were only taken because the planetary signal was already detected. The situation is different with the additional observations that are triggered by a high-magnification alert. Many high-magnification events receive intensive observations from follow-up groups because they have much more sensitivity to planetary signals than low-magnification events. The crucial difference is that the additional observations for a high-magnification events are triggered by this high sensitivity and not an actual planetary signal. So, very high-cadence observations of high-magnification events are often planned before any planetary signal is detected, so it is possible to use the follow-up observations to set detection thresholds for a statistical analysis (Gould et al. 2010), as long as any observations taken in response to the detection of a planetary signal are not included.

Most planetary signals are presently found by the high-cadence surveys, and MOA-II, as the first high-cadence survey, began finding planetary signals in survey data back in 2007 (or 2006 if non-alerted discoveries are included; Bennett et al. 2012). But, with only a single high-cadence observing site until 2010, and only partial coverage of the MOA high-cadence fields afterwards, we have found that planetary signals identified in the MOA-II can often be characterized much better when follow-up data, as well as OGLE and the WISE survey, are included in the light-curve modeling. Therefore, in our analysis, we use only the MOA-II survey data for detection of potential planetary light-curve anomalies, but we use all available data for characterization of the binary light-curve parameters.

3. EVENT SELECTION

In 2007–2012, a total of 3300 microlensing alerts18 were issued by MOA. Most of the alerts are microlensing events, but some other types of transient events are also included. To determine the exoplanet mass-ratio function, it is necessary to construct a sample of well-characterized single lens and planetary events with good coverage by MOA-II survey observations. While it is possible to detect planets in events with a strong stellar binary signal (Gould et al. 2014; Poleski et al. 2014), the analysis of such events is substantially more complicated, so we postpone such an effort to a later analysis. Therefore, events that are strong stellar binary events or are obviously not due to microlensing are removed from our sample. However, we do not rely upon simple light-curve inspection to distinguish between stellar binary and planetary-binary microlensing events. Instead, we have done a systematic analysis of the binary microlensing events in our sample to determine the best-fit mass ratio, q. After confirming that an event passes our binary or planetary detection criteria that the binary or planetary model should improve the fit χ2 by Δχ2 > 100 using the data on the MOA alert page, we collect re-reduced data from all the observatories that have observed each event. This data is not used for determining which events pass our selection criteria, but it is used to find the best-fit models that are used for event classification.

The modeling with these expanded data sets was performed with a grid search using the image-centered ray-shooting method (Bennett & Rhie 1996; Bennett 2010). As discussed by Suzuki et al. (2014), to find the best-fit model, we searched a wide range of the parameter space using a variation of the Markov Chain Monte Carlo (MCMC) algorithm (Verde et al. 2003), starting from 9680 grid points with fixed q, s, and α (s: the projected separation normalized to the angular Einstein radius, α: the angle between the source trajectory and the binary-lens axis). Then, we searched χ2 minima of the best 100 models in order of χ2 in the first grid search with letting the each parameter free. We also included the microlensing parallax effect (e.g., Muraki et al. 2011) to the model in this systematic analysis; 678 events were analyzed with this modeling effort. Independent modeling of all unpublished planetary events and possible planetary events was also performed using the initial condition grid method of Bennett (2010).

Following Bond et al. (2004), we define planetary events to be ones with a best-fit mass ratio of q < 0.03. We also check for degenerate solutions that have best-fit χ2 values that differ by Δχ2 < 25 (after the error rescaling as discussed below), and in these cases, we include all the models as possible models for the event. While a threshold of Δχ2 = 10 would be sufficient for Gaussian statistics, since it would imply a probability of <1%, we use a higher threshold in order to enable the possibility of testing our results for sensitivity to non-Gaussian effects. We consider all possible solutions in Bayesian analysis to determine the exoplanet mass-ratio function. For any binary events that might be potentially planetary events, we attempt to improve the MOA photometry with a variety of methods in order to get the best possible event classification, and, as mentioned in the last section, we include data from other groups. But, the anomaly detection is based on the data directly from the MOA alert pages.

The systematic modeling has recovered all of the previously known planetary events and detected four new planetary events. One of these four events was previously published as a stellar binary event, but our systematic analysis showed that it is clearly a planetary event. The other three planetary events were newly discovered by the systematic modeling effort, and their planetary signals were missed in the years in which they occurred. Our treatment of all these planetary events is discussed in Section 3.2. We find 23 planetary events (including one ambiguous event) and 1451 single-lens events in our six-year MOA-II survey data. This compares to six planets (in 5 events) and 8 single-lens events in the high-magnification sample of Gould et al. (2010) and three planets and 196 single-lens events in six years of PLANET follow-up photometry of events alerted by OGLE. One planetary event, OGLE-2007-BLG-349/MOA-2007-BLG-379, is common to all three analyses, but there is no other overlap in these planetary event samples.

Here, we briefly summarize our classification of alert events. We classify the events into six categories: single-lens, planetary, stellar binary, ambiguous, cataclysmic variable (CV) stars, and unknown events. Single-lens events have simple light curves, known as Paczyński (1986) curves, and we use this simple light-curve form to identify well-characterized single-lens events with the selection criteria given in Table 1 as described in the following subsection. For planetary (stellar binary) events, following three conditions are required: the MOA light curves are better fit by a binary-lens model than a single-lens model by Δχ2 > 100; the best-fit mass ratio is smaller (larger) than 0.03 using all available data including OGLE and other follow-up data, and the best-fit model is better than stellar binary (planetary) solutions by Δχ2 > 25. If the χ2 difference between the best planetary and stellar binary solution is smaller than Δχ2 = 25, then the event is classified as an ambiguous event. The light curves of CVs are not well fit by both single-lens and binary-lens model, and their light curves are so characteristic that they can be classified by visual inspection. The events classified as "other" are ones that fail the selection criteria (described in the next subsection) due to such problems as poor light-curve coverage or large systematic errors in the photometry. We avoid possible misclassification by setting our χ2 threshold to be Δχ2 > 100, which is much higher than would be necessary if we were only concerned with Gaussian random noise. There are at least two planetary candidates that have signals below this detection threshold, OGLE-2007-BLG-292/MOA-2007-BLG-190 and OGLE-2012-BLG-0724/MOA-2012-BLG-323 (Hirao et al. 2016). For the purpose of this paper, these are treated as non-detections, and they are the price we pay for setting a threshold high enough to avoid misclassifications.

Table 1.   Selection Criteria for the Well-monitored Events

Cut-1 2 < tE < 300 and u0 < 2
Cut-2 σu0/u0 < 0.40 or σu0 < 0.02, σtE/tE < 0.25 and σtE < 20
Cut-3 $\sum f/\mathrm{fe}\gt 500$ within $| t-{t}_{0}| \lt {t}_{{\rm{E}}}$
Cut-4 the number of data points, Ndata > 30 within $| t-{t}_{0}| \lt {t}_{{\rm{E}}}$
Cut-5 $\sum f/\mathrm{fe}\gt 200$ within $0\lt t-{t}_{0}\lt {t}_{{\rm{E}}}$ and $-{t}_{{\rm{E}}}\lt t-{t}_{0}\lt 0$
Cut-6 Ndata > 10 within $t-{t}_{0}\lt -20$ and $t-{t}_{0}\gt 20$
  Ndata > 3 within $| t-{t}_{0}| \lt {u}_{0}\times {t}_{{\rm{E}}}$
Cut-7 Ndata > 20 within $| t-{t}_{0}| \lt 0.1\times {t}_{{\rm{E}}}$
  $\sum f/\mathrm{fe}\gt 2000$ within $| t-{t}_{0}| \lt 0.2\times {t}_{{\rm{E}}}$

Note. The single-lens model used for the event selection is described by three parameters: the time of the peak magnification, t0, the Einstein radius crossing time (or event timescale), tE, and the impact parameter in units of the angular Einstein radius, u0. The parameters with σ in the table are 1σ uncertainties for the fit parameters. The parameter f and fe are the delta-flux (see the text) and its error bar.

Download table as:  ASCIITypeset image

3.1. Selection Criteria for Single-lens Events

Once the planetary, stellar binary events and CVs are removed from the sample, the remaining events are the events classified as single-lens or "other" events. In order to consider only those events that allow the detection of planets with well-defined parameters, we apply a set of cuts, listed in Table 1 on the remaining events. These cuts rely upon single-lens microlensing fits to the online MOA-II data. In those cases where a finite-source microlensing model provides a better fit than a point-source model, we use the finite-source model for cuts given in Table 1. Note that none of these criteria are related to the presence or absence of planetary light-curve signals.

In order to search for planetary signals in the microlensing light curves, we require well-constrained single-lens microlensing fit parameters. Otherwise, we cannot tell what types of planets the light curve might be sensitive to. Cut-1 and Cut-2 impose these constraints and also exclude events that may not be due to lensing by stellar or brown dwarf primaries. The flux, f, in the MOA online data is the delta-flux, which is the measured flux in the difference image, as computed by the difference image method (Bond et al. 2001). The flux error, fe, is the error bar computed by the photometry code, so Cut-3 is a measure of the S/N of the primary microlensing signal above the baseline brightness. Cut-4 requires a sufficient number of data points in the magnified part of the light curve to detect a planetary signal, and Cut-5 ensures that there is reasonable coverage and a reasonable S/N on both the rising and falling sides of the light curve. We also want to avoid microlensing events that reach their single-lens peak outside of or too close to the end of the season, and this is the rationale for Cut-6. Cut-7 consists of three conditions that must be satisfied for the event to be accepted. All three conditions ensure that the event does not have large gaps in data coverage around the peak. The events failing Cut-7 also have very low planet detection efficiency, so this cut has little effect on the total detection efficiency.

The result of applying these cuts on our sample of 3300 microlensing events is that 1451 events pass all of the cuts. These 1451 events consist of 226, 230, 247, 288, 184, and 276 events, respectively, in each year of the range 2007–2012. Since the selection criteria we use are not related to the presence or absence of planetary signals, the rejection of poorly sampled events induces no bias. Our selection criteria do influence the properties of the distribution of stars that are searched for planetary signals, but this distribution is already influenced by unavoidable factors, such as observing conditions. Some readers might be surprised at that only half of the observed events passed the event selection criteria, but this ratio is comparable to previous studies. Cassan et al. (2012) and Shvartzvald et al. (2016) used 45% and 49% of observed events for their analysis. Gould et al. (2010) used only 4% of the observed events due to the strict cuts they used.

The cumulative tE distribution for these 1451 events is shown in the top panel of Figure 1, along with the distribution of events with planets. As we show in Section 5.2, the difference between these distributions is largely, but perhaps not entirely, a selection effect. The lower panel of Figure 1 shows these same distributions for the μFUN (Gould et al. 2010) and PLANET (Cassan et al. 2012) statistical analyses in red and blue, respectively. Both show a similar bias toward planetary signals in long events, while the μFUN distribution shows a slight bias toward longer events in the full sample, including single-lens events. The full tE distribution for our sample is similar to that of the PLANET sample (Cassan et al. 2012), which is only slightly skewed toward longer duration events, compared to our sample.

Figure 1.

Figure 1.  Cumulative distribution of the event timescale tE, for the 1474 single-lens (thin lines) and planetary events (thick lines) that pass cuts 1–7. The top panel compares our full sample to the planetary events in the sample. The bottom panel shows the same tE distributions for the previous statistical samples of Gould et al. (2010) and Cassan et al. (2012), respectively.

Standard image High-resolution image

3.2. Planetary Events

The systematic modeling of all the anomalous events in our sample was conducted using all available data, including data from the OGLE survey and follow-up groups. This effort found four new planetary events, as well as the known planetary events that have either been published or are in preparation for publication. In order to determine the significance of the planetary signals in the MOA data, all of the planetary events were fit with both single-lens and planetary models using only the online MOA data. For possible planetary signals discovered while they are in progress, MOA often takes additional data in order to better characterize the parameters of candidate planetary microlensing events. We must be careful not to let these additional data affect the planetary detection statistics, because data taken in response to the detection of the planetary signal cannot be used to determine the detection efficiency. So, the additional data taken in response to the anomaly detection are removed from the fits used to determine which events pass the anomaly detection threshold. The detection threshold fits use only the data taken at the standard MOA survey cadence. We calculate the fit χ2 for the best-fit single-lens and binary-lens models, and we renormalize the error bars to set χ2/dof ≡ 1 for the best-fit model. This process is necessary before using Δχ2 values as detection thresholds because the error bars provided by crowded field photometry codes are more accurate measures of the relative rather than absolute photometric accuracy. This procedure is commonly used in most of the previous microlensing analyses (e.g., Muraki et al. 2011). Then, we calculate the fit χ2 for the best-fit single-lens and planetary-binary-lens models again with these renormalized error bars. If the difference between the renormalized Δχ2 values for the best-fit single-lens and the best-fit planetary models is Δχ2 > 100, then the event passes our detection threshold. There are four planetary events in the MOA data set that do not pass this detection threshold cut because the planetary signals for these events were observed in other data sets, so that a significant planetary signal is not seen in the MOA data. These events are MOA-2007-BLG-400 (Dong et al. 2009a), MOA-2008-BLG-310 (Janczak et al. 2010), MOA-2011-BLG-293 (Yee et al. 2012), and OGLE-2012-BLG-0724/MOA-2012-BLG-323 (Hirao et al. 2016), and they are examples of lens systems with planets that are not detected because of the low planet detection efficiency of the MOA data for these events.

We require that our sample of planetary events should have that same good light-curve coverage that we require for non-planetary (single-lens) events. For each possible planetary event, an artificial single-lens light curve is generated using the parameters of the best-fit planetary model and the Gaussian noise based on the residuals of the best-fit model. Then, our selection criteria (see Section 3.1 and Table 1) are applied to each artificial single-lens light curve for each candidate planetary event. We find that 23 events that are best fit by a binary-lens model with a planetary mass ratio pass our cuts. However, for one of these events, there is a competing model with a stellar binary mass ratio that fits the data nearly as well. We classify this event as an "ambiguous event", which we treat in detail in Section 3.3. This leaves the 22 clear planetary events, which we list in Table 2. It is the Einstein radius crossing times, tE, of these 22 clear and 1 ambiguous planetary events that are plotted in Figure 1. The median tE value for the planetary sample is 42 days, compared to 25 days for the full sample. Thus, the typical planetary event is 1.7 times longer than the typical well-sampled event. This is primarily due to selection effects, as we discuss in Section 5.2.

Table 2.   List of Planetary Events

Event Name tE u0 q s C.C. OGLE ID Reference
  days θE ×10−3 θE
MOA-2007-BLG-192 71.3 0.004 0.179 1.00 ? (a)
MOA-2007-BLG-308 55.4 0.079 0.095 0.926 yes 368 (b)
MOA-2007-BLG-379 121 0.002 0.316 1.26 yes 349 (c)
MOA-2008-BLG-288 34.0 0.27 11.80 0.877 yes 355 (d)
MOA-2008-BLG-379 42.5 0.0060 6.85 0.903 yes 570 (e)
  42.1 0.0060 6.99 1.119 yes
MOA-2009-BLG-266 61.45 0.13 0.058 0.91 yes (f)
MOA-2009-BLG-319 16.6 0.0062 0.395 0.975 yes (g)
MOA-2009-BLG-387 40.1 0.089 13.35 0.913 yes (h)
MOA-2010-BLG-117 116.3 0.279 0.75 0.86 yes (i)
MOA-2010-BLG-328 62.6 0.072 0.26 1.15 yes (j)
MOA-2010-BLG-353 11.1 0.5 1.38 1.46 no (k)
MOA-2010-BLG-477 46.7 0.0034 2.18 1.12 yes (l)
MOA-2011-BLG-028 34.4 0.908 0.181 1.673 yes 0203 (m)
MOA-2011-BLG-197 53.63 0.130 3.954 1.04 yes 0265 (n)
MOA-2011-BLG-262 3.85 0.015 0.454 1.01 yes 0703 (o)
MOA-2011-BLG-291 23.6 0.004 0.409 1.21 yes (p)
MOA-2011-BLG-322 23.17 0.046 28.4 1.822 no 1127 (q)
MOA-2012-BLG-006 21.13 1.3 16.14 4.32 yes 0022 (r)
MOA-2012-BLG-288 77.5 0.0014 1.09 0.41 no 0563 (s)
  77.7 0.0014 1.09 2.43 no
MOA-2012-BLG-355 22.2 0.19 0.70 0.76 yes (t)
MOA-2012-BLG-505 9.9 0.0076 0.19 0.91 yes (u)
  10.0 0.0075 0.20 1.12 yes
MOA-2012-BLG-527 67.0 0.103 0.24 0.890 no 0950 (v)
  66.5 0.104 0.22 1.004 no

Note. The C.C. column shows caustic crossing. The OGLE ID column shows the last part of OGLE event name, i.e., OGLE-2007-BLG-368 for 368.

References. (a) Bennett et al. (2008), (b) Sumi et al. (2010), (c) Bennett et al. (2016), (d) Koshimoto et al. (2014), (e) Suzuki et al. (2014), (f) Muraki et al. (2011), (g) Miyake et al. (2011), (h) Batista et al. (2011), (i) D. P. Bennett et al. (2016, in preparation), (j) Furusawa et al. (2013), (k) Rattenbury et al. (2015), (l) Bachelet et al. (2012b), (m) Skowron et al. (2016), (n) Skowron et al. (2015), (o) Bennett et al. (2014), (p) D. P. Bennett et al. (2017, in preparation), (q) Shvartzvald et al. (2014), (r) R. P. Poleski et al. (2016, in preparation); (s) Fukui et al. (2015), (t) OGLE et al. (2016, in preparation); (u) M. N. Nagakane et al. (2016, in preparation), (v) Koshimoto et al. (2016).

Download table as:  ASCIITypeset image

Our sample of 22 clear planetary events includes four that were found for the first time by our systematic modeling effort. These events are OGLE-2008-BLG-355/MOA-2008-BLG-288 (Koshimoto et al. 2014), MOA-2008-BLG-379/OGLE-2008-BLG-570 (Suzuki et al. 2014), MOA-2010-BLG-353 (Rattenbury et al. 2015), and MOA-2011-BLG-291 (D. P. Bennett et al. 2017, in preparation). Here are short descriptions of these four events.

MOA-2008-BLG-288: This event was initially announced as a stellar binary event in an analysis of OGLE-III archival data using only OGLE survey data (Jaroszyński et al. 2010). But, our systematic modeling of all the anomalous events using all available data found that a planetary model is strongly favored for this event, by Δχ2 = 810. The MOA data cover a caustic exit that is critical for the correct classification of this event.

MOA-2008-BLG-379: This is a high-magnification event with a faint source star, Is = 21.3. The faintness of the source meant that the strong anomaly at the peak dominated the apparent brightness variation, and it was not immediately recognized as a high-magnification event. As a result, the planetary nature of the event was not recognized in real time. The planetary feature is characterized by the survey data alone.

MOA-2010-BLG-353: This event has a relatively weak anomaly caused by the approach of the source to a cusp of a major image planetary caustic. The peak magnification is low and the source star is not so bright, and, as a result, the planetary anomaly for this event has the smallest significance of all planetary signals that pass our threshold. The OGLE data does not cover the anomaly, but it helps to constrain the non-planetary light-curve parameters.

MOA-2011-BLG-291: This is a high-magnification event that was alerted and observed by the μFUN follow-up group, as well as the MOA, OGLE, and Wise survey telescopes. Nevertheless, the planetary signal occurs almost entirely in the MOA and OGLE data sets. This event was identified as a possible planetary event shortly after the peak, but it was not seriously investigated prior to our systematic analysis.

High-magnification events (Amax ≳ 50) account for 41% of our clear planetary event sample, compared to two of the three (67%) planetary events in the PLANET sample (Cassan et al. 2012) and all five of the planetary events (100%) in the μFUN sample (Gould et al. 2010) that are high-magnification events. (But note that there is one high-magnification event common to all three samples.) It is natural that analyses based on follow-up data will be dominated by high-magnification events, because the planet detection efficiency per event is much larger for high-magnification events than for low-magnification events (Griest & Safizadeh 1998; Rhie et al. 2000). For the microlensing follow-up results that have been published to date, the total detection efficiency in the sample (total sensitivity) has been dominated by the high-magnification events.

High-magnification events are relatively rare and the planetary caustics are usually larger than the central caustics, so signals from planetary caustics in low-magnification events are expected to dominate the planetary discoveries made in survey mode, and this is what we find with our sample. Most high-magnification planetary events suffer from close/wide separation degeneracy (Griest & Safizadeh 1998), because the shape of central caustic that depends on the mass ratio, q and the absolute value of the logarithm of the star–planet separation, $| \mathrm{log}\,s| $ (Chung et al. 2005). So, the central caustics are nearly identical for separations of s and 1/s. There is a similar degeneracy with the location of the caustics for low-magnification events, but the shapes of the caustics for s < 1 and s > 1 are very different. Thus, observations of low-magnification events can generally break this $s\leftrightarrow 1/s$ degeneracy, and the separation is uniquely determined. This means that low-magnification planetary events are suitable for study of planet distribution as a function of separation beyond the snow line. More planetary signals in low-magnification events will be expected in the future survey observation.

The planetary signals that we detect in microlensing light curves are due to the source star crossing or approaching a caustic or cusp. The caustic crossing events are of special interest because caustic crossing features in the light curves allow us to measure the angular Einstein radius, θE, as mentioned in Section 1. We find that 23% of our planetary events do not include measured caustic crossings. This percentage is larger than the ∼12% of published microlensing planet discoveries without measured caustic crossings, but it is substantially smaller than the ∼55% of the planet detections predicted by simulation (Zhu et al. 2014) of a three-site high-cadence survey like the Korean Microlensing Telescope Network (Park et al. 2012, KMTNet). This simulation makes a number of optimistic assumptions, but it illustrates the point that many of the planetary signals detected by higher sensitivity surveys will be in non-caustic crossing signals from higher mass-ratio planets (q ≳ 10−3). This is borne out by the non-caustic-crossing planetary deviations in our sample. But MOA-2012-BLG-527 (Koshimoto et al. 2016) is different from the other non-caustic-crossing planetary signals, because the source star crossed the demagnified region between the central caustic and two minor image, triangular planetary caustics. This region provides a strong negative perturbation that makes it easier to detect low-mass-ratio (q ∼ 10−4) planets.

In these non-caustic-crossing events, the detection of finite-source effects is difficult. So, the source radius crossing time is not well constrained for MOA-2010-BLG-353 (with a 2.3σ signal), MOA-2011-BLG-322 (Shvartzvald et al. 2014), and MOA-2012-BLG-527. However, the mass ratios are measured precisely in these events. This means that we must assume the source radius crossing times when calculating the detection efficiency for these events, but this assumption has no significant effect on the results we present in this paper.

In our planet sample, five planetary events have not been published yet: MOA-2010-BLG-117 (D. P. Bennett et al. 2016, in preparation), MOA-2011-BLG-291 (which is described above), MOA-2012-BLG-006/OGLE-2012-BLG-0022 (R. P. Poleski et al. 2016, in preparation), MOA-2012-BLG-355 (OGLE et al. 2016, in preparation), and MOA-2012-BLG-505 (M. N. Nagakane et al. 2016, in preparation). Here, we briefly describe each event.

MOA-2010-BLG-117: This is a low-magnification event with a binary-source star whose trajectories crossed almost perpendicular to the minor image caustics. The best-fit planet model that assumes a single-source star failed to explain the demagnified part between the two triangular caustics. The observed flux at the trough was brighter than the model. This is well fit by adding a bright companion to the source, and the binary-source model fits the data much better than a single-source, triple-lens model, which could also explain the excess flux at the time of the trough. This event was also observed by the OGLE-IV survey, but OGLE did not announce any alerts in the first year of the OGLE-IV survey.

MOA-2012-BLG-006: This event was detected as an apparent short timescale event, but a few month later the magnification rose again with a much longer duration, but a lower amplitude. Light-curve modeling (Poleski et al. 2016, in preparation) indicates that the initial magnification was caused by a planet located at s ∼ 4, which is the largest separation among the discovered microlensing planets, except for OGLE-2008-BLG-092LAc (Poleski et al. 2014).

MOA-2012-BLG-355: This is a low-magnification planetary event, but an RR Lyrae star is located very close to this event on the sky. It contaminates the light curve, which slowed the identification of the planetary signal by a couple days. The proximity to a known variable prevented its detection in the OGLE EWS system. We removed the contamination from the MOA and OGLE data, and the best-fit model was found as shown in Table 2.

MOA-2012-BLG-505: This event is a short, high-magnification event that was only observed by MOA. An anomaly lasting only a few hours over the light-curve peak was well covered due to MOA's high-cadence observing strategy. As often occurs with high-magnification events, there is a close-wide degeneracy, but it does not imply a large uncertainty in s because both values are close to 1.

False positives are generally not a problem for planetary microlensing signals, unlike the cases of planetary transit (Mullally et al. 2016) and RV (Robertson et al. 2014) signals. The only types of false positive that have been identified for planetary microlensing events are binary-source stars (Gaudi 1998) and variable-source stars microlensed at high magnification (Gould et al. 2013). Both of these types of false positives are intrinsically rare because a binary companion to the source must be microlensed at high magnification to mimic a planetary microlensing event and because variable-source stars are themselves intrinsically rare. Furthermore, most planetary microlensing light curves have shapes that are not compatible with false positives. In our sample, only MOA-2010-BLG-353 has a light curve that might be compatible with a false positive. It is a low-magnification event, so source variability with an amplitude large enough to mimic the microlensing signal would easily be seen in the baseline data, but a binary source might be a possibility. We fit a binary-source model to this event and find a best-fit χ2 that it is worse than the best-fit planetary model by Δχ2 = 20.3. The error bars during the MOA-2010-BLG-353 planetary signal are ≳4%, which is well above the usual level of the correlated errors that are sometimes seen in microlensing light curves. Thus, it is safe to exclude the binary-source model for this event.

Errors in light-curve modeling might also be considered to be false positives, but these are not a problem for our analysis because of our systematic modeling procedure. All anomalous events have been systematically modeled by at least two independent investigations, and potential planetary events have been modeled by at least three independent efforts. The recently published OGLE event, OGLE-2013-BLG-0723, provides a good example of the effectiveness of the MOA light-curve modeling program. The initial paper (Udalski et al. 2015a) claimed the discovery of a Venus-mass planet orbiting a brown dwarf in a binary system. After this paper was published, the photometric data for this event was obtained from the authors, and one of us (D.P.B.) began systematic modeling of the event, using the procedure described in Bennett (2010). This systematic modeling effort quickly revealed the correct model, which includes only a stellar binary with no planet and a much smaller microlensing parallax signal. The authors of the original paper were contacted, and this resulted in a paper with the correct model (Han et al. 2016).

Photometry for published events should be available permanently on the NASA Exoplanet Archive (http://exoplanetarchive.ipac.caltech.edu/), while the data for all planetary events in this sample, including those that have yet to be published can be obtained directly from the first or second author.

3.3. Ambiguous Events

We define the mass-ratio threshold q < 0.03 as the upper limit for planetary events. Binary-lens systems with larger q values are considered to be stellar binary events. However, for a number of events, there are several models that can fit the data. In most cases, these degenerate models are caused by a well-known degeneracy between different planetary models, such as the $s\leftrightarrow 1/s$ degeneracy discussed in Section 3.2, but there can also be degeneracies between planetary and stellar binaries, as discussed by Choi et al. (2012a). In fact, one of the events analyzed by Choi et al. (2012a), OGLE-2011-BLG-0950/MOA-2011-BLG-336, is in our sample. We cannot arbitrarily ignore or accept such an event, because either choice would bias our measurement of the exoplanet mass-ratio function. So, we include both possibilities in a Bayesian analysis, as discussed below.

Our systematic analysis of all binary events turned up a number of other potentially ambiguous events, but further light-curve analysis with optimized MOA photometry revealed that most of these were unambiguous stellar binary events or contaminated by systematic photometry errors. Events were classified as ambiguous based upon the χ2 values of the best-fit stellar binary and planetary models (${\chi }_{\mathrm{SB}}^{2}$ and ${\chi }_{\mathrm{PL}}^{2}$, respectively). These χ2 values are calculated with renormalized error bars to give χ2/dof = 1, and we classify an event as ambiguous if $| {\chi }_{\mathrm{SB}}^{2}-{\chi }_{\mathrm{PL}}^{2}| \lt 25$. With this definition, the only ambiguous event is the one identified by Choi et al. (2012a), OGLE-2011-BLG-0950/MOA-2011-BLG-336. The other possible ambiguous events did not pass our event selection cut. Choi et al. (2012a) found that the planetary model was favored by Δχ2 = 105, but considered the event to be ambiguous because of evidence for systematic photometry errors. Our optimized photometry has removed much of these systematic photometry problems, but this also reduced the χ2 difference between the best stellar binary and planetary models to Δχ2 = 19.80, favoring the planetary model. This is the difference between the best-fit stellar binary and planetary models, but there are two models in each category, due to the close-wide degeneracy for high-magnification events. The wide model is favored in each case, as can be seen in Table 3, which lists the light-curve model parameters and fit χ2 values for each of the four (approximately) degenerate for this event.

Table 3.   Parameters of the Ambiguous Event

Event Name Model t0 tE u0 q s α t* χ2 C.C. OGLE ID
    HJD' days 10−3 10−3   rad days
MOA-2011-BLG-336 PL(close) 5786.3980 −69.62 8.05 0.611 0.688 4.741 0.296 8176.29 no 0950
  σ 0.0006 1.47 0.17 0.038 0.011 0.003 0.006
  PL(wide) 5786.3974 70.70 −7.95 0.551 1.409 1.600 0.301 8175.21 no
  σ 0.0007 1.83 0.21 0.035 0.021 0.003 0.005
  SB(close) 5786.3975 −65.53 −8.45 471.197 0.075 3.874 0.052 8196.07 no
  σ 0.0011 1.21 0.18 141.309 0.001 0.007 0.069
  SB(wide) 5786.3940 −123.21 −4.49 2548.603 24.131 3.875 0.052 8195.01 no
  σ 0.0007 7.70 0.32 429.180 0.746 0.004 0.056

Note. The C.C. column shows caustic crossing. The first lines show the parameters of each model, and the second lines show the uncertainties. PL is planetary model, and SB is stellar binary model. The OGLE ID column shows the last part of OGLE event name, i.e., OGLE-2011-BLG-0950.

Download table as:  ASCIITypeset image

The probability that this ambiguous event is due to a planetary system rather than a stellar binary depends both on the χ2 values for the respective fits and the prior probabilities of lens systems with the parameters of the best-fit models. We define the relative fit probability of the stellar binary model to be ${e}^{-{\rm{\Delta }}{\chi }^{2}/T}$, where T = 2 is the correct value if the observational uncertainties follow an uncorrelated Gaussian distribution. When T = 2, this is equivalent to the likelihood ratio that is sometimes used for similar model comparisons (Bozza et al. 2016). We can adjust T in order to test the effects of non-Gaussian uncertainties.

In order to determine the prior probability for the stellar binary model, we need to estimate the mass of the primary star, and we do this with a Bayesian analysis using constraints from the observed tE values. We conduct the analysis with the same manner as used in the previous planetary events (e.g., Beaulieu et al. 2006). The Galactic model of Han & Gould (2003) was used for the prior. This yields a median primary star mass of ∼0.6 M. The period and mass-ratio distribution for stellar binaries have been determined for solar-type stars by Raghavan et al. (2010) and for low-mass stars by Delfosse et al. (2004). Since the likely primary mass for the OGLE-2011-BLG-0950/MOA-2011-BLG-336 lens system is intermediate between the solar-type and low-mass stars, we use the average value of the probability for solar-type and low-mass stars. For solar-type stars, the stellar multiplicity is ∼41%, while for low-mass stars it is ∼26%. For the probability of a planetary companion with the parameters indicated by each fit, we use our mass-ratio function. The error bars in Table 3 are also important, as they serve to determine the range in the mass ratio, q, and separation, s, that are consistent with the light-curve data, and these factors tend to enhance the probability of the stellar binary interpretation, because the stellar binary models are consistent with a wider range of separations and mass ratios than the planetary models. The priors that we derived from these considerations were 0.196294 and 0.072682 for the wide and close planetary models and 0.530297 and 0.200727 for the wide and close stellar binary models, respectively.

In order to obtain the final probabilities for the different models, we must multiply these priors by the relative probabilities, ${e}^{-{\rm{\Delta }}{\chi }^{2}/T}$, for the different models. With the standard, T = 2, value for uncorrelated Gaussian uncertainties, we find final probabilities of 0.825276 and 0.174628 for the wide and close planetary models and 0.000078 and 0.000017 for the wide and close stellar binary models, respectively. These values are used for our final results, but we have also considered the possibility of strong correlations or non-Gaussianity in the photometric uncertainties. If we increase T to T = 20, ten times the Gaussian random value, then the relative probabilities become 0.375144 and 0.131346 for the wide and close planetary models and 0.363331 and 0.130179 for the wide and close stellar binary models, respectively. So, this would imply roughly equal probabilities for the planetary and stellar binary models. Fortunately, even such an extreme assumption of correlated or non-Gaussian errors, does not have a significant effect on our final results.

4. PLANET DETECTION EFFICIENCY

The detection efficiency for a planetary signal in a microlensing event is denoted by $\epsilon (\mathrm{log}s,\mathrm{log}q)$. This is the probability that a planet with its mass ratio, q, at the projected separation (in units of the Einstein radius), s, is detected, assuming a random orientation for the planetary system. Following Rhie et al. (2000), we use the best-fit single-lens model for each event, while for events with planetary signals we generate an artificial single-lens light curve using the parameters of the best-fit planetary model. We inject planetary signals with a complete range of parameters ($\mathrm{log}s,\mathrm{log}q,\alpha $) covering all of the potentially detectable signals to each single-lens light curve. The generated light curves have same observation times as actual MOA survey data. The error bars in each data point are renormalized to be χ2/dof ≡ 1 for the best single-lens or planetary model, as discussed in Section 3.2. Note that any additional MOA observations taken in response to a planetary signal are removed at this stage, as described in Section 3.2. For each of these generated light curves, we search for the best single-lens model, we then calculate ${\rm{\Delta }}{\chi }^{2}={\chi }_{\mathrm{single}}^{2}-{\chi }_{\mathrm{planetary}}^{2}$ to determine if the planetary parameters give a large enough signal to pass our detection threshold. We do not need to calculate ${\chi }_{\mathrm{planetary}}^{2}$ because it remains the same as the original χ2 for single-lens fit without a binary signature added. We use 360 grid points to cover the angle between the source trajectory and lens axis, α, and define $\epsilon (\mathrm{log}s,\mathrm{log}q)$ as the fraction of grid points for which ${\rm{\Delta }}{\chi }^{2}(\mathrm{log}s,\mathrm{log}q,\alpha )\gt {\chi }_{\mathrm{thresh}}^{2}$. Of course, we must use the same ${\chi }_{\mathrm{thresh}}^{2}$ as we used as the detection threshold for the planetary signals, discussed in Section 3.2. So, we set ${\chi }_{\mathrm{thresh}}^{2}=100$. Microlensing is insensitive to inner planets with s < 0.1 and has a very low sensitivity for outer planets with s > 10, except in events where the microlensing signal of the host star is not directly detected (Bennett et al. 2012). So, for the purpose of the detection efficiency calculations, we sample s with 40 logarithmically spaced bins spanning the range 0.1 ≤ s ≤ 10. Since we know that the detection efficiency depends smoothly on q, we use 9 logarithmically spaced bins covering the range $-5\leqslant \mathrm{log}q\leqslant -1.5$. Thus, the total number of artificial light curves simulated for one event is 129,600. For the handful of events sensitive to planets with smaller mass ratios, we extend the grid down to $\mathrm{log}q=-6$.

Formally, if we assume that the photometric errors follow an uncorrelated Gaussian distribution, then our detection threshold of ${\chi }_{\mathrm{thresh}}^{2}=100$ corresponds to a false alarm probability of ∼2 × 10−23. However, we do not believe that our error bars follow an uncorrelated Gaussian distribution. Also, we want to ensure that the planetary signal is detected with enough S/N to provide reliable determinations of the planetary light-curve parameters, q and s. Note that the analysis of high-magnification events by Gould et al. (2010) used a much higher threshold, ${\chi }_{\mathrm{thres}}^{2}=500$. However, this threshold was applied to data that was reprocessed after a planetary signal was suspected. Since some of the original data was reduced by DoPHOT (Schechter et al. 1993) instead of more precise difference imaging photometry, the reprocessed data is likely to have a substantially larger Δχ2 than the original data. Also, much of the data in the Gould et al. (2010) sample was taken with unfiltered telescopes. This means that these data are vulnerable to systematic photometry errors due to the color dependence of atmospheric transparency (Dong et al. 2009b; Bennett et al. 2010). Finally, the high-magnification events are often observed continuously all night, and this means that the effect of correlated photometric errors are much more significant than with the MOA data, with a cadence of at most one observation every 15 minutes. For these reasons, Gould et al. (2010) thought that their higher threshold was reasonable for their analysis.

4.1. Source Star Radius

The size of a source star is important for the detection of planets, because planetary light-curve features are broadened and depressed with increasing source star size. The finite-source effect is characterized by the source radius crossing time, t*, which is the time needed for the lens to cross the angular source star radius, θ*. A larger source star has more chance to cross caustics, but it usually weakens the deviation from the single-lens light curve. For the detection of planets with giant source stars and q ≲ 10−4 or main sequence sources and q ≲ 10−5, the source star size is critical because the planet signals can be smoothed out by finite-source effects (Bennett & Rhie 1996). Thus, finite-source effects can have important effects on the detection efficiency, so we need to estimate t* for each event.

For the planetary event, we use the t* value from each paper or the best-fit planetary model. Although the t* is not well constrained for MOA-2010-BLG-353, MOA-2011-BLG-322, and MOA-2012-BLG-527, we used the t* of the best-fit model. For the high-magnification events with lens stars passing over source stars MOA-2007-BLG-176, MOA-2007-BLG-233/OGLE-2007-BLG-302, OGLE-2008-BLG-290/MOA-2008-BLG-241, MOA-2010-BLG-436, MOA-2011-BLG-093, and OGLE-2011-BLG-1101/MOA-2011-BLG-325, we use the t* value for each event from Choi et al. (2012b). Also, we use the t* of Batista et al. (2009), Gould et al. (2009, 2013), Yee et al. (2009), Bachelet et al. (2012a), Yee et al. (2013), Freeman et al. (2015), and Park et al. (2014) for OGLE-2007-BLG-050/MOA-2007-BLG-103, OGLE-2007-BLG-224/MOA-2007-BLG-163, OGLE-2008-BLG-279/MOA-2008-BLG-225, MOA-2009-BLG-411, MOA-2010-BLG-311, MOA-2010-BLG-523, MOA-2011-BLG-274, and OGLE-2012-BLG-0455/MOA-2012-BLG-206, respectively. Note that MOA-2009-BLG-411 and MOA-2012-BLG-206 are stellar binary events with a brown dwarf companion or possible planetary events, but the anomaly signals are found in the follow-up data, instead of MOA data. So, they are treated as single-lens events in this analysis. In addition to the above events, we use the best-fit t* value estimated with all available data for the following events because they are probably showing the finite-source effect: MOA-2008-BLG-383, MOA-2011-BLG-393, MOA-2012-BLG-278/OGLE-2012-BLG-1430, OGLE-2012-BLG-0617/MOA-2012-BLG-285, OGLE-2012-BLG-1098/MOA-2012-BLG472, and OGLE-2012-BLG-1274/MOA-2012-BLG-538.

For the other single-lens events, we need to estimate t* by estimating both θ* and the relative lens-source proper motion, μrel, since t* = θ*/μrel. In most planetary microlensing events, the light curves are observed with different filters to estimate θ* using the empirical relation of the θ*, color and magnitude of the source star (Kervella et al. 2004; Kervella & Fouqué 2008; Boyajian et al. 2014). For the 2007–2011 observing seasons, we do not have source color information for the single-lens events, because the only filter used for regular MOA observations was the MOA-Red wide-band filter. But, the source star magnitude in the MOA-Red passband, RM,S, is determined by the light-curve model for each event. The absolute magnitude of the source star, MI,S, can then be estimated with the formula ${M}_{I,{\rm{S}}}={R}_{{\rm{M}},{\rm{S}}}-({R}_{{\rm{M}},\mathrm{RC}}-{M}_{I,\mathrm{RC}})$, where ${R}_{{\rm{M}},\mathrm{RC}}$ is the apparent magnitude of red clump giant (RCG) centroid in the MOA-Red passband, and ${M}_{I,\mathrm{RC}}$ is the absolute magnitude of RCG centroid in the I-band. ${R}_{{\rm{M}},\mathrm{RC}}$ is estimated from the color–magnitude diagram of the reference images in the RM and V-bands in each subframe (1k × 1k section of each CCD). We assume ${M}_{I,\mathrm{RC}}=-0.12$ following Nataf et al. 2013. Using the calculated ${M}_{I,{\rm{S}}}$ value with the 10 Gyr isochrone model from the PARSEC isochrones, version 1.2S (Bressan et al. 2012; Chen et al. 2014, 2015; Tang et al. 2014), the source star radius can be estimated. The source star distance depends on Galactic longitude, and we use the values from (Nataf et al. 2013, Table 1) to determine θ*.

We need to estimate the lens-source relative proper motion, μrel to determine t*. μrel ranges over about an order of magnitude, with $1\,\mathrm{mas}\,{\mathrm{yr}}^{-1}\lesssim {\mu }_{\mathrm{rel}}\lesssim 10\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$. As discussed in (Bennett et al. 2014), for events with both the source and lens located in the bulge, the typical value is $\langle {\mu }_{\mathrm{rel}}\rangle \simeq 5.6\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$. For bulge-source disk-lens events, we have $\langle {\mu }_{\mathrm{rel}}\rangle \simeq 8.2\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$. We assume the typical value for our microlensing events is $\langle {\mu }_{\mathrm{rel}}\rangle \simeq 5.6\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$, because the typical lens star is located at the bulge. Of course, this is an underestimate of μrel for cases when the lens star is located at the disk, which means it overestimates t* for these events. But, these overestimates are by only about 50%, and this has a very small effect on our final results.

To check whether our results are sensitive to our treatment of finite-source effects, we consider the uncertainty in t*, which is a quadrature sum of the uncertainties in θ* and μrel. We assume a 30% error for the estimated θ* values. Our planetary event sample yields an average relative proper motion of ${\mu }_{\mathrm{rel}}={5.74}_{-2.84}^{+2.94}$ mas yr−1, and we assume that the single-lens events have the same μrel distribution as the planetary events. In Figure 2, we compare the estimated uncertainties in t* to the measured t* values for planetary events and finite-source events. The dimensionless source star size is given by $\rho ={\theta }_{* }/{\theta }_{{\rm{E}}}={t}_{* }/{t}_{{\rm{E}}}$, and its uncertainty, σρ, for each event is determined from the uncertainties in θ*, μrel, and tE. The detection efficiency will be computed in the next section with three different values of the dimensionless source size: our estimated value, ρ, as well as ρminand ρmax, where ${\rho }_{\min }=\rho -{\sigma }_{\rho }$ and ${\rho }_{\max }=\rho +{\sigma }_{\rho }$.

Figure 2.

Figure 2.  Source radius crossing time t* vs. the relative magnitude of source stars to red clump giant (RCG) centroid in MOA-Red filter. The gray lines indicate single-lens events that need the estimated t* to compute the detection efficiency. The length of the gray lines indicate 1σ uncertainties of estimated t*. The gold filled stars show the planetary events with measured t* in our sample, whereas the gold open stars are the planetary events with poorly determined t*. The blue stars are also planetary events with measured t*, but they are not included in our planet sample because they do not show enough anomaly signals in the MOA data. The red circles indicate the events with a clear finite-source effect. Most of them are high magnified single-lens events, but some events are stellar binary events that do not show the anomaly signals in the MOA data. The red filled and open circles are published and not-yet-published events, respectively.

Standard image High-resolution image

4.2. Detection Efficiency as a Function of Separation and Mass Ratio

We have computed the detection efficiency $\epsilon (\mathrm{log}s,\mathrm{log}q)$ for each of the 1474 events that passed the event selection criteria discussed in Section 3. Figures 35 show the detection efficiency for the planetary and ambiguous events. Typically, each figure has roughly triangular sensitivity contours as first shown by Rhie et al. (2000), and subsequently shown by a number of subsequent studies (Dong et al. 2006; Batista et al. 2009; Yee et al. 2009; Gould et al. 2010; Choi et al. 2012b). But, low-magnification events tend to show a double peak toward lower mass ratios at s > 1 and s < 1 and reduced detection efficiency at the Einstein ring radius (s = 1). These features were also seen in Bennett & Rhie (1996) and Gaudi et al. (2002).

Figure 3.

Figure 3.  Detection efficiencies $\epsilon (\mathrm{log}s,\mathrm{log}q)$ for the planetary events. The filled circles indicate the locations of the detected planets and the open circles are used for the events with the close/wide separation degeneracy.

Standard image High-resolution image
Figure 4.

Figure 4.  Same as Figure 3.

Standard image High-resolution image
Figure 5.

Figure 5.  Same as Figure 3, but the last panel (with the gray boxes) is the detection efficiencies for the ambiguous event, where the planetary model is used for the calculation of detection efficiency.

Standard image High-resolution image

When we sum the detection efficiency for all the events in our analysis, we obtain the survey sensitivity,

Equation (3)

which is plotted in Figure 6. The survey sensitivity $S(\mathrm{log}s,\mathrm{log}q)$ is the number of planet detections that we expect with our 1474 events sample,19 if we assume the uniform planet distribution in $\mathrm{log}s$ and $\mathrm{log}q$. If our survey was dominated by well-sampled, high-magnification events, then the survey sensitivity contours would be nearly symmetric in $\mathrm{log}s$ (Gould et al. 2010; Choi et al. 2012b). However, most of the sensitivity in our sample comes from more modest magnification events. This leads to higher sensitivity for planets with s > 1 than for planets with s < 1. This is due to the fact that outer planets perturb the higher magnification "major" image instead of the lower magnification "minor" image. For a single-lens event with magnification A, the major image has a magnification of $(A+1)/2$ while the minor image has a a magnification of $(A-1)/2$, so when the magnification approaches A = 1, we can only expect to detect planets that perturb the major image. As shown in Figure 6, most of the planets discovered in our data have separations close to the Einstein radius, s ∼ 1, and the distribution is pretty flat in $\mathrm{log}q$ with slightly more planets with 10−4 < q < 10−3 than between ${10}^{-3}\lt q\lt {10}^{-2}$. Since we are less sensitive to planets with small q values, this implies that lower mass planets are more common than more massive planets, at least down to q ∼ 10−4.

Figure 6.

Figure 6.  Sum of the detection efficiencies, or survey sensitivity $S(\mathrm{log}s,\mathrm{log}q)$, plotted as a function of $\mathrm{log}s$ and $\mathrm{log}q$. The contours indicate the number of planets detected if each lens star has a planet at the specified s and q values. The average detection efficiency is obtained by dividing the survey sensitivity by the total number of events, 1474. The 23 planets in the sample are plotted with the filled and open circles for the planets without and with the separation degeneracy, respectively. The close and wide solutions for each separation degeneracy event are connected with a line, and the ambiguous event is highlighted with light red halos.

Standard image High-resolution image

We can also define the survey sensitivity as a function of $\mathrm{log}q$, by integrating $S(\mathrm{log}s,\mathrm{log}q)$ over $\mathrm{log}s$,

Equation (4)

where ${s}_{-}\leqslant s\leqslant {s}_{+}$ is 0.1 ≤ s ≤ 10. Figure 7 shows the survey sensitivity, $S(\mathrm{log}q)$, where we have linearly interpolated the sensitivity between each calculated grid point. The statistical uncertainties in the sensitivity are negligible. The survey sensitivities estimated using the ρmin and ρmax, which are 1σ lower and upper values of ρ, are also plotted in Figure 7. We use the two-dimensional $S(\mathrm{log}s,\mathrm{log}q)$ to determine the exoplanet mass-ratio function in the following section, instead of the integrated, one-dimensional survey sensitivity, $S(\mathrm{log}q)$.

Figure 7.

Figure 7. Survey sensitivity as a function of $\mathrm{log}q$, derived by integrating $S(\mathrm{log}s,\mathrm{log}q)$ over $\mathrm{log}s$ over the range 0.1 ≤ s ≤ 10, plotted along with the left Y-axis scale. The red points show the survey sensitivity, and the black solid lines are the linearly interpolated sensitivity between the calculated grid points. The statistical uncertainties are plotted with error bars and dotted lines. Also, the sensitivity estimated with the 1σ lower and upper values of the finite-source size, ρmin and ρmax, are plotted in the orange and purple lines, respectively. The uncertainty of the finite-source size is negligible above q ∼ 10−3, so they are not plotted. If the survey sensitivity is divided by the total number of events 1474, then we obtain averaged detection efficiency, which is shown on the right Y-axis.

Standard image High-resolution image

5. PLANET FREQUENCY

The planet–star mass ratio, q, is measured in virtually all planetary microlensing events, while the masses can be sometimes determined by subtle microlensing parallax effects or follow-up observations with high angular resolution, but the masses of the lens stars and planets have not been measured for a large fraction of events. So, it is natural to determine the exoplanet frequency as a function of q. Similarly, the projected separation, s, is readily determined in units of the angular Einstein radius, and it can generally only be determined in physical units for the same events that allow the determination of lens star and planet masses. Our choice of parameters differs from that of Cassan et al. (2012), who attempted to determine the exoplanet frequency as a function of planet mass without the use of any mass measurements. So, the masses could only be estimated with a Bayesian analysis that assumed that the planet frequency does not depend on the host-star mass. We feel that simply reporting the mass-ratio function is a better representation of the constraints of the exoplanet population from microlensing observations.

5.1. As a Function of Separation and Mass Ratio

Figure 8 shows histograms of the the distribution of planet mass ratios in this analysis. The black histogram is the number of detected planets as a function of mass ratio, and the red histogram includes a detection efficiency correction that considers a Poisson probability. So, the red histogram represents our estimate of the true planetary distribution, with 1σ and 2σ upper limits in the non-detection mass-ratio bins. Note that these histograms are only for visualization. The unbinned data is used for the calculations that determine the mass-ratio function. There are several features that are apparent from this distribution. First, the detected distribution is pretty flat down to q ∼ 10−4 or 5 × 10−5. This implies that the frequency of planets increases at small q with a slope that nearly cancels the slope of the survey sensitivity function shown in Figure 7. However, it is also clear that this trend does not continue much below q ∼ 10−4. An extrapolation of the slope of the sensitivity-corrected mass-ratio function at q > 10−4 toward q < 10−4 passes well above the 2σ upper limits for the bins with q < 3 × 10−5. This indicates that the power-law mass-ratio functions used by Gould et al. (2010) and Cassan et al. (2012) will not be able to explain our data.

Figure 8.

Figure 8. Mass-ratio distribution of the planets in this analysis. The black histogram shows the number of the detected planets including the one ambiguous event. The number of detected planets corrected by the survey sensitivity in Figure 7 is plotted in red, with 1σ error bars. No planets were found in the bins with $\mathrm{log}q\lt -4.5$, so the upper limits with 1σ and 2σ confidence level are plotted with red and pink arrows, respectively. The best-fit broken power-law mass-ratio function, evaluated at s = 1, is shown in black, and the gray shaded region covers the area mapped out by the broken power-law mass-ratio function models that are consistent with the data at the 68% confidence level. The blue point and error bar indicate the median and 1σ range of the exoplanet mass function at q = 3 × 10−4.

Standard image High-resolution image

Another shortcoming of the previous microlensing studies of Gould et al. (2010) and Cassan et al. (2012) is that they did not have enough events to measure the separation dependence of the mass-ratio function. This difficulty is exacerbated by the close/wide degeneracy for the high-magnification events that dominate the Gould et al. (2010) and Cassan et al. (2012) samples. In particular, only two of the six planets in the Gould et al. (2010) sample and three of the eight planets in the Cassan et al. (2012) had the close/wide degeneracy resolved. In contrast, the close-wide degeneracy is resolved in 17 of the 23 events in our sample due to the fact that most of the planetary signals in our sample were discovered through the planetary-caustic channel as expected for a high-cadence survey. Therefore, our sample can be used to measure the exoplanet separation distribution.

The considerations above lead us to consider an exoplanet mass-ratio function consisting of a power law in projected separation, s, and a broken power law in the mass ratio, q,

Equation (5)

where Npl is an average number of planets around a lens star, n and p are the power-law indices for the mass ratio above and below the break in the mass-ratio function at q = qbr. The power-law index for the separation is denoted by m, and Θ is the (Heavyside) step function.

For events with multiple degenerate models that fit the data, we weight the different models with ${e}^{-{\rm{\Delta }}{\chi }^{2}/2}$, where Δχ2 is the χ2 difference between the model in question and the best-fit model. Mostly, these degenerate models are close and wide models, but we also use a similar procedure for our ambiguous planetary event OGLE-2011-BLG-0950/MOA-2011-BLG-336, as discussed in Section 3.3.

We use a likelihood analysis to determine the exoplanet mass-ratio function parameters, A, qbr, n, p, and m. It is straight-forward to derive the likelihood function from Poisson statistics (Alcock et al. 1996; Gould et al. 2010), and it is given by

Equation (6)

where Nexp is the number of planets detected for a mass-ratio function model, $f(q,s;A,{q}_{\mathrm{br}},n,p,m)$, with the measured survey sensitivity, $S(\mathrm{log}{s}_{i},\mathrm{log}{q}_{i})$, defined in Equation (3). The parameters qi and si refer to the mass ratio and separation (in Einstein radius units) for each of the Nobs observed planets. The uncertainties in qi and si are included as a two-dimensional Gaussian distribution. (This modified the model parameters by <0.1σ compared to a calculation that ignored the error bars.) The best-fit exoplanet mass-ratio function parameters are the ones that maximize the likelihood function, and they are shown in Table 4.

Table 4.  Best-fit Exoplanet Mass-ratio Function Parameters

Parameter qbr Free qbr Fixed
  MOA-only All MOA-only All
A 0.50 0.62 0.56 0.61
qbr × 104 1.97 1.65 ≡1.7 ≡1.7
n −0.96 −0.92 −0.94 −0.92
p 0.57 0.47 0.71 0.44
m 0.65 0.50 0.65 0.50

Download table as:  ASCIITypeset image

In order to derive confidence intervals for the exoplanet mass-ratio function parameters using a Bayesian approach, we must chose "prior distributions" of these parameters, which are equivalent to the integration measure that we use to integrate over the the likelihood function, ${ \mathcal L }(A,{q}_{\mathrm{br}},n,p,m)$. We chose uniform priors for n and m, while for A and qbr we chose priors uniform in $\mathrm{log}A$ and $\mathrm{log}{q}_{\mathrm{br}}$. The power index, p, for planets with small mass ratios, q < qbr, presents a more difficult choice of a prior. We only have a small number of planet detections for $q\lesssim {10}^{-4}$, so mass-ratio functions with qbr close to the q value of the lowest mass-ratio planet in our sample are not strongly excluded by our data set. Such models are consistent with very large p values, which would imply that planets with q ≲ 10−5 are extremely rare. But, we know that such planets exist from Kepler and solar system examples, and we also expect that it is relatively common for planets to change orbits due to gravitational interactions with other bodies. Similarly, values of p ≪ −1 are also unlikely due to planetary scattering considerations. Planet distributions determined from the Kepler survey indicate that small planets are quite common in shorter period orbits (Howard et al. 2012; Dong & Zhu 2013; Foreman-Mackey et al. 2014), so there should be some small planets in wider orbits, if only due to planet–planet scattering. Thus, it is sensible to chose a prior that disfavors priors with $| p| \gg 1$. Hence, we chose a prior distribution that is uniform in hyperbolic arcsine. This prior reduces $| p| $ by about a factor of two, and changes the other parameters by much less than their error bars. The uncertainties in q and s for each planetary event are included as the noise model for this analysis, as discussed above.

The maximization of the likelihood function was conducted with the MCMC algorithm to determine the best-fit parameters and uncertainties for each parameter using the Affine Invariant MCMC Ensemble sampler of Foreman-Mackey et al. (2013). The best-fit values are shown as crosses in Figure 9, and they are also listed in Table 4. Figure 9 also shows scatter plots and 68% and 95% confidence levels of the MCMC results, and Table 5 shows the median value and 68% confidence level for each parameter. Although the uncertainty of the power-law index, p is large, the power-law break at q ∼ 10−4 is clearly detected. The best-fit model with a broken power-law function is favored over the single power-law model by Δχ2 = 11.6 with two additional model parameters. The χ2 distribution with two degrees of freedom indicates that the single power-law model can be excluded with p-value of 0.0030. However, this use of the p-value to estimate the evidence for the additional model parameters has been criticized as being insufficiently conservative. Trotta (2008) argue that it is better to use the Bayesian evidence or Bayes factor, which includes some accounting of the Occam's razor preference for a simpler model. Following Trotta (2008), we can estimate the Bayes factor as ${ \mathcal B }=-1/({ \mathcal P }e\mathrm{ln}{ \mathcal P })$, where ${ \mathcal P }$ is p-value, and this yields a Bayes factor of ${ \mathcal B }=21$, implying that the broken power-law model is favored by a factor of 21 over the single power-law model.

Figure 9.

Figure 9. Likelihood contours for the broken double power-law function, $f=A\left[{\left(\tfrac{q}{{q}_{\mathrm{br}}}\right)}^{n}{\rm{\Theta }}(q-{q}_{\mathrm{br}})+{\left(\tfrac{q}{{q}_{\mathrm{br}}}\right)}^{p}{\rm{\Theta }}({q}_{\mathrm{br}}-q)\right]{\left(\tfrac{s}{{s}_{0}}\right)}^{m}$ for the MOA sample. The contours indicate the 68% and 95% confidence intervals, and the crosses indicate the best-fit parameters (which are also listed in Table 4).

Standard image High-resolution image

Table 5.  Exoplanet Mass-ratio Function Parameters from MCMC

Parameter qbr Free qbr Fixed
  MOA-only All MOA-only All
A ${0.65}_{-0.26}^{+0.41}$ ${0.95}_{-0.39}^{+0.48}$ ${0.56}_{-0.17}^{+0.22}$ ${0.61}_{-0.16}^{+0.21}$
B ${0.28}_{-0.08}^{+0.11}$ ${0.29}_{-0.07}^{+0.09}$ ${0.31}_{-0.08}^{+0.09}$ ${0.34}_{-0.08}^{+0.09}$
qbr × 104 ${1.19}_{-0.58}^{+1.14}$ ${0.67}_{-0.18}^{+0.90}$ ≡1.7 ≡1.7
n $-{0.91}_{-0.17}^{+0.15}$ $-{0.85}_{-0.13}^{+0.12}$ $-{0.96}_{-0.15}^{+0.14}$ −0.93 ± 0.13
p ${1.4}_{-1.1}^{+3.4}$ ${2.6}_{-2.1}^{+4.2}$ ${1.0}_{-0.6}^{+0.9}$ ${0.6}_{-0.4}^{+0.5}$
m ${0.62}_{-0.58}^{+0.55}$ 0.46 ± 0.49 ${0.62}_{-0.58}^{+0.55}$ ${0.49}_{-0.49}^{+0.47}$

Note. Note that B = f(q = 3 × 10−4, s = 1) is not an independent parameter, but it has a smaller dispersion than A.

Download table as:  ASCIITypeset image

Figure 8 shows the behavior of the exoplanet mass-ratio function at s = 1. The uncertainties of the estimated finite-source sizes for events without finite source are also statistically negligible. If we use the survey sensitivity estimated with ρmin values, the changes of the best-fit model parameters are less than 0.1σ. If we use the ρmax values, which reduce the sensitivity for low-mass-ratio planets, the normalization, A, changes by 0.1σ.

It is clear from Figures 8 and 9 that the amplitude and slope of the mass-ratio function above the break at q = qbr are pretty well constrained. But, the value of qbr and the amplitude A of the mass ratio (measured at qbr) are not determined so precisely. The power-law index, p, of the mass-ratio function below the break is even more poorly determined. This is due to the small number of planets found with q < 10−4 and the lack of planets with q < 3 × 10−5. The choice to normalize the mass-ratio function at q = qbr allows us to write the mass-ratio function, Equation (5), in a more compact way, but it also increases the apparent uncertainty in the mass-ratio function at larger q values. Therefore, we introduce the parameter $B=f(q=3\times {10}^{-4},s=1)$, shown in Table 5 and Figure 8 with the blue point, to indicate the uncertainty in the normalization at the mass ratios where microlens planets are more typically found.

5.2. As a Function of Event Timescale

Microlensing has the potential to determine the exoplanet mass function beyond the snow line as a function of host mass and galactocentric distance, but this requires additional information such as microlensing parallax (Bennett et al. 2008, 2010; Gaudi et al. 2008; Muraki et al. 2011; Skowron et al. 2015) or detection of the lens star (Bennett et al. 2006, 2007, 2015; Dong et al. 2009b; Batista et al. 2014, 2015). While such measurements exist for some of the planets in our sample, our analysis does not use this information. Nevertheless, we can still explore the dependence of the exoplanet mass-ratio function on the one microlensing parameter, tE, that depends on the mass and galactocentric distance of the host star.

Shorter tE events can be caused by low-mass lens stars or a combination of a high lens-source relative velocity and a small lens-source distance. These later conditions hold for lens systems in the Galactic bulge, and a systematic trend of planet properties with tE would likely indicate a trend in planet properties with host-star mass or galactocentric distance, or some combination of the two effects. (The mass and distance estimates alluded to in the previous paragraph are needed to break this degeneracy between lens mass and distance.)

In this section, we investigate the planet frequency as a function of tE. The median of tE is ∼20 days for microlensing events toward the Galactic bulge, but the median of tE for planetary events is ∼40 days, as shown in Figure 1. This effect is largely explained by the higher planet detection efficiency for long tE events, as shown in the left panel of Figure 10. But, it is possible that planets favor long tE events, even after correction for the detection efficiency. If planets favor long tE events, it could imply that cold planets are more common orbiting massive stars, Galactic disk stars, or both. Note that it is important that we have selected only events with robust tE measurements in order to avoid systematic errors associated with large tE uncertainties.

Figure 10.

Figure 10. Left: detection efficiencies as a function of mass ratio and tE. For each 0.5 dex × 0.5 dex bin, detection efficiencies are averaged. The filled circles show the 22 planetary events in this analysis, and the open circle indicates the ambiguous event. The contours show the averaged detection efficiencies of 40%, 20%, 10%, 5%, 2%, 1%, and 0.5%, as labeled from top to bottom. Right: the survey sensitivity, ${N}_{\mathrm{sens}}(\mathrm{log}{q}_{\mathrm{bin}},{t}_{{\rm{E}}})={\sum }_{i\,=\,1}^{{N}_{\mathrm{bin}}}\epsilon (\mathrm{log}{s}_{i},\mathrm{log}{q}_{i})$, which is given by the summation of the detection efficiencies for each observed event in 0.5 dex × 0.5 dex bins. While the detection efficiency is highest for very long events, there are few such events, so the survey sensitivity is highest for events in the tE ∼ 20–50 days range. The contours show the averaged number of expected planets of 100, 50, 20, 10, 5, 2, and 0.5 from top to bottom, as the labels indicate.

Standard image High-resolution image

The left panel of Figure 10 shows the detection efficiencies $\mathrm{log}\varepsilon (\mathrm{log}\,{t}_{{\rm{E}}},\mathrm{log}\,q)$ in the tEq plane. The detection efficiency is averaged over 0.5 dex bins in both $\mathrm{log}\,{t}_{{\rm{E}}}$ and $\mathrm{log}q$ for visualization. The detection efficiency is larger for long tE and higher q values. But, the number of events drops for tE ≳ 60 days due to the small number of long tE events in our sample. So, our sensitivity to planets in longer events drops, as shown in the right panel of Figure 10. This panel shows the survey sensitivity, which is the sum of the detection efficiency for each event, in 0.5 dex bins in both $\mathrm{log}\,{t}_{{\rm{E}}}$ and $\mathrm{log}q$. The survey sensitivity is highest for events with tE ∼ 20–50 days where both the detection efficiency and the number of events are relatively high.

We expect any trend with tE to be relatively weak, because tE depends on the lens primary mass and distance, as well as the lens-source relative transverse velocity. Thus, even a strong trend with host mass of galactocentric distance could be partially washed out by the variation of the other variables that tE depends on. With only 22 or 23 events, we do not have a great deal of sensitivity to subtle variations in the planet frequency with tE, so we divide the planetary events into two mass-ratio ranges: 14 planets with $-4.5\lt \mathrm{log}\,q\lt -3$ and 9 planets with $-3\lt \mathrm{log}\,q\lt -1.5$. The planets in the lower mass range have masses less than that of Jupiter mass, and are mostly cold Neptunes, while the planets in the later mass range are Jupiters and super-Jupiters. As in Section 5.1, we define ${f}_{{t}_{{\rm{E}}}}$ as a number of planets per star per logarithmic event timescale,

Equation (7)

where Npl is an averaged number of planets around a lens star for the two mass-ratio ranges. The parameters C0 and γ are the normalization and slope for the assumed power-law function, respectively. We select a pivot point of ${t}_{{\rm{E}}0}=30\,$ days. We use a likelihood analysis to determine C0 and γ with an assumption of uniform priors. The likelihood function is just the Poisson probability of finding the observed number of events, Nobs, times the product of the probability of finding events with each of the observed event timescale, tE i,

Equation (8)

where $\varepsilon ({t}_{{\rm{E}}i})$ is the detection efficiency in each 0.5 dex tE bin, and Nexp is the number of events expected for the given C0 and γ values, integrated over $0.5\lt \mathrm{log}{t}_{{\rm{E}}}\lt 2.5$. We find that the power-law index γ is flat for both mass-ratio ranges: γ = 0.14 ± 0.39 for q > 10−3 and γ = 0.07 ± 0.30 for the range with q < 10−3, as shown in Figure 11. So, the mass-ratio function is flat or slightly increasing with tE. The normalizations are ${C}_{0}={0.018}_{-0.006}^{+0.007}$ for q > 10−3 and ${C}_{0}={0.28}_{-0.07}^{+0.09}$ for q < 10−3. A larger survey will be evidently be needed to measure the tE dependence of the cold-exoplanet mass-ratio function.

Figure 11.

Figure 11. Planet occurrence frequency as a function of tE for high mass-ratio planets ($-3\lt \mathrm{log}q\lt -1.5$) on the left and low-mass-ratio planets ($-4.5\lt \mathrm{log}q\lt -3$) on the right.

Standard image High-resolution image

6. DISCUSSION

6.1. Comparison with Previous Microlensing Results

Previous attempts to measure the exoplanet mass ratio or mass function using microlensing data have used more limited samples than we have used in this paper. The first attempt was the measurement by Sumi et al. (2010) of the slope of the exoplanet mass-ratio function using only relative and not absolute detection efficiencies. Sumi et al. (2010) used the 10 microlensing discoveries known at that time to estimate logarithmic mass-ratio function slope of n = −0.7 ± 0.2, which is consistent with our value of $n=-{0.90}_{-0.17}^{+0.15}$ for the mass ratios >qbr. But, without the break in the mass-ratio function, it appears that the events with q < qbr have flattened their mass function slope.

Gould et al. (2010) used 13 high-magnification events containing 6 planetary signals for the first microlensing analysis with absolute detection efficiencies. These events were selected based on their observation by the μFUN microlensing follow-up collaboration They did not solve for a slope with mass ratio, and they found ${d}^{2}{N}_{\mathrm{pl}}/(d\mathrm{log}q\ d\mathrm{log}s)=0.36\pm 0.15$ at an average mass ratio of q = 5 × 10−4, assuming a distribution that is flat in both $\mathrm{log}q$ and $\mathrm{log}s$. As shown in Figure 12, the μFUN result is a factor of ∼2.0 higher than our result when evaluated at the central q for the μFUN analysis. Our value is $1.2{\sigma }_{\mu \mathrm{FUN}}$ below the μFUN result, and less than 1σ when the combined MOA and μFUN error bars are used. However, our mass-ratio function has a different functional form than theirs. So, a better comparison would be to integrate our mass-ratio functions over a range of mass ratios centered at q = 5 × 10−4. In order to avoid the complication associated with the break, we select the range extending from q1 = 2 × 10−4 to q2 = 1.25 × 10−3, a factor of 2.5 in each direction from the central value. Our mass-ratio function average over this q1q2 range is 0.226, which compares to the μFUN value of 0.36 ± 0.15. So, our value is 37% lower.

Figure 12.

Figure 12. Planet frequencies at the mass ratio of q = 5 × 10−4 as a function of semimajor axis, which is normalized by the snow line ∼2.7(M/M) au. Microlensing results are compared with the radial velocity result of Cumming et al. (2008). We assume that the exoplanet frequency scales with the mass ratio, q, and the separation compared to the snow line. The Gould et al. (2010) and our pivot point (s = 1) is plotted at three times the snow line.

Standard image High-resolution image

Our mass-ratio function is fully consistent with the earlier μFUN measurement, so there is no need to identify a bias to explain the fact that our value is lower. However, we have identified such a bias, nonetheless. This is a "publication date" bias. That is, the publication date for the μFUN sample was not decided in advance. Instead, it was selected, in part, because the authors judged that they had discovered enough planets for such a paper. Such a decision is more likely after a few years or two with a higher-than-average planet detection rate. However, we have identified a possible bias, nonetheless. Gould et al. (2010) found six planets in four years of data for a rate of 1.5 per year, but the following six years of observations yielded only three planets with μFUN data (Bachelet et al. 2012b; Yee et al. 2012; Fukui et al. 2015) that pass the criteria of Gould et al. (2010), when we would have expected nine. Event OGLE-2013-BLG-0341 does not qualify for this sample because of its stellar binary signal is much stronger than the signal of the planet (Gould et al. 2014). There are also several recent strong planetary events with important μFUN data that are excluded from this list of planets because they have $| {u}_{0}| \gt 0.005$ (Miyake et al. 2011; Han et al. 2013; Yee et al. 2014). So, it is likely that an update of the Gould et al. (2010) analysis will yield a result closer to our MOA result. The possibility of a similar bias affecting our own analysis seems to be smaller because our sample is larger.

While the previous analyses of Sumi et al. (2010) and Gould et al. (2010) dealt with exoplanet mass ratio, q, as our analysis does, Cassan et al. (2012) attempted to derive the exoplanet mass function. They found ${d}^{2}{N}_{\mathrm{pl}}/(d\,\mathrm{log}a\ d\,\mathrm{log}M)\,={{0.24}_{-0.10}^{+0.16}(M/{M}_{\mathrm{Sat}})}^{-0.73\pm 0.17}$, where a is semimajor axis, M is the planetary mass, and MSat is the mass of Saturn, MSat = 95.2 M. The PLANET Collaboration (Cassan et al. 2012) analysis includes constraints from Sumi et al. (2010) and Gould et al. (2010), and it uses a Bayesian analysis, rather than mass and distance measurements to convert from mass ratios to masses. This can be somewhat misleading, as it makes the assumption that exoplanet properties do not depend on the mass or distance of the host star. Since no mass measurements are used, we feel that it is more sensible to convert the Cassan et al. (2012) result back to mass ratios instead of masses. This can be done trivially because a power-law mass function transforms into a power-law mass-ratio function. (This applies for every mass, so it applies for all masses.) However, we must note that Cassan et al. (2012) do not correctly estimate the median mass that their survey probes for planets. They estimate the median mass corresponding the event with the median tE in the detection efficiency sample. However, the planet detection efficiency is higher for longer ${t}_{{\rm{E}}}$ events, so the sensitivity is increased for longer events. Since the planet occurrence rate is only weakly sensitive on tE, as shown in Figure 11, we can use the median value of tE for the detected planets, which is tE = 42.5 days for the 23 events in the MOA sample, as well as the combined MOA+μFUN+PLANET sample discussed below. With the Galactic model of (Bennett & Rhie 2002), this corresponds to a median mass of 0.59 M. This is more appropriate median mass for the PLANET sample, so we use it to convert the Cassan et al. (2012) mass function back into a mass-ratio function. This converted result (which is used for only this comparison, not for the following subsection) is compared to our result and the μFUN analysis in Figure 13. The PLANET result is obviously consistent with our mass-ratio function.

Figure 13.

Figure 13.  MOA exoplanet frequency result from Figure 8 as a function of mass ratio is compared with previous microlensing and RV results. The reference, type of primary star, and semimajor axis are labeled in the same color as the points, error bars, and shaded regions in the main body of the figure.

Standard image High-resolution image

The most recent statistical analysis of microlensing exoplanet data is the Shvartzvald et al. (2016) analysis of the combined MOA, OGLE, and Wise 2011–2014 microlensing survey data set. This analysis is rigorous in the sense that the detection efficiency for binary and planetary microlensing events is calculated, but the analysis of the detected binary and planetary microlensing events is incomplete. Some of the binary and planetary events have had a complete analysis done in other papers, and the results of these published analyses are used when they are available. But, for the other binary and planetary events, the light-curve analysis has been limited to a coarse grid in q and s. This results in light-curve models that obviously, by eye, do not fit the data for many of their binary and planetary events, including OGLE-2011-BLG-0481, OGLE-2011-BLG-0974, OGLE-2012-BLG-0442, OGLE-2014-BLG-0257, OGLE-2014-BLG-0572, OGLE-2014-BLG-0676, OGLE-2014-BLG-1255, and OGLE-2014-BLG-1720, which could lead to incorrect statistical inferences about planetary system properties. For example, OGLE-2011-BLG-0481 is a stellar binary, but the model presented in Shvartzvald et al. (2016) is a planetary model with q ∼ 0.014. OGLE-2014-BLG-0676 is correctly classified as a planetary event, but Shvartzvald et al. (2016) estimate a mass ratio of q ∼ 1.4 × 10−3 whereas the full analysis by Rattenbury et al. (2016) finds q = 4.78 ± 0.13 × 10−3, a difference of a factor of 3.4. They used nine planets including the above one to find a mass-ratio slope of ${dN}/d\,\mathrm{log}q={q}^{-0.50\pm 0.17}$, which is only marginally consistent with our slope of ${dN}/d\,\mathrm{log}q={q}^{n}$ with $n=-{0.91}_{-0.17}^{+0.15}$. They also found that the average number of planets per star is ${0.55}_{-0.22}^{+0.34}$ with 10−4.9 < q < 10−1.4 and 0.5 < s < 2. This compares to 0.43 planets per star with the same range of q and s values for our combined mass-ratio function, discussed in the next section.

6.2. Combined Microlensing Mass-ratio Function Analysis

Our analysis includes 22 planetary events and one likely planetary event, which compares to six and eight events used in the previous statistical analyses of exoplanet microlensing data sets (Gould et al. 2010; Cassan et al. 2012). However, seven of the eight events used in these previous analyses are not included in our sample. (The planetary event OGLE-2007-BLG-349 is included in all three samples.) Thus, we can increase our sample to 29–30 planets by adding these samples to the MOA sample. This is straight-forward to do, because Gould et al. (2010) and Cassan et al. (2012) have made their survey sensitivity functions, $S(\mathrm{log}s,\mathrm{log}q)$ and detection efficiencies, $\epsilon (\mathrm{log}s,\mathrm{log}q)$, available to us. (To compute the survey sensitivity as a function of the semimajor axis and planet mass, Cassan et al. (2012) converted their detection efficiencies as a function of $\mathrm{log}s$ and $\mathrm{log}q$ to those of physical parameters. We use their detection efficiencies as a function of $\mathrm{log}s$ and $\mathrm{log}q$ (A. Cassan, private communication).) We then simply multiply the likelihood functions, Equation (6), for the three different surveys together to get our final likelihood function. The only complication is the overlap between the three surveys. We deal with this by removing the 2007 season from the PLANET sample and removing the events in the μFUN analysis from the MOA sample. We make these choices because the MOA-2007 season has a higher sensitivity than the PLANET 2007 season and because the inclusion of the follow-up data gives the μFUN analysis higher sensitivity for the 13 events in their sample, despite their higher detection threshold.

The results of this combined analysis are summarized in Tables 4 and 5, as well as Figures 14 and 15. The mass-ratio function parameters from this combined sample are very similar to the parameters from the MOA-only sample. Again, the uncertainties in mass ratio and separation for the additional planets from Gould et al. (2010) and Cassan et al. (2012) are negligible. They change the best-fit model by 10% of the uncertainties for each parameter. With this combined sample, the best-fit broken power-law model is a better fit than the best single power-law model by Δχ2 = 12.2 (with two additional parameters). This implies a p-value of 0.0022, and a Bayes factor of 27, implying that the broken power-law model is favored over the single power-law model by a factor of 27. The normalization is increased by 5% or ∼0.13σ when qbr is fixed, and the qbr value is decreased somewhat when qbr is set as a free parameter. As expected, the larger sample decreases the error bars slightly.

Figure 14.

Figure 14. Mass-ratio distribution, as in Figure 8, for the combined microlensing data set. The black histogram is the number of detected planets per bin, and the red histogram is corrected for detection efficiencies.

Standard image High-resolution image

We can also use this combined sample to look at the number of planets implied by our mass-ratio function. Microlensing has been traditionally thought to be sensitive to planets in a so-called "lensing zone" (Gould & Loeb 1992; Bennett & Rhie 1996), which roughly spans the range 0.5 < s < 2. If we integrate our best-fit mass-ratio function for the combined sample (shown in Figure 15 and Table 4), we find 0.48 planets per stars inside the "lensing zone." However, some of these planets have q < 5 × 10−5, which is below the lowest mass ratio for any planet in our combined sample. If we consider only planets in the range of our known microlens planets, with q > 5 × 10−5, we find 0.33 planets per star. On the other hand, 10% of the combined microlens exoplanet sample are located outside the nominal 0.5 < s < 2 "lensing zone." To include all the planets in our sample, we must use an extended "lensing zone" spanning the range 0.3 < s < 5. Our mass-ratio function predicts 0.79 cold planets per star with q > 5 × 10−5 and 1.26 cold planets per star if we extrapolate our mass function down to q = 0.

Figure 15.

Figure 15. Likelihood contours for the broken double power-law function $f=A\left[{\left(\tfrac{q}{{q}_{\mathrm{br}}}\right)}^{n}{\rm{\Theta }}(q-{q}_{\mathrm{br}})+{\left(\tfrac{q}{{q}_{\mathrm{br}}}\right)}^{p}{\rm{\Theta }}({q}_{\mathrm{br}}-q)\right]{\left(\tfrac{s}{{s}_{0}}\right)}^{m}$ for the full MOA+μFUN+PLANET sample. The contours indicate the 68% and 95% confidence intervals, and the crosses indicate the best-fit parameters (which are also listed in Table 4).

Standard image High-resolution image

6.3. Comparison with RV Results

Our microlensing mass-ratio function results are more readily compared with the results of RV surveys than transit surveys because both microlensing and the RV method detect exoplanet masses rather than the radius, which is detected by the transit method. We consider five different RV studies in refereed journals that report statistical results (Cumming et al. 2008; Mayor et al. 2009; Howard et al. 2010; Johnson et al. 2010; Bonfils et al. 2013; Montet et al. 2014). All of these studies have smaller exoplanet samples than our combined samples, except for Cumming et al. (2008) and possibly Mayor et al. (2009), and none of them have an exoplanet sample more than a factor of two larger than our samples.

Two of these studies (Mayor et al. 2009; Howard et al. 2010) focus on planets orbiting stars of ∼1 M in short period orbits (period <50 days), and so they have no overlap with our planetary microlens sample. However, their sensitivity does extend down to q ∼ a few × 10−5, where these results are consistent with the mass-ratio function measured with our microlensing at much larger orbital separations, as shown in Figures 13 and 17.

The study by Cumming et al. (2008) had the largest sample (48 planets) and covered periods of P < 2000 days (corresponding to a ≲ 3.1 au for solar mass host stars), but they have little sensitivity to planets less massive than Jupiter in wide orbits. So, they were sensitive primarily to planets with shorter period orbits than the planets probed by the microlensing method, and their host stars are more massive than most of the planets found by microlensing. Nevertheless, our results do agree well with their slope with separation, as indicated in Figures 12 and 16. While the Cumming et al. (2008) slope with semimajor axis was consistent with the microlensing result of Gould et al. (2010), the agreement with our result is somewhat better. However, our constraint on the slope of the mass-ratio function is rather weak, m ≃ 0.5 ± 0.5 for the complete sample, so it is consistent with a logarithmically flat distribution (Öpik's law) at 1σ. This is due to the relatively narrow range of planetary separations that most microlensing observations are sensitive to, as well as the close-wide degeneracy that occurs for high-magnification events that are sensitive to planets at a wide range of separations.

Figure 16.

Figure 16. Same as Figure 12, but the result of the full microlensing sample is compared to the Cumming et al. (2008) separation compared to the snow line.

Standard image High-resolution image

In contrast, microlensing is sensitive to a large range of exoplanet mass ratios. Our detected planets span a range of 500 in mass ratio, and we have significant sensitivities at masses below the planet with the lowest observed mass ratio. The Cumming et al. (2008) result does not match our more precise measurement of the dependence of the mass-ratio function on q at q ≲ 10−3, where their value is lower than our measurement, as our mass-ratio function rises more steeply toward lower mass ratios than theirs does, as indicated in Figures 13 and 17. This may be an indication that there are more cold Neptunes beyond the snow line that have undergone little or no migration than there are at ∼1 au.

Figure 17.

Figure 17.  Combined MOA+PLANET+μFUN exoplanet frequency result from Figure 14 as a function of mass ratio is compared with previous microlensing and RV results. The reference, type of primary star, and semimajor axis are labeled in the same color as the points, error bars and shaded regions in the main body of the figure.

Standard image High-resolution image

The RV studies of Johnson et al. (2010), Bonfils et al. (2013), and Montet et al. (2014) focused on M-dwarf primaries, although with rather small samples of 5, 11, and 14 planets, respectively. Johnson et al. (2010) found a significantly lower occurrence frequency than we find, but this is largely due to the limited mass and separation sensitivity of their RV survey. In contrast, using some of the same data, Montet et al. (2014) found an occurrence rate consistent with our measurements (see Figures 13 and 17). They have been much more careful to consider the RV selection effects and the possible signals of planets with periods longer than the survey. Bonfils et al. (2013) also searched for the planets orbiting M-dwarfs using the HARPS instrument, and found planets spanning a range of 300 in mass and 150 in semimajor axis, but only two of these planets are in the range of periods, 103–104 days, where microlensing is sensitive. However, only one of the two planets were confirmed by the HARPS survey (as pointed out in Clanton & Gaudi 2014b). So, we consider only one planet discovery to plot the Bonfils et al. (2013) results in Figures 13 and 17.

A recent study by Clanton & Gaudi (2014a, 2014b), simulated RV observations of the planet population found by microlensing showed that the number of RV detections expected by the simulation agrees with the actual RV observation of Bonfils et al. (2013). However, this was a comparison to the Cassan et al. (2012) mass function, so it is not a direct comparison with our results.

High-contrast imaging is a promising method to determine the frequency of exoplanets in wide orbits, but their sensitivity is currently limited to self-luminous planets, and their published statistical results tend to be upper limits (Rameau et al. 2013). So, this method does not yet have a great deal of overlap with the primary sensitivity region for the microlensing method at 2–3× the snow line.

6.4. Comparison with Kepler Results

By far, the largest statistical samples of planets with well-characterized detection efficiencies come from the Kepler mission (Petigura et al. 2013a; Foreman-Mackey et al. 2014; Dressing & Charbonneau 2015), but Kepler provides only a direct measurement of the planet radius and not mass for the vast majority of planets discovered. Not only that, but the Kepler planets for which both the radius and mass are measured are biased on samples of planets. Masses are more readily determined by RVs for more massive planets in short period orbits, but the smaller planets are considered to be more interesting. So, much observing time is devoted to planets that are thought to have barely detectable RV signals. Transiting planet masses can also be determined by transit-timing variations, but many of the transit-timing variation measurements are made for an unusual class of low-density planets that are tightly packed in highly co-planar orbits (Lissauer et al. 2013; Jontof-Hutter et al. 2014). These planets might have different compositions and therefore a different mass–radius distribution from planets in more sparsely packed, less co-planar systems. Also, attempts to determine the mass–radius relation have generally not spanned the full range of planet masses. At masses near a Jupiter mass and above, the physics of electron degeneracy pressure conspires to make the planet radius nearly independent of the planet mass (but still dependent on the planet composition).

As a result of these considerations, we have limited our comparison to Kepler to a simple comparison of the masses of the mass function breaks or peaks seen in each data set. For Kepler, we use the exoplanet radius functions of Petigura et al. (2013a) for G and K host stars and Dressing & Charbonneau (2015) for M-dwarf hosts. To convert to masses, we use the probabilistic mass–radius relation of Wolfgang et al. (2016) in both cases. The exoplanet mass functions derived from Petigura et al. (2013a) and Dressing & Charbonneau (2015) are shown in Figure 18. Only the red and blue colored part of the curves in this figure are reliable. Both of these exoplanet mass functions have a similar shape in the region where both are reliable, 2.5 M < Mpl < 13 M. The mass functions for M, and GK host stars have apparent peaks at ∼6 M and ∼8 M, respectively.

Figure 18.

Figure 18. Planet frequencies for M- and GK-dwarf hosts as a function of mass from the Kepler mission, based upon the Dressing & Charbonneau (2015) and Petigura et al. (2013a, 2013b) exoplanet radius functions for M and GK stars, respectively. These have been converted to masses using the probabilistic mass–radius relation of Wolfgang et al. (2016). The red and blue colored curves indicate the reliable region of these curves, and the gray colored regions indicate areas where the input radius functions or the mass–radius relation are thought (by the respective authors) to be unreliable.

Standard image High-resolution image

Because there are few microlens planets below the mass-ratio function break at qbr ∼ 10−4, we do not have a very precise measure of qbr. From Table 5, we see that the 1σ range is roughly (0.5–2) × 10−4. Since the typical host-star mass is ∼0.6 M, for our detected planets (in both the MOA and Combined samples), qbr = 10−4 corresponds to ∼20 M, which is just above the mass of Neptune. Thus, the break in the mass function beyond the snow line is likely to occur at a mass in the range 10–40 M. This is likely to be somewhat higher than the apparent mass function peaks in shorter period orbits, as estimated from Kepler data, but we will need a larger sample of low-mass microlens planets to confirm that the break in the mass function at higher masses beyond the snow line.

We can also compare the absolute numbers of planets found by Kepler at separations well inside the snow line to our values well outside the snow line. We focus on the planets near the peak of the mass-ratio function and the Kepler radius functions. From the Petigura et al. (2013a) sample, we consider planets in the period range 25 days < P < 200 days. This spans a factor of eight in period, which is equivalent to a factor of four in separation, by Kepler's third law. Using the mass–radius relation, we select lower limits for Rp > 2 R for the Petigura et al. (2013a) sample as equivalent to a lower limit of q > 2.6 × 10−5 in our analysis. We find 0.27 planets per star when we sum the Petigura et al. (2013a) bins between 25 days < P < 200 days with Rp > 2 R. The equivalent range for our exoplanet mass-ratio function (Equation (5) with the parameters from the qbr free-All column from Table 4) is 2.6 × 10−5 < q < 0.03 and 0.5 < s < 2. Integrating over this range yields 0.38 planets per star, suggesting that the peak of the exoplanet mass function of cold planets might have ∼41% more planets than an equivalent size region near at the peak of the warm planet distribution for G and K stars. We can also do a similar calculation with the tabulated radius function of Dressing & Charbonneau (2015) for early M-dwarf host stars. To match the bins of their tabulated radius distribution, we select 0.45 < s < 2.22 and 2.9 × 10−5 < q < 0.03 to span the same logarithmic range in separation and mass ratio as planets at the peak of the M-dwarf exoplanet radius function with 10 days < P < 110 days with Rp > 2 R. Integrating over this (q, s) gives 0.43 planets per star found by microlensing, as compared to a whopping 3.03 planets per star near the peak of the early M-dwarf radius function. This dramatic difference between the height of the peak in the radius function is also seen by Burke et al. (2015), but they also point out that some potentially large systematic errors may still exist in these analyses.

6.5. Masses of Host Stars

The masses of the lens stars are not known for the vast majority of microlensing events in our sample, but the situations is somewhat better for the lens stars found to host planets. For most planetary events, finite-source effects are seen, and this allows the angular Einstein radius to be determined. When the angular Einstein radius is known, the lens mass can be determined in those cases where the microlensing parallax effect can be measured. This has been done for six events in our sample, MOA-2007-BLG-192 (Bennett et al. 2008), OGLE-2007-BLG-349/MOA-2007-BLG-379 (Bennett et al. 2016), MOA-2009-BLG-266 (Muraki et al. 2011), MOA-2010-BLG-117 (D. P. Bennett et al. 2016, in preparation), MOA-2010-BLG-328 (Furusawa et al. 2013), and MOA-2011-BLG-197/OGLE-2011-BLG-0265 (Skowron et al. 2015), as well as one event with two planets, OGLE-2006-BLG-109 (Gaudi et al. 2008; Bennett et al. 2010), from the Gould et al. (2010) sample. Alternatively, it is also possible to determine lens mass with a measurement of the angular Einstein radius or microlensing parallax along with the measurement of the host-star brightness. With the help of an empirical or theoretical mass–luminosity relation, the host-star brightness can be turned into a lens mass measurement. This has been done using microlensing parallax for OGLE-2012-BLG-0950/MOA-2012-BLG-527 in our sample (Koshimoto et al. 2016). It also has been done using the angular Einstein radius for OGLE-2005-BLG-169 (Bennett et al. 2015; Batista et al. 2015) and OGLE-2006-BLG-109 (Gaudi et al. 2008; Bennett et al. 2010) in the Gould et al. (2010) sample and for OGLE-2005-BLG-071 (Dong et al. 2009b) in the Cassan et al. (2012) sample. The observations that detect the lens are often taken many years after the event, as this allows us to check that the putative lens star motion is actually consistent with being the lens. We plan to include these mass constraints in a future analysis.

The host-star masses in the planetary events in our sample mainly range from G-dwarfs to late M-dwarfs. The primary lens masses in MOA-2011-BLG-262 (Bennett et al. 2014) could be a Jupiter-mass object, although a low-mass star is more likely. Also, the primary lenses for OGLE-2007-BLG-224/MOA-2007-BLG-163 (Gould et al. 2009) and MOA-2011-BLG-274 (Freeman et al. 2015), both of which are single-lens events and used to estimate the detection efficiency, are very likely to be brown dwarfs or possibly a planetary mass object in the case of MOA-2011-BLG-274. Thus, conservatively, the host-star mass in our analysis is ranging over late-type stars extending to both G-dwarfs and brown dwarfs, or even planetary mass objects.

7. SUMMARY

We have reported the planet frequency as a function of planet-to-star mass ratio, separation, and event timescale, derived from the microlensing survey data by MOA-II telescope in 2007–2012. This sample included 1474 well-sampled microlensing light curves, and we find signals of 22 planets and one likely, but not certain, planet (see Section 3). Our sample probes the planetary systems of stars that can serve as lenses for gravitational microlensing surveys toward the Galactic bulge, so it is an average over stars that orbit inside the solar circle. It is possible that this population is somewhat different from the neighbors that can be probed by the other detection methods.

Our results are consistent with the previous statistical microlensing analyses of Sumi et al. (2010), Gould et al. (2010), and Cassan et al. (2012), but we find that, contrary to previous analyses, the data are not consistent with a single power-law mass function in the mass ratio, q. Instead, we find a strong change in the slope of the mass-ratio function at qbr ∼ 10−4, which we model as a broken power-law mass function, given by Equation (5) or

Equation (9)

with the best-fit parameters given in Table 4 and the Markov Chain medians and error bars given in Table 5. We present results for the 23-planet MOA sample, and the full sample of 30 planets (including one somewhat ambiguous event). The p values for the detection of the break in the mass-ratio function are 0.0030 and 0.0022 for the MOA and full samples, respectively. A more conservative estimate of the need for the power-law break in the mass function ratio is given by the Bayes factors of 21 and 27 favoring the broken power law for the MOA and full samples. These Bayes factors imply that the broken power-law model is favored over the single power-law model by these factors.

We present our results in terms of the planetary mass ratio, q, instead of mass because the mass ratio is measured robustly in virtually all planetary microlensing events. Converting mass ratios to masses requires the use of a Galactic model and an assumption about how the prevalence of planets might depend on Galactocentric distance and host-star mass. The common assumption is that there is no dependence on either, but this is motivated by our nearly complete ignorance of these dependencies. By presenting our results in terms of the mass ratio, we avoid invalidating them by making incorrect assumptions regarding the distance and host-mass dependence of the mass-ratio function. Nevertheless, it is sometimes useful to consider what mass might correspond to a particular mass ratio, such as the mass ratio of the break in the mass-ratio function. Since our typical planet host star has a mass of ∼0.6 M, qbr ∼ 10−4 corresponds to ∼20 M, with an uncertainty of about a factor of two. This break in the mass-ratio function seems to be at a somewhat larger mass than the breaks or peaks that appear in Kepler populations, with planetary radii converted to masses following Wolfgang et al. (2016), but a larger sample of low-mass microlens planets is needed to confirm this.

Our analysis is the most comprehensive systematic statistical analysis of the exoplanet occurrence frequency using planetary signals discovered in a microlensing survey. The previous analyses relied on planet detections by follow-up groups, which have a lower detection rate, and therefore, smaller samples. In the years since the MOA-II high-cadence survey started, several other high-cadence surveys have commenced, including the OGLE-IV survey (Udalski et al. 2015b), the Wise survey (Shvartzvald et al. 2016), and the KMTNet survey (Park et al. 2012).

An even more powerful improvement to the statistical study of planetary systems found by microlensing will be to include all the existing constraints on the mass and distance to the lens systems. In particular, the lens system masses can often be measured when the microlensing parallax effect is detected (Bennett et al. 2008, 2010; Gaudi et al. 2008; Muraki et al. 2011; Skowron et al. 2015), which can be done much more easily with follow-up telescopes located far from the Earth (Udalski et al. 2015c; Street et al. 2016). For events without giant source stars, it is also possible to determine the lens system mass and distance by detecting the lens star in high angular resolution follow-up observations (Bennett et al. 2006, 2007, 2015; Dong et al. 2009b; Batista et al. 2014, 2015). The inclusion of these constraints will allow us, for the first time, to determine the dependence of the properties of cold exoplanets as a function of host-star mass and Galactocentric distance.

The ultimate word on the statistical properties of planetary systems will be achieved from the space-based exoplanet survey (Bennett & Rhie 2002) of the WFIRST (Spergel et al. 2015) mission, and hopefully also the Euclid (Penny et al. 2013) mission. The high angular resolution of these space telescopes will allow mass and distance determinations of thousands of exoplanets because it will be possible to detect the lens star and measure the lens-source relative proper motion with the high resolution survey data itself. This will give us the same comprehensive picture of the properties of cold exoplanets that Kepler is providing for hot planets.

D.S. and D.P.B. acknowledge support from NASA grants NNX13AF64G, NNX14AG49G, and NNX15AJ76G. T.S. acknowledges the financial support from the JSPS, JSPS23103002, JSPS24253004, and JSPS26247023. D.S. was supported by Grant-in-Aid for JSPS Fellows. A.Y. acknowleges the financial support from the JSPS, JSPS25870893. The MOA project is supported by the grant JSPS25103508 and 23340064. This work was performed in part under contract with the California Institute of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. The authors thank OGLE, μFUN, MiNDSTEp, PLNAET, and RoboNet collaborations for letting us use their data to characterize the alerted events.

Footnotes

  • 16 

    They used 13 high-magnification events from four years of observations.

  • 17 

    They used 199 events from six years of observations.

  • 18 
  • 19 

    The survey sensitivity and detection efficiency can be obtained by contacting the first or second author.

Please wait… references are loading.
10.3847/1538-4357/833/2/145