Free-floating Planet Mass Function from MOA-II 9 yr Survey toward the Galactic Bulge

Gravitational microlensing observations toward the Galactic bulge enable exoplanet searches (Mao 1991; Gaudi et al. 2008; Bennett et al. 2010; Suzuki et al. 2016; Koshimoto et al. 2021b) and the measurement of the stellar and substellar mass functions (MFs; Paczyński 1991; Sumi et al. 2011; Mróz et al. 2017, 2019, 2020).

Sumi et al. (2011) first interpreted the detection of short Einstein radius crossing time (0.5 < t_E/day < 2) microlensing events as evidence for the existence of a population of free-floating planets (FFPs) and/or wide-orbit planets. While that analysis was limited by the small number of events found in a 2 yr subset of the survey by the Microlensing Observation in Astrophysics (MOA) group (Sumi et al. 2003) in collaboration with the Optical Gravitational Lensing Experiment (OGLE; Udalski et al. 1994), it opened up the field of FFP studies using microlensing.

Mróz et al. (2017) extended the work by using a larger sample from 5 yr of the OGLE survey. They discovered six events with timescales shorter (t_E ∼ 0.2 days) than those in the previous work. These events are separated from the longer events by a gap around t_E ∼ 0.5 days, which implies the possibility of a population of several Earth-mass FFPs.

These studies are based on the distribution of t_E, in which t_E is proportional to the square root of the lens mass M as follows:

$$\begin{eqnarray}\displaystyle \begin{array}{rcl}{t}_{{\rm{E}}} & = & \frac{\sqrt{\kappa M{\pi }_{\mathrm{rel}}}}{{\mu }_{\mathrm{rel}}}\\ & = & 0.1\,\mathrm{days}{\left(\frac{M}{5{M}_{\oplus }}\right)}^{1/2}{\left(\frac{{\pi }_{\mathrm{rel}}}{18\mu \mathrm{as}}\right)}^{1/2}{\left(\frac{{\mu }_{\mathrm{rel}}}{5\mathrm{mas}\,{\mathrm{yr}}^{-1}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 1 }$$

Here κ = 4G/(c²au) = 8.144 masM_⊙ ⁻¹ and we expect t_E ∼ 0.1 days assuming the typical value of the lens−source relative parallax, ${\pi }_{\mathrm{rel}}={\pi }_{l}^{-1}-{\pi }_{s}^{-1}=1$ au $({D}_{l}^{-1}-{D}_{s}^{-1})=18$ μas, for the bulge lens and a typical value of the lens−source relative proper motion in the direction of the Galactic center of μ_rel = 5 mas yr⁻¹. The lens mass M, the distance D _l to the lens, and the relative proper motion μ_rel are degenerate in the observable t_E. (D_s is the distance to the source star.) This means that the MF of the lens population has to be determined statistically, assuming a model of the star population density and velocities in the galaxy.

Mróz et al. (2018) found the first short-timescale (t_E = 0.32 days) event showing the finite-source (FS) effect, i.e., an FS and a single point lens (FSPL), in which one can measure an FS parameter ρ = θ_*/θ_E. Here θ_* is the angler source radius that can be estimated from an empirical relation with the source magnitude and color. θ_E is the angular Einstein radius, given by

$$\begin{eqnarray}&&{\theta }_{{\rm{E}}}=\displaystyle \frac{{\mu }_{\mathrm{rel}}}{{t}_{{\rm{E}}}}=\sqrt{\kappa M{\pi }_{\mathrm{rel}}}.\end{eqnarray} \tag{ 2 }$$

This value of θ_E can give us an inferred mass of the lens with better accuracy, as we can eliminate one of the threefold degenerate terms that affect t_E, namely, μ_rel:

$$\begin{eqnarray}&&M=\displaystyle \frac{{\theta }_{{\rm{E}}}^{2}}{\kappa {\pi }_{\mathrm{rel}}}=5{M}_{\oplus }{\left(\displaystyle \frac{{\theta }_{{\rm{E}}}}{1.5\mu \mathrm{as}}\right)}^{2}{\left(\displaystyle \frac{{\pi }_{\mathrm{rel}}}{18\mu \mathrm{as}}\right)}^{-1}.\end{eqnarray} \tag{ 3 }$$

While the inclusion of the angular Einstein radius, θ_E, enables tighter constraints on the lens masses, it adds a complication to a statistical analysis of FFP properties because the microlensing event detection efficiency depends on both t_E and θ_E (or equivalently t_E and ρ).

So far, six short FSPL events have been discovered (Mróz et al. 2018; Mróz et al. 2019, 2020, 2020; Kim et al. 2021; Ryu et al. 2021). All of these have θ_E < 10 μas, implying that their lenses are most likely of planetary mass. All of these sources are red giants with the exception of the subgiant source for OGLE-2016-BLG-1928 because their angular radii, i.e., cross section, are significantly larger than those of main-sequence (MS) stars.

Mróz et al. (2020) found the short FSPL event, OGLE-2016-BLG-1928, with the smallest value of θ_E = 0.842 ± 0.064 μas to date. Its lens is the first terrestrial-mass FFP candidate and the first evidence of such a population.

Kim et al. (2021) began a new approach to probing the FFP population by focusing on analyzing the θ_E distribution in events with giant sources. Ryu et al. (2021) found a gap at 10 < θ_E/μas < 30 in the cumulative θ_E distribution, which suggests a separation between the planetary-mass population and other known populations, like brown dwarfs (BDs).

Gould et al. (2022) completed the analysis of 29 FSPL giant source events found in the 2016–2019 KMTNet survey. They presented the θ_E distributions down to θ_E = 4.35 μas and confirmed that there is a clear gap in the distribution of θ_E at 9 < θ_E/μas < 26. They note that it is consistent with the gap in the t_E distribution shown by Mróz et al. (2017), indicating the existence of the low-mass FFP population. They used what they refer to as a "relative detection efficiency" that depends only on θ_E, but not t_E, to model the θ_E distribution with a power-law MF for the FFP and found ${{dN}}_{\mathrm{FFP}}/d\mathrm{log}M={(0.4\pm 0.2)(M/38{M}_{\oplus })}^{-p}$ dex⁻¹ star⁻¹, using a power law with 0.9 ≲ p ≲ 1.2. This range of the power, p, was estimated based on consideration of possible formation mechanisms, rather than a measurement. This would imply that the number of FFPs is at least an order of magnitude larger than the number of known bound planets.

We note that the Gould et al. (2022) result cannot be considered a measurement for a the following reasons. First, the true detection efficiency depends on both t_E and θ_E, and it is difficult to see how any selection criteria could remove the t_E dependence. As we discuss below in Section 4.1.1 and in Koshimoto et al. (2023), one can integrate over the t_E dependence of the detection efficiency to obtain an integrated detection efficiency. However, the integration over short t_E values depends on the FFP MF. However, Gould et al. (2022) seem to avoid this difficulty by simply adopting an analytic formula for the "relative detection efficiency" depending only on θ_E. The Gould et al. (2022) paper gives no justification for this analytic formula.

In this paper, we present the distributions of θ_E and t_E values for the microlensing events toward the Galactic bulge from 9 yr of the MOA-II survey. We also present the first measurement of MF of the planetary-mass objects using the t_E distribution. We describe the data in Section 2. We show the θ_E distribution in Section 3. We present the t_E distribution and the best-fit MF in Section 4. The discussion and conclusions are given in Section 5, and we compare the integrated detection efficiency in the Appendix.

We use the microlensing sample selected from the MOA-II high-cadence photometric survey toward the Galactic bulge in the 2006–2014 seasons (Koshimoto et al. 2023). MOA-II uses the 1.8 m MOA-II telescope, which has a 2.18 deg² field of view and is located at the Mt. John University Observatory, New Zealand. ¹⁴

Koshimoto et al. (2023) used an analysis method similar to what was used by Sumi et al. (2011, 2013), but including a correction of systematic errors and taking into account the FS effect. They applied a detrending code to all light curves to remove the systematic errors that correlate with seeing and air mass owing to differential refraction, differential extinction, and relative proper motion of stars in the same way as in Bennett et al. (2012) and Sumi et al. (2016). These corrections are important, as they result in higher confidence in the light-curve fitting parameters.

Koshimoto et al. (2023) selected light curves with a single instantaneous brightening episode and a flat constant baseline, which can be well fit with a point-source point-lens (PSPL) microlensing model (Paczyński 1986). In addition to PSPL, they modeled the events with an FSPL model (Bozza et al. 2018), which is especially important for short events. These are the major improvements compared to the previous analysis in Sumi et al. (2011, 2013), in addition to the extension of the survey duration.

Although they identified 6111 microlensing candidates, they selected only 3554 and 3535 objects as the statistical sample using the two relatively strict criteria CR1 and CR2, respectively. Here CR2 was defined as the stricter criteria compared to their nominal criteria CR1 to check the effect of the choice of the criteria on a statistical study. These strict criteria ensure that t_E is well constrained for each event and reject any contamination.

Sumi et al. (2011) reported 10 short events with t_E < 2 days in the 2006–2007 data set. Only five and four events survived following the application of CR1 and CR2, respectively. This is because the fitting results changed as a result of the re-reduction of the data set. On the other hand, two events are newly found, resulting in seven and six events following the application of CR1 and CR2, respectively. As a result, the excess at t_E = 0.5–2 days in the t_E distribution is not significant anymore; however, an even shorter event, MOA-9y-6057 (t_E = 0.22 ± 0.06 days), is added.

There are 13 FSPL events with θ_E measurements in the sample, including two FFP candidates, MOA-9y-5919 and MOA-9y-770, that have terrestrial and Neptune masses, respectively. See Koshimoto et al. (2023) for the light curves and detailed parameters of the 13 events.

The red line in Figure 1 indicates the cumulative distribution of θ_E from Table 7 of Koshimoto et al. (2023). The black line indicates the distribution of 29 FFPs by Gould et al. (2022) normalized to 13 events as a comparison. Although these cannot be directly compared because these are not corrected for detection efficiencies, the general trends seen in Figure 1 may give us some insights.

Zoom In Zoom Out Reset image size

The distributions are consistent for θ_E > 30 μas, where the effects of the detection efficiencies are likely small. There is a gap around 5 < θ_E/μas < 70, which is roughly consistent with the gap at 10 < θ_E/μas < 30 found by Ryu et al. (2021) and Gould et al. (2022). This gap confirmed the existence of the planetary-mass population as distinct and separated from the stellar/BD population as indicated by Gould et al. (2022).

The MOA cumulative distribution shows fewer events over 30 < θ_E/μas < 70 compared to Gould et al. (2022). This may be just due to the small number of statistics. But note that Koshimoto et al. (2023) found a BD candidate MOA-9y-1944 with θ_E = 46.1 ± 10.5 μas, although this is not in the final sample for statistical analysis because the source magnitude of I _s = 21.91 mag is fainter than the threshold of I _s < 21.4 mag.

In our sample, there is one event with a very small value of θ_E of 0.90 ± 0.14 μas. This confirms the existence of the terrestrial-mass population, which gives rise to events such as OGLE-2016-BLG-1928, which has θ_E = 0.842 ± 0.064 (Mróz et al. 2020). These values are significantly smaller than the lower edge of θ_E ∼ 4.35 μas as reported in Gould et al. (2022). This is partly a result of selection bias given that Gould et al. (2022) focused on the sample with supergiant sources; see Figure 2.

Zoom In Zoom Out Reset image size

We compare the parameters of these events to six known FFP candidates with θ_E measurements in Table 1. The sources of all known FFP candidates except OGLE-2016-BLG-1928are red clump giants (RCGs) or red supergiants that have large θ_* = 5.4, 7.1, 11.9, 15.1, 19.5, and 34.9 μas. The magnification tends to be suppressed by large θ_* with small θ_E, i.e., large ρ as ${A}_{\mathrm{FS},\max }=\sqrt{1+4/{\rho }^{2}}$ (ρ > 1) (Maeder 1973; Agol 2003; Riffeser et al. 2006). For example, in the case of the terrestrial-mass lens with θ_E ∼ 1 μas, the maximum magnification will be only ${A}_{\mathrm{FS},\ \max }=1.066,1.039,1.014,1.009,1.005$, and 1.002 for the above values of θ_*, respectively. Note that the source of the terrestrial FFP candidate event OGLE-2016-BLG-1928S is a subgiant with θ_* = 2.37 μas. It is important to search for short FSPL with subgiants and dwarf sources to find low-mass FFPs. There is no FSPL event with a red supergiant source in our sample because these are saturated in MOA image data.

Table 1. Comparison of Parameters of Short FS Events with Known FFP Candidates

field-chip-sub-ID	tE	ρ	Is,0	θ*	θE	Reference
	(days)		(mag)	(μas)	(μas)
MOA-9y-5919	0.057 ± 0.016	1.40 ± 0.46	17.23	1.26 ± 0.48	0.90 ± 0.14	Koshimoto et al. (2023)
MOA-9y-770	0.315 ± 0.017	1.08 ± 0.07	14.71	5.13 ± 0.86	4.73 ± 0.75	Koshimoto et al. (2023)
OGLE-2016-BLG-1928	0.0288 ${}_{-0.0016}^{+0.0024}$	${3.39}_{-0.11}^{+0.10}$	15.78	2.85 ± 0.20	0.842 ± 0.064	Mróz et al. 2020
KMT-2019-BLG-2073	0.272 ± 0.007	1.138 ± 0.012	14.45	5.43 ± 0.17	4.77 ± 0.19	Kim et al. (2021)
KMT-2017-BLG-2820	0.288 ± 0.015	1.096 ± 0.079	14.31	7.05 ± 0.44	5.94 ± 0.37	Ryu et al. (2021)
OGLE-2012-BLG-1323	0.155 ± 0.005	5.03 ± 0.07	14.09	11.9 ± 0.5	2.37 ± 0.10	Mróz et al. (2019)
OGLE-2016-BLG-1540	0.320 ± 0.003	1.65 ± 0.01	13.51	15.1 ± 0.8	9.2 ± 0.5	Mróz et al. (2018)
OGLE-2019-BLG-0551	0.381 ± 0.017	4.49 ± 0.15	12.61	19.5 ± 1.6	4.35 ± 0.34	Mróz et al. (2020)
MOA-9y-1944 a	1.594 ± 0.136	0.00928 ± 0.00032	20.14	0.43 ± 0.10	46.1 ± 10.5	Koshimoto et al. (2023)
OGLE-2017-BLG-0560 a	0.905 ± 0.005	0.901 ± 0.005	12.47	34.9 ± 1.5	38.7 ± 1.6	Mróz et al. (2019)

Notes.

^a Likely BD lens.

Download table as:  ASCIITypeset image

In the final sample of Koshimoto et al. (2023), there are 10 (12) short-timescale events with t_E < 1 day after applying CR2 (CR1). Figure 3 shows the t_E distribution of the CR2 sample. The distribution is roughly symmetric in $\mathrm{log}{t}_{{\rm{E}}}$, with a tail at t_E < 0.5. This confirmed the existence of such short-timescale events with t_E < 0.5 days as reported by Mróz et al. (2017). In this section, we perform a likelihood analysis on each of the 3554 (CR1) and 3535 (CR2) events using a Galactic model to constrain the MF of lens objects.

Zoom In Zoom Out Reset image size

We define the likelihood, ${ \mathcal L }$, in Section 4.1. In Sections 4.2 and 4.3, we determine the MF without and with a planetary-mass population, respectively, by minimizing ${\chi }^{2}\equiv -2\mathrm{ln}{ \mathcal L }$. Although the absolute value of χ² is not meaningful owing to its dependence on an arbitrary normalization associated with our likelihood calculation, the fitting procedure is still statistically valid, as the relative likelihood between two models, represented by Δχ², is independent of the normalization.

Note that results of the likelihood analysis for samples CR1 and CR2 are very similar. In the following sections, we show only the results for CR2 as our final results except in the tables.

4.1. Likelihood

Although our sample contains more than 3500 events, the MF of planetary-mass objects is largely determined by the events with t_E < 1 day, which account for about 0.3% of these events. We define two likelihoods: ${{ \mathcal L }}_{\mathrm{short}}$ for short-timescale events with the best-fit t_E < 1 day, and ${{ \mathcal L }}_{\mathrm{long}}$ for events with the best-fit t_E ≥ 1 day. In our likelihood analysis, we use the combined likelihood ${ \mathcal L }={{ \mathcal L }}_{\mathrm{short}}{{ \mathcal L }}_{\mathrm{long}}$.

For ${{ \mathcal L }}_{\mathrm{long}}$, we simply use the best-fit t_E values provided by Koshimoto et al. (2023), which is similar to the approach by previous studies (Sumi et al. 2011; Mróz et al. 2017). This is because of (i) the relatively small uncertainties in t_E, (ii) the fact that the effect of individual t_E uncertainties is statistically marginalized by the large number of events, (iii) the limited sensitivity to θ_E, and (iv) the minimal impact on our primary goal of measuring the MF of planetary-mass objects.

On the other hand, the situation is the opposite for the short events, ${{ \mathcal L }}_{\mathrm{short}}$. That is, (i) the t_E uncertainties are relatively large owing to their shorter magnification period, but they must be smaller than the event selection threshold listed in Table 2 of Koshimoto et al. (2023); (ii) the number of events is very limited (12 for CR1 and 10 for CR2) and the t_E < 1 day range is only sparsely covered in Figure 3, and thus the number of t_E < 1 day events may not be sufficient to statistically marginalize the effect of t_E uncertainties of individual events in the likelihood analysis; (iii) because the ρ = θ_*/θ_E values are generally much larger than those of longer-timescale events, one may get beneficial constraints on θ_E even when the θ_E values are not well determined; and (iv) they play a crucial role in determining the MF of planetary-mass objects. Therefore, we must use the joint probability distribution of (t_E, θ_E) for each event derived by Koshimoto et al. (2023) using the Markov Chain Monte Carlo (MCMC) method for ${{ \mathcal L }}_{\mathrm{short}}$. However, the (t_E, θ_E) probability distributions for each event depend on the FFP MF that we are trying to measure, while the event detection efficiency also depends on both t_E and θ_E. Hence, the probability distribution for the t_E and θ_E values for each event depends on both the light-curve data and the FFP MF. Rather than running our light-curve model MCMC calculations for the short events separately for every MF model we consider, we simplify our calculations by using the "importance sampling" method of Monte Carlo integration (Press et al. 1992). This means that we run the MCMC light-curve models with weighting of the $\mathrm{log}{t}_{{\rm{E}}}$ and $\mathrm{log}{\theta }_{{\rm{E}}}$ distributions given by an uninformative (and incorrect) "prior" ${p}_{0}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})$ that is uniform in both $\mathrm{log}{t}_{{\rm{E}}}$ and $\mathrm{log}{\theta }_{{\rm{E}}}$. A function like ${p}_{0}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})$ is sometimes called an "interim prior" (Foreman-Mackey et al. 2014), but we have not used it as a Bayesian prior. Instead, we replace ${p}_{0}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})$ with the correct distribution over $\mathrm{log}{t}_{{\rm{E}}}$ and $\mathrm{log}{\theta }_{{\rm{E}}}$ for each MF model in our FFP MF likelihood calculation. The only Bayesian prior assumptions assumed in this analysis are the Galactic model assumptions discussed in Section 4.2 and the MF model priors discussed in Section 4.3.

We describe the simpler likelihood function for the long-duration events, ${{ \mathcal L }}_{\mathrm{long}}$, in Section 4.1.1, and then we describe ${{ \mathcal L }}_{\mathrm{short}}$ in Section 4.1.2.

4.1.1. Likelihood for Events with t_E ≥ 1 day

We define the likelihood for events with t_E ≥ 1 day by

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{long}}\propto \displaystyle \prod _{i=1}^{{N}_{\mathrm{long}}}{ \mathcal G }({t}_{{\rm{E}},i};{\rm{\Gamma }}),\end{eqnarray} \tag{ 4 }$$

where i runs over all the N_long events that have the best-fit t_E ≥ 1 day in our sample (N_long = 3542 for CR1 and N_long = 3525 for CR2), and t_E,i is the best-fit t_E value for the ith event given by Koshimoto et al. (2023).

The function ${ \mathcal G }({t}_{{\rm{E}}};{\rm{\Gamma }})$ is the model's detectable event rate as a function of t_E with given model event rate Γ, combined for the 20 survey fields, given by

$$\begin{eqnarray}&&{ \mathcal G }({t}_{{\rm{E}}};{\rm{\Gamma }})=\displaystyle \sum _{j}{w}_{j}\,{g}_{j}({t}_{{\rm{E}}};{{\rm{\Gamma }}}_{j}).\end{eqnarray} \tag{ 5 }$$

Here j takes field index values gb1 to gb21, except for gb6. See Table 1 of Koshimoto et al. (2023) for the location and properties of each field. The weight w_j for the jth field is given by

$$\begin{eqnarray}&&{w}_{j}=\displaystyle \sum _{k\in j}{n}_{\mathrm{RC},k}^{2}{f}_{\mathrm{LF},k},\end{eqnarray} \tag{ 6 }$$

where k indicates a 1024 pixel ×1024 pixel subframe in the jth field (k = 1, 2,..., 80), n_RC,k is the number density of RCGs in the kth subfield, f_LF,k is the fraction of stars with magnitude I < 21.4 mag in the kth subfield, and w_j is thus proportional to the expected event rate in the jth field. To calculate f_LF,k , we used a combined luminosity function that uses the OGLE-III photometry map (Szymański et al. 2011) for bright stars and the Hubble Space Telescope data by Holtzman et al. (1998) for faint stars.

The function g_j is the model's detectable event rate as a function of t_E for field j as given by

$$\begin{eqnarray}&&{g}_{j}({t}_{{\rm{E}}};{{\rm{\Gamma }}}_{j})={\tilde{\epsilon }}_{j}({t}_{{\rm{E}}};{{\rm{\Gamma }}}_{j})\,{{\rm{\Gamma }}}_{j}({t}_{{\rm{E}}}),\end{eqnarray} \tag{ 7 }$$

where $\tilde{\epsilon }({t}_{{\rm{E}}};{\rm{\Gamma }})$ is the integrated detection efficiency of the survey as a function of t_E. Koshimoto et al. (2023) demonstrated that when FS effects are important, the detection efficiency, (t_E, θ_E), is a function of two variables, t_E and θ_E. Therefore, we must integrate over θ_E to obtain the integrated detection efficiency, $\tilde{\epsilon }({t}_{{\rm{E}}};{\rm{\Gamma }})$, which now depends on the event rate and the MF of the lens objects. This gives

$$\begin{eqnarray}&&{\tilde{\epsilon }}_{j}({t}_{{\rm{E}}};{\rm{\Gamma }})={\int }_{{\theta }_{{\rm{E}}}}{\epsilon }_{j}({t}_{{\rm{E}}},{\theta }_{{\rm{E}}})\,{{\rm{\Gamma }}}_{j}({\theta }_{{\rm{E}}}| {t}_{{\rm{E}}})d{\theta }_{{\rm{E}}},\end{eqnarray} \tag{ 8 }$$

where _j (t_E, θ_E) is the detection efficiency for events with t_E and θ_E for the ith field. We use the detection efficiency _j (t_E, θ_E) estimated by the image-level simulations in Koshimoto et al. (2023) for the 20 fields of the MOA-II 9 yr survey.

We consider the model event rate as a function of t_E and (t_E, θ_E), denoted by Γ(t_E) and Γ(t_E, θ_E), respectively. These are normalized functions so that their integrations give 1, i.e., these are probability density functions of t_E and (t_E, θ_E), respectively. Γ(θ_E∣t_E) = Γ(t_E, θ_E)/Γ(t_E) is the probability density of events with θ_E given t_E. Thus, the calculation of $\tilde{\epsilon }({t}_{{\rm{E}}};{\rm{\Gamma }})$ in Equation (8) has to be done for every proposed MF during the fitting procedure because Γ(θ_E∣t_E) depends on the MF.

The function Γ_j (t_E, θ_E) for jth field can be separated from the MF (Han & Gould 1996),

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{j}({t}_{{\rm{E}}},{\theta }_{{\rm{E}}})=\displaystyle \int {\gamma }_{j}({t}_{{\rm{E}}}{M}^{-1/2},{\theta }_{{\rm{E}}}{M}^{-1/2}){\rm{\Phi }}(M)\sqrt{M}{dM},\end{eqnarray} \tag{ 9 }$$

where γ_j (t_E, θ_E) is the event rate for lenses with mass 1 M_☉ and Φ(M) is the present-day MF (expressed as dN/dM). Although substituting Equations (8) and (9) makes the calculation of g_j (t_E; Γ_j ) in Equation (7) a double integral over M and θ_E, Koshimoto et al. (2023) showed that the integration over θ_E is largely avoidable during a fitting procedure by switching the order of the integrals and calculating the integral over θ_E before the fitting.

We calculate γ_j (t_E, θ_E) for each field using the density and velocity distribution of stars from the latest parametric Galactic model toward the Galactic bulge based on Gaia and microlensing data (Koshimoto et al. 2021a).

Figure 4 shows the integrated detection efficiencies $\tilde{\epsilon }({t}_{{\rm{E}}};{\rm{\Gamma }})$ for the event rate calculated with the best-fit MF model with the criteria CR2. The curve for CR1 is similar. This detection efficiency is about a factor two lower than that of Mróz et al. (2017) at the low end around t_E = 0.1 days even for the similar cadence of the survey. The main reason is likely that Mróz et al. (2017) did not include the FS effect in their simulation. Koshimoto et al. (2023) confirmed that this difference is about a factor two at t_E = 0.1 days by the simulation in their Figure 7.

Zoom In Zoom Out Reset image size

Note that detection efficiencies at short t_E with ρ > 1 may be improved in a future analysis. For events with ρ > 1, the magnification can be significant with the minimum impact parameter up to u₀ ≲ ρ. This is likely to be more important for bright giant sources because these have higher signal-to-noise ratio even at low magnification (see also the Appendix). However, such events are rejected by the criterion u₀ ≤ 1 in Koshimoto et al. (2023). This criterion is applied because it is useful to robustly remove the various artifacts and keep the sample as clean as possible. This may be improved in a future analysis with a more careful investigation.

4.1.2. Likelihood for Short-timescale (t_E < 1 day) Events

We follow the importance sampling method used by Hogg et al. (2010) to convert the uninformative "interim prior" ${p}_{0}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})$ used for the light-curve MCMC for each event to a probability distribution for event i, ${ \mathcal G }(\mathrm{log}{t}_{{\rm{E}},i},\mathrm{log}{\theta }_{{\rm{E}},i};{\rm{\Gamma }})$, that depends on the event rate for each MF model, Γ. However, while Hogg et al. (2010) characterized their calculation as a modification of the assumed prior, based on the data, this is not the case for our analysis. Instead, we are replacing p₀ with the probability distribution implied by our MF model, using the importance sampling Monte Carlo integration method (Press et al. 1992). We use the probability distribution for each event from its MCMC analysis to calculate the likelihood for the short-timescale events, ${{ \mathcal L }}_{\mathrm{short}}$. Given the output MCMC samples of posterior distributions for individual events by Koshimoto et al. (2023), the likelihood is given by

$$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{short}}\propto \displaystyle \prod _{i=1}^{{N}_{\mathrm{short}}}\displaystyle \sum _{k=1}^{{K}_{i}}\displaystyle \frac{{ \mathcal G }(\mathrm{log}{t}_{{\rm{E}},{ik}},\mathrm{log}{\theta }_{{\rm{E}},{ik}};{\rm{\Gamma }})}{{p}_{0}(\mathrm{log}{t}_{{\rm{E}},{ik}},\mathrm{log}{\theta }_{{\rm{E}},{ik}})},\end{eqnarray} \tag{ 10 }$$

where i runs over all the N_short events that have the best-fit t_E < 1 day (N_short = 12 for CR1 and N_short = 10 for CR2), k runs over all the K_i samples in the MCMC sample of the probability distribution for the ith event, and ${p}_{0}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})$ is the uninformative prior distribution used for these MCMC calculations. The model's detectable event rate as a function of $(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})$ is given by

$$\begin{eqnarray}&&{ \mathcal G }(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}};{\rm{\Gamma }})=\displaystyle \sum _{j}{w}_{j}\,{g}_{j}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}};{{\rm{\Gamma }}}_{j})\end{eqnarray} \tag{ 11 }$$

with

$$\begin{eqnarray}&&{g}_{j}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}};{{\rm{\Gamma }}}_{j})={\epsilon }_{j}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})\,{{\rm{\Gamma }}}_{j}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}}),\end{eqnarray} \tag{ 12 }$$

where we represented it as a function of $(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})$ rather than (t_E, θ_E) because the MCMC calculations of Koshimoto et al. (2023) provide the probability distributions based on the uninformative uniform prior in $(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})$, i.e., ${p}_{0}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})=\mathrm{const}.$.

Equation (10) calculates the likelihood by summing the ratio of ${ \mathcal G }(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}};{\rm{\Gamma }})$ to ${p}_{0}(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})$ to replace the uniform prior (i.e., p₀), used for the MCMC calculations with the new probability distribution (i.e., ${ \mathcal G }$) that depends on our MF model. This method, which uses all the MCMC samples, allows ${{ \mathcal L }}_{\mathrm{short}}$ to account for the uncertainty of the parameters, unlike ${{ \mathcal L }}_{\mathrm{long}}$ given in Equation (4).

Despite the significant computational cost of Equation (10) associated with performing a summation over K_i (typically ∼5 × 10⁵) samples for each proposed MF during the fitting process, we addressed this by implementing a binning strategy for the MCMC sample using grids of $(\mathrm{log}{t}_{{\rm{E}}},\mathrm{log}{\theta }_{{\rm{E}}})$ with a size of (0.05 dex × 0.05 dex), which significantly increased the computational efficiency.

4.2. Mass Function of Known Population

First, we perform the likelihood analysis without the short events with t_E < 1 day using the Galactic model with the MF of known populations, i.e., stellar remnants (black holes (BHs), neutron stars (NSs), and white dwarfs (WDs)), MS stars, and BDs. We use a broken power-law MF given by

$$\begin{eqnarray}\displaystyle \frac{{dN}}{d\mathrm{log}M}\propto \left\{\begin{array}{ll}{M}^{-{\alpha }_{1}} & ({M}_{1}\lt M/{M}_{\odot }\lt 120)\\ {M}^{-{\alpha }_{2}} & (0.08\lt M/{M}_{\odot }\lt {M}_{1})\\ {M}^{-{\alpha }_{3}} & (3\times {10}^{-4}\lt M/{M}_{\odot }\lt 0.08).\end{array}\right.\end{eqnarray} \tag{ 13 }$$

We adopt the values of parameters α₁ = 1.32, α₂ = 0.13, α₃ = −0.82, and M₁=0.86 from the E + E_X model of Koshimoto et al. (2021a) by default unless specified as fitting parameters in the following three models. The minimum mass 3 × 10⁻⁴ M_☉ is taken to be smaller than the theoretical minimum mass of the gas cloud, ∼Jupiter mass, that collapses to form a BD (Boss et al. 2003). During our fitting procedure, a proposed initial mass function (IMF) is converted into a present-day MF following the procedure used by Koshimoto et al. (2021a) that combines their stellar age distribution and the initial−final mass relation by Lam et al. (2020) to evolve stars into stellar remnants.

We consider three models here: BD1, BD2, and BD3. In BD1, we fit only α₃ as a fitting parameter, while fixing α₁, α₂, and M₁. Similarly, in BD2, we fit α₃ and α₂, and in BD3, we fit α₁, α₂, α₃, and M₁. To perform the fitting, we use the MCMC method (Metropolis et al. 1953) and assign uniform distributions as priors for all the parameters.

The best-fit models BD1, BD2, and BD3 are almost indistinguishable from the blue dotted line in Figure 3. One can see that the models fit the data with t_E > 1 day very well. The best-fit parameters and χ² values are listed in Table 2. There is no significant difference in the resultant parameters between different selection criteria or among the BD1, BD2, and BD3 models.

Table 2. Best-fit Parameters of the Mass Function for Known Populations

Model	BD1		BD2		BD3		Koshimoto et al. (2021a) a
	CR1	CR2	CR1	CR2	CR1	CR2
M1	(0.86)	(0.86)	(0.86)	(0.86)	${0.97}_{-0.34}^{-0.04}$	${0.99}_{-0.37}^{-0.06}$	${0.86}_{-0.10}^{+0.09}$
α1	(1.32)	(1.32)	(1.32)	(1.32)	${1.33}_{-0.17}^{+0.21}$	${1.34}_{-0.18}^{+0.18}$	${1.32}_{-0.10}^{+0.14}$
α2	(0.13)	(0.13)	${0.20}_{-0.05}^{+0.07}$	${0.20}_{-0.05}^{+0.07}$	${0.23}_{-0.19}^{+0.04}$	${0.24}_{-0.21}^{+0.04}$	${0.13}_{-0.12}^{+0.11}$
α3	$-{0.60}_{-0.13}^{+0.08}$	$-{0.62}_{-0.14}^{+0.09}$	$-{0.74}_{-0.30}^{+0.13}$	$-{0.76}_{-0.30}^{+0.14}$	$-{0.76}_{-0.26}^{+0.19}$	$-{0.79}_{-0.25}^{+0.22}$	$-{0.82}_{-0.51}^{+0.24}$
χ2	35,919.4	35,722.6	35,918.2	35,721.5	35,918.0	35,721.3

Notes. Some of the upper errors of M₁ are negative because the best-fit value is outside of the 68% range. This is because M₁ is restricted to be less than 1 M_☉.

^a Results of fitting to various bulge data, including the OGLE-IV t_E distribution of t_E > 1 day (Mróz et al. 2017, 2019). The representative values are shifted to the ones for the E + E_X model from their original ones for the G + G_X model.

Download table as:  ASCIITypeset image

All of the parameters are consistent with those of Koshimoto et al. (2021a) within 1σ. This indicates that our data set confirmed the Galactic model and MF of known objects by Koshimoto et al. (2021a). This also indicates that our data set is consistent with the OGLE-IV t_E distribution for t_E > 1 day (Mróz et al. 2017, 2019) that is fitted by Koshimoto et al. (2021a).

In the following analysis, we fit only α₃ and fix all other parameters for the known populations. Note that, in Koshimoto et al. (2021a), the Galactic model and MF are constrained to satisfy the microlensing t_E distribution, stellar number counts, and the Galactic bulge mass from other observations simultaneously. In principle, the MF should not be changed alone because it is related to other parameters of the Galactic model. However, the contribution of objects with M/M_☉ < 0.08 is negligible in stellar number counts and as a fraction of the Galactic bulge mass. Thus, we assume that a model with a different slope at lower masses with M/M_☉ < 0.08 is still valid.

4.3. Mass Function of Planetary-mass Population

If the candidates with t_E < 0.5 days are really due to microlensing, they cannot be explained by known populations, i.e., stellar remnants, MS stars, or BDs. To explain the tail for short values of t_E, we defined a new model "PL" that introduces a planetary-mass population by the following power law in addition to known populations (Equation (13)):

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{4}}{d\mathrm{log}M}=Z{\left(\displaystyle \frac{M}{{M}_{\mathrm{norm}}}\right)}^{-{\alpha }_{4}},({M}_{\min }\lt M{/M}_{\odot }\lt 0.02).\end{eqnarray} \tag{ 14 }$$

Here Z is a normalization factor and M_norm is a reference mass whose inclusion allows Z to have a unit of (dex)⁻¹. Although M_norm can be an arbitrary zero-point, we found that the uncertainty in Z is minimized when we adopt M_norm = 8 M_⊕, which is recognized as a pivot point.

In the model PL, we use α₃, α₄, and Z as fitting parameters and fix parameters α₁ = 1.32, α₂ = 0.13, and M₁ = 0.86 (Koshimoto et al. 2021a). We assign uniform distributions as priors for α₃, α₄, and $\mathrm{log}Z$ in our MCMC run. We found that the fitting result does not depend on ${M}_{\min }$ at all when ${M}_{\min }\lt 3\times {10}^{-7}\,{M}_{\odot }$, which indicates that our data sensitivity is down to ∼3 × 10⁻⁷ M_☉. Thus, we decided to use ${M}_{\min }={10}^{-7}$.

The red solid line in Figure 3 represents the best-fit model for all populations with the CR2 sample. This figure indicates that the model represents the observed t_E distribution well. Note that although the observed t_E distribution shown in black in Figure 3 does not include error bars along the t_E axis, the best-fit line is derived from our likelihood analysis that takes into account the t_E errors, as well as the θ_E constraints for the short events with t_E < 1 day. Figure 5 shows the posterior distributions of the parameters of PL model. The best-fit parameters and χ² are listed in Table 3.

Zoom In Zoom Out Reset image size

Table 3. Best-fit Parameters of the Mass Function for the Planetary-mass Population

	CR1	CR2		Gould et al. (2022)
(Mnorm)	(8 M⊕)	(8 M⊕)	(38 M⊕)	(38 M⊕)
M1	(0.86)	(0.86)
α1	(1.32)	(1.32)
α2	(0.13)	(0.13)
α3	$-{0.55}_{-0.17}^{+0.13}$	$-{0.58}_{-0.16}^{+0.12}$
α4	${0.90}_{-0.27}^{+0.48}$	${0.96}_{-0.27}^{+0.47}$		Fixed at 0.9 or 1.2
Z	${2.08}_{-1.33}^{+0.54}$	${2.18}_{-1.40}^{+0.52}$	${0.49}_{-0.37}^{+0.17}$
ZMS+BD	${2.27}_{-1.46}^{+0.60}$	${2.38}_{-1.53}^{+0.58}$	${0.53}_{-0.40}^{+0.19}$	0.39 ± 0.20 ± ?
${Z}^{{M}_{\odot }}$	${5.33}_{-3.40}^{+1.26}$	${5.48}_{-3.50}^{+1.18}$	${1.22}_{-0.91}^{+0.35}$
${Z}_{\mathrm{MS}+\mathrm{BD}}^{{M}_{\odot }}$	${10.63}_{-6.78}^{+2.52}$	${10.95}_{-6.97}^{+2.36}$	${2.44}_{-1.82}^{+0.71}$	1.96 ± 0.98 ± ?
f a	${17}_{-11}^{+20}$	${21}_{-13}^{+23}$
fMS+BD a	${19}_{-12}^{+22}$	${23}_{-15}^{+25}$
${f}^{{M}_{\odot }}$ a	${45}_{-30}^{+54}$	${53}_{-34}^{+59}$
${f}_{\mathrm{MS}+\mathrm{BD}}^{{M}_{\odot }}$ a	${89}_{-59}^{+107}$	${106}_{-68}^{+117}$
m b	${89}_{-56}^{+96}{M}_{\oplus }$	${80}_{-47}^{+73}{M}_{\oplus }$
mMS+BD b	${98}_{-61}^{+107}{M}_{\oplus }$	${88}_{-51}^{+81}{M}_{\oplus }$
${m}^{{M}_{\odot }}$ b	${229}_{-140}^{+219}{M}_{\oplus }$	${202}_{-114}^{+166}{M}_{\oplus }$
${m}_{\mathrm{MS}+\mathrm{BD}}^{{M}_{\odot }}$ b	${457}_{-279}^{+439}{M}_{\oplus }$	${404}_{-228}^{+333}{M}_{\oplus }$
χ2	36,273.0	36,024.1

Notes. We adopt the model for CR2 as the final result.

^a Number of planetary-mass objects per BD+MS+WD (f), per MS+BD (f_MS+BD), per solar mass of BD+MS+WD (${f}^{{M}_{\odot }}$), or per solar mass of MS+BD (${f}_{\mathrm{MS}+\mathrm{BD}}^{{M}_{\odot }}$) when MFs down to 10⁻⁶ M_☉ are integrated. These vary depending on the minimum mass. ^b Total mass of planetary-mass objects per BD+MS+WD (m), per MS+BD (m_MS+BD), per solar mass of BD+MS+WD (${m}^{{M}_{\odot }}$), or per solar mass of MS+BD (${m}_{\mathrm{MS}+\mathrm{BD}}^{{M}_{\odot }}$) when MFs down to 10⁻⁶ M_☉ are integrated.

Download table as:  ASCIITypeset image

The best-fit power-law index for BD is ${\alpha }_{3}=-{0.58}_{-0.16}^{+0.12}$, which is consistent with the model without the planetary-mass population.

The best-fit MF of the planetary-mass populations with the normalization Z relative to stars (MS+BD+WD) (integrated IMF over 3 × 10⁻⁴ < M/M_☉ < 8) can be expressed as

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{4}}{d\mathrm{log}M}=\displaystyle \frac{{2.18}_{-1.40}^{+0.52}}{\mathrm{dex}\times \mathrm{star}}{\left(\displaystyle \frac{M}{8\,{M}_{\oplus }}\right)}^{-{\alpha }_{4}},\end{eqnarray} \tag{ 15 }$$

where ${\alpha }_{4}={0.96}_{-0.27}^{+0.47}$. Figure 6 shows the IMF of the best-fit PL model. This α₄ is consistent with the corresponding power-law index of 0.9 ≲ p ≲ 1.2 suggested by Gould et al. (2022).

Zoom In Zoom Out Reset image size

This can be translated to the normalization per stellar mass of stars, ${Z}^{{M}_{\odot }}$, as

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{4}}{d\mathrm{log}M}=\displaystyle \frac{{5.48}_{-3.50}^{+1.18}}{\mathrm{dex}\times {M}_{\odot }}{\left(\displaystyle \frac{M}{8\,{M}_{\oplus }}\right)}^{-{\alpha }_{4}}.\end{eqnarray} \tag{ 16 }$$

This implies that the number of FFPs per star is $f\,={21}_{-13}^{+23}$ star⁻¹ over the mass range 10⁻⁶ < M/M_☉ < 0.02 (0.33 < M/M_⊕ < 6660). Note that this value varies depending on the minimum mass. The total mass of FFPs per star is $m={80}_{-47}^{+73}{M}_{\oplus }({0.25}_{-0.15}^{+0.23}{M}_{{\rm{J}}})$ star⁻¹. This is less dependent on the minimum mass. The total mass of FFPs per M_☉ is ${m}^{{M}_{\odot }}={202}_{-114}^{+166}$ ${M}_{\oplus }({0.64}_{-0.11}^{+0.19}{M}_{{\rm{J}}}){M}_{\odot }^{-1}$. This is a more robust value that is less dependent on uncertainty in the abundances of the low-mass objects for both FFP and BD.

The normalization, number, and total mass of FFPs relative to MS+BD (3 × 10⁻⁴ < M/M_☉ < 1.1) are also shown in Table 3. These normalizations can be translated to ${Z}_{\mathrm{MS}+\mathrm{BD}}={0.53}_{-0.40}^{+0.19}$ dex⁻¹ star⁻¹ and ${Z}_{\mathrm{MS}+\mathrm{BD}}^{{M}_{\odot }}={2.44}_{-1.82}^{+0.71}$ dex⁻¹ ${}^{}{M}_{\odot }^{-1}$ with M_norm = 38 M_⊕. These are almost the same as Z_MS+BD = 0.39 ± 0.18 dex⁻¹ star⁻¹ and Z_MS+BD = 1.96 ± 0.98 dex⁻¹ ${}^{}{M}_{\odot }^{-1}$ with M_norm = 38 M_⊕ by Gould et al. (2022).

Note that the lenses for these short events could be either FFPs or planets with very wide separations of more than about 10 au from their host stars, for which we cannot detect the host star in the light curves.

4.4. Broken Power-law Mass Function for the Planetary-mass Population

In order to demonstrate the FFP MF uncertainty at low masses, we have also modeled the planetary-mass population with a broken power-law MF given by

$$\begin{eqnarray}\displaystyle \frac{{{dN}}_{4}}{d\mathrm{log}M}\propto \left\{\begin{array}{ll}Z{\left(\displaystyle \frac{M}{{M}_{\mathrm{norm}}}\right)}^{-{\alpha }_{4}}, & ({M}_{\mathrm{br}}\lt M/{M}_{\odot }\lt 0.02)\\ {M}^{-{\alpha }_{5}}, & ({M}_{\min }\lt M/{M}_{\odot }\lt {M}_{\mathrm{br}}).\end{array}\right.\end{eqnarray} \tag{ 17 }$$

Here M_br is a break mass and α₅ is a power below M_br. ${M}_{\min }={10}^{-7}\,{M}_{\odot }$ is the same as in the previous section.

In Figure 6, we show the 1σ range of the broken power-law PL model along with the best-fit single power-law MF given in the previous section for comparison. The median and 1σ range of the parameters and χ² are listed in Table 4. The resultant broken power-law MF is consistent with the single power-law model, while the uncertainty is larger. Although the MF is relatively well constrained down to an Earth mass, the uncertainty is much larger below an Earth mass. This is as expected because of our low sensitivity below to planets of less than an Earth mass.

Table 4. Median and Uncertainty of Parameters of the Broken Power-law Mass Function for the Planetary-mass Population

	CR1	CR2
(Mnorm)	(8 M⊕)	(8 M⊕)
α3	$-{0.54}_{-0.17}^{+0.12}$	$-{0.58}_{-0.19}^{+0.12}$
α4	${1.07}_{-0.49}^{+0.93}$	${1.14}_{-0.54}^{+0.97}$
α5	${0.13}_{-3.07}^{+1.33}$	${0.13}_{-3.10}^{+1.32}$
$\mathrm{log}{M}_{\mathrm{br}}$	$-{5.35}_{-1.02}^{+1.35}$	$-{5.27}_{-1.05}^{+1.28}$
Z	${1.79}_{-1.08}^{+2.91}$	${1.85}_{-1.17}^{+3.14}$
ZMS+BD	${1.96}_{-1.18}^{+3.19}$	${2.03}_{-1.28}^{+3.43}$
${Z}^{{M}_{\odot }}$	${4.57}_{-2.77}^{+7.54}$	${4.62}_{-2.91}^{+7.92}$
${Z}_{\mathrm{MS}+\mathrm{BD}}^{{M}_{\odot }}$	${9.12}_{-5.52}^{+15.05}$	${9.22}_{-5.82}^{+15.79}$
f a	${15}_{-11}^{+36}$	${17}_{-12}^{+39}$
fMS+BD a	${17}_{-12}^{+40}$	${18}_{-13}^{+42}$
${f}^{{M}_{\odot }}$ a	${39}_{-28}^{+96}$	${42}_{-30}^{+98}$
${f}_{\mathrm{MS}+\mathrm{BD}}^{{M}_{\odot }}$ a	${79}_{-57}^{+191}$	${85}_{-61}^{+196}$
m	${73}_{-40}^{+119}$	${69}_{-36}^{+107}$
mMS+BD b	${80}_{-44}^{+131}$	${75}_{-39}^{+118}$
${m}^{{M}_{\odot }}$ b	${192}_{-103}^{+275}$	${175}_{-89}^{+246}$
${m}_{\mathrm{MS}+\mathrm{BD}}^{{M}_{\odot }}$ b	${384}_{-206}^{+551}$	${349}_{-178}^{+493}$
χ2	36,271.6	36,022.9

Notes. The median and 1σ ranges are shown for understanding the uncertainty.

^a Same as Table 3. ^b Same as Table 3.

Download table as:  ASCIITypeset image

This model implies that the number of FFPs per star is $f={17}_{-12}^{+39}$ star⁻¹ over the mass range 10⁻⁶ < M/M_☉ < 0.02 (0.33 < M/M_⊕ < 6660). The total mass of FFPs per star is $m={69}_{-36}^{+107}$ ${M}_{\oplus }({0.22}_{-0.11}^{+0.33}{M}_{{\rm{J}}})$ star⁻¹. The total mass of FFPs per M_☉ is ${m}^{{M}_{\odot }}={175}_{-89}^{+246}$ ${M}_{\oplus }\,({0.55}_{-0.28}^{+0.77}{M}_{{\rm{J}}}){M}_{\odot }^{-1}$. These numbers are also consistent with those for the single power-law model but have larger uncertainties. This result is useful to see the conservative uncertainty of the MF. In the following discussion, although we use only the results for the single power-law model, the discussion is qualitatively the same for the broken power law.

4.5. Comparison to Sumi et al. (2011)

We derived the MF of lens objects from the 9 yr MOA-II survey toward the Galactic bulge. The 3535 high-quality single-lens light curves used in our statistical analysis include 10 very short (t_E < 1 day) events and 13 events with strong FS effects that allow the determination of the angular Einstein radius, θ_E.

The cumulative θ_E histogram for these 13 events reveals an "Einstein gap" at 5 < θ_E/μas < 70, which is roughly consistent with the gap at 10 < θ_E/μas < 30 found by the KMTNet group (Ryu et al. 2021; Gould et al. 2022). This gap indicates that there is a distinct planetary-mass population separated from the known populations of BDs, stars, and stellar remnants.

We constructed the t_E distribution of all selected samples including both PSPL and FSPL. We calculated the integrated detection efficiency $\tilde{\epsilon }({t}_{{\rm{E}}};{\rm{\Gamma }})$ of the survey by integrating the two-dimensional detection efficiency, (t_E, θ_E), measured from image-level simulations that included the FS effect, and convolving this with the event rate Γ(t_E, θ_E) given by a Galactic model and MF. We found that the t_E distribution has an excess at short t_E values that cannot be explained by known populations.

We then adopted the single power-law MF for the planetary-mass population. We found that these short events can be well modeled by ${{dN}}_{4}/d\mathrm{log}M={({2.18}_{-1.40}^{+0.52})\times (M/8\,{M}_{\oplus })}^{-{\alpha }_{4}}$ dex⁻¹ star⁻¹ with ${\alpha }_{4}={0.96}_{-0.27}^{+0.47}$ at 10⁻⁷ < M/M_⊙ < 0.02 (or 0.033 < M/M_⊕ < 6660).

This can also be expressed by the MF per stellar mass as ${{dN}}_{4}/d\mathrm{log}M={5.48}_{-3.50}^{+1.18}\times {(M/8\,{M}_{\oplus })}^{-{\alpha }_{4}}$ dex⁻¹ ${}^{}{M}_{\odot }^{-1}$. We showed that the number of FFPs or distant planets is $f\,={21}_{-13}^{+23}$ star⁻¹. Note that we found $f={17}_{-12}^{+39}$ FFPs per star for the broken power-law model, which is consistent with our result for the single power-law model, with a larger uncertainty. In the following discussion, we only use the results for the single power-law model; the conclusions are qualitatively the same for the broken power-law model.

It is well known that planet–planet scattering during the planet formation process is likely to produce a population of unbound or wide-orbit planetary-mass objects (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997). The probability of planet scattering likely increases with declining mass because planets usually require more massive planets to scatter. Hence, we expect the power-law index of the MF of bound planets α _b to be smaller than that of α₄ for unbound or large-orbit planets, i.e., α₄ > α _b .

One can compare our FFP result to the MF of known bound planets. At present, microlensing surveys have only measured the mass ratio function, rather than the MF, of the bound planets. Currently, the most sensitive study of the bound planet mass ratio function, Suzuki et al. (2016), found that the mass ratio function can be well explained by the broken power law with α _b = 0.93 ± 0.13 for q > q_br = 1.7 × 10⁻⁴ and ${\alpha }_{b}=-{0.6}_{-0.5}^{+0.4}$ for q < q_br = 1.7 × 10⁻⁴. While the Suzuki et al. (2016) data could establish the existence of the power-law break with reasonably high confidence (a Bayes factor of 21), there was a large, correlated uncertainty in the mass ratio of the break and slope of the mass ratio function below the break. Hence, we chose to fix the mass ratio of the break at q_br = 1.7 × 10⁻⁴ in order to estimate the power law below the break.

More recently, several papers have attempted to improve on this estimate by including a heterogeneous set of lower mass ratio planets found by a number of groups without a calculation of the detection efficiency. These efforts included attempts to estimate the effect of a "publication bias" that might cause planets deemed to be of greater interest to be published much more quickly, leading to a biased, inhomogeneous sample of planets. This "publication bias" is caused by the decision to publish some planet discoveries at a higher priority than others. With such an analysis, Udalski et al. (2018) reported ${\alpha }_{b}=-{1.05}_{-0.78}^{+0.68}$ with their sample and ${\alpha }_{b}=-{0.73}_{-0.34}^{+0.42}$ when combined with the Suzuki et al. (2016) result for q < 1 × 10⁻⁴ < q_br. A similar analysis by Jung et al. (2019) attempted a new measurement of the location of the break and found α _b = −4.5 for q < q_br = 0.55 × 10⁻⁴, which is consistent with the Suzuki et al. (2016) result when q_br is not fixed. However, a more recent paper (Zang et al. 2022) by many of the same authors reported a number of planetary microlensing events that were missed by the analyses described in Udalski et al. (2018) and Jung et al. (2019). This casts some doubt on the validity of some of the assumptions in these papers. This later paper also suggests that planets with mass ratios of q < q_br = 1.7 × 10⁻⁴ may be more common than previously thought, although a more definitive claim awaits a detection efficiency calculation. In addition, the Suzuki et al. (2016) analysis does not imply that there is a peak in the mass ratio. Instead, it concludes that the slope does not rise as steeply toward low mass ratios as it does for q > 1.7 × 10⁻⁴.

The broken power-law model of Suzuki et al. (2016) is consistent with the hypothesis that these unbound or wide-orbit planetary-mass objects are the result of scattering from bound planetary systems. It is the lower-mass planets that are preferentially removed by planet–planet scattering interactions, so the initial planetary MF may have been closer to a single power law with α _b ∼ 0.9, but planet–planet scattering has likely depleted the numbers of low-mass planets at separations beyond the snow line where microlensing is most sensitive. Thus, planet–planet scattering may be responsible for the mass ratio function "break" observed in the Suzuki et al. (2016) sample.

This idea that planet–planet scattering is responsible for a FFP MF slope that is steeper than the slope of the mass ratio function for low-mass bound planets is also consistent with the single power-law models that were found in smaller data sets (Sumi et al. 2010). The best-fit single power-law model for the Suzuki et al. (2016) sample gives ${{dN}}_{\mathrm{bound}}/d\mathrm{log}q\,={0.068}_{-0.014}^{+0.016}\,{\mathrm{dex}}^{-2}\,{\mathrm{star}}^{-1}\times {(q/0.001)}^{-{\alpha }_{{\rm{b}}}}$ with α _b = 0.58 ± 0.08 for 3 × 10⁻⁶ < q < 3 × 10⁻², but the broken power law is a significantly better fit to the Suzuki et al. (2016) data. Note that this single power-law model with α _b = 0.58 satisfies α₄ > α _b , for our value of ${\alpha }_{4}={0.96}_{-0.27}^{+0.47}$, implying that unbound (or very wide orbit) planets increase more rapidly than bound planets at low masses. Thus, our main conclusion discussed below with the broken power-law model, in which the lower-mass planets are increasingly scattered, is not specific to the Suzuki et al. (2016) broken power-law model.

As a comparison, we transformed the bound planet's mass "ratio" function of Suzuki et al. (2016) to an MF by using the estimated average mass of their hosts of ∼0.56 M_☉ as shown ¹⁵ in Figure 6. We estimate the abundance of the wide-orbit bound planets to be ${f}_{\mathrm{wide}}={1.1}_{-0.3}^{+0.6}$ planets star⁻¹ in the mass range 10⁻⁶ < M/M_☉ < 0.02 (0.33 < M/M_⊕ < 6660) and separation range 0.3 < s < 5, which corresponds to a semimajor axis of roughly 0.7 < a/au < 12. This indicates that the abundance of FFPs, $f={21}_{-13}^{+23}$ planets star⁻¹, is ${19}_{-13}^{+23}$ times more than wide-orbit bound planets in this mass range.

This is because the number of wide-orbit bound planets decreases at lower masses than the break at M_break ≈ 1.0 × 10⁻⁴ M_☉, while the number of high-mass bound planets is larger than that for FFPs. Again, this is consistent with the hypothesis that the low-mass planets are more likely to be scattered. Note that there is still large uncertainty in the MFs at low masses for both bound and unbound planets. It is very important to constrain these MFs at low masses.

We can also compare our number for the FFP abundance with the abundance of the bound planets with short-period orbits of P = 0.5–256 days and planetary radii of R _p = 0.5–4 R_⊕ found by Kepler. Hsu et al. (2019) find ${f}_{\mathrm{FGK}}={3.5}_{-0.6}^{+0.7}$ for FGK dwarfs, and Hsu et al. (2020) find ${f}_{{\rm{M}}}={4.2}_{-0.6}^{+0.6}$ for M dwarfs. Because the typical spectral types of their samples are G2 (M = 1 M_☉) and M2.5 (M = 0.4 M_☉), their typical semimajor axes are 0.012 ≲ a/au ≲ 0.79 and 0.009 ≲ a/au ≲ 0.58, respectively. The fractions of FGK and M dwarfs relative to all populations except BHs and NSs are 0.157 and 0.465 in our best-fit MF. By weighting with these stellar type fractions, the abundance of the known close-orbit bound planets is about ${f}_{\mathrm{close}}={2.5}_{-0.3}^{+0.3}$ star⁻¹. (This ignores the relatively small number of gas giant planets in short-period orbits; Bryant et al. 2023.)

The total abundance of the wide-orbit and known close-orbit bound planets is about ${f}_{\mathrm{bound}}={3.6}_{-0.4}^{+0.7}$ star⁻¹. This indicates that the abundance of FFPs, $f={21}_{-13}^{+23}$ planets star⁻¹, is ${5.8}_{-3.8}^{+6.4}$ times more than known bound planets in this mass range.

We found that the total mass of FFPs or distant planets per star is $m={80}_{-47}^{+73}{M}_{\oplus }({0.25}_{-0.15}^{+0.23}{M}_{{\rm{J}}})$ star⁻¹ in this 10⁻⁶ < M/M_☉ < 0.02 (0.33 < M/M_⊕ < 6660) mass range. This is comparable to the value of ${91}_{-22}^{+33}$ M_⊕ star⁻¹ for wide-orbit bound planets with separations of 0.3 < s < 5 in the same mass range. It is not straightforward to estimate the total mass of inner planets found by Kepler because only a small, and somewhat biased, sample of Kepler planets have mass measurements. The total masses of FFPs and bound planets are less dependent on the uncertainty of the number of low-mass planets than the total numbers of FFPs and bound planets are.

These comparisons indicate that ${19}_{-13}^{+23}$ times more planets than the ones currently in wide orbits have been ejected to unbound or very wide orbits. These comparisons also suggest that the total mass of scattered planets is of the same order as those remaining bound in wide orbits (beyond the snow line) in their planetary systems. The low-mass bound planets in wide orbits are much less abundant than those orbiting closer to their host stars. This may be explained by the fact that planets in wide orbits are more easily ejected than those in close orbits.

The power-law index of the IMF of planets formed in wide orbits in protoplanetary disks is likely to be α₄ ∼ 0.9 with an abundance of ${22}_{-13}^{+23}$ planets star⁻¹ or ${171}_{-52}^{+80}{M}_{\oplus }({0.54}_{-0.16}^{+0.25}{M}_{{\rm{J}}})$ star⁻¹.

Various formation mechanisms of FFPs from bound planetary systems have been proposed. Planets can be ejected from their hosts by dynamical interactions with other (mostly giant) planets (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997), by stellar flybys (Malmberg et al. 2011), or by the post-MS evolution of their hosts (Adams et al. 2013). Coleman et al. (2023) simulated the circumbinary planetary systems for Kepler-16 and Kepler-34 and found that such systems may eject 6.3 and 9.3 planets on average, respectively, and most of these have masses smaller than Neptune. However, there are very few or almost no studies on the prediction for the number of the ejection of the Earth-to-Neptune-mass planet population because the abundance of such planets in less tightly bound wide orbits is not well known. The results of our study may shed light on this area.

Another, rather speculative possibility is that most of the low-mass objects found by microlensing are primordial BHs (PBHs; Niikura et al. 2019, 2019). Hashino et al. (2022) predicted PBHs generated at a first-order electroweak phase transition have masses of about 10⁻⁵ M_☉. They found that, depending on parameters of the phase transition, a sufficient number of PBHs can be produced to be observed by current and future microlensing surveys. The mass of such PBHs is a function of the time of their generation, i.e., the electroweak phase transition, and is expected to be a delta-function distribution. To differentiate PBHs from FFPs, we need to measure the shape of the MF accurately. This can be done by current (MOA, OGLE, KMTNet) surveys, a near-future (PRIME) ground telescope, and the Roman Space Telescope.

For the first time, we have determined the detection efficiency as a function of both the Einstein radius crossing time and the angular Einstein radius, because FS effects have a large influence on the detectability of microlensing events owing to low-mass planets. This method is necessary for reliable results for low-mass FFPs, and it should be very useful for the analysis of these future surveys that will detect many short events.

A precise measurement of the FFP MF will require a microlensing survey that can obtain precise photometry of MS stars with relatively low magnification because the small angular Einstein radii, θ_E, of low-mass planetary lenses prevent high magnification. The exoplanet microlensing survey of the Roman Space Telescope is such a survey, and it should provide the definitive measurement of the FFP MF. Johnson et al. (2020) predicted ∼250 FFPs with masses down to that of Mars (including ∼25 with masses of 0.1 ≤ M/M_⊕ ≤ 1 and ∼48 with 0.316 ≤ M/M_⊕ ≤ 3.16) assuming the fiducial MF of cold, bound planets adapted from Cassan et al. (2012). Our FFP MF results imply a large increase in the number of FFP events that should be detected by Roman. We predict ${988}_{-566}^{+1848}$ FFPs with masses down to that of Mars (including ${575}_{-424}^{+1733}$ with 0.1 ≤ M/M_⊕ ≤ 1 and ${391}_{-259}^{+344}$ with 0.316 ≤ M/M_⊕ ≤ 3.16) for our single power-law model. The broken power-law model predicts ${699}_{-418}^{+1424}$ FFPs down to that of Mars (including ${303}_{-271}^{+1268}$ with 0.1 ≤ M/M_⊕ ≤ 1 and ${261}_{-213}^{+436}$ with 0.316 ≤ M/M_⊕ ≤ 3.16).

The Earth 2.0 (ET) mission is a proposed space telescope to conduct the transit and microlensing exoplanet surveys. One of seven 30 cm telescopes will be used for the microlensing survey toward the GB. The ET is planning to measure the masses of FFPs by the space parallax in collaboration with ground-based telescopes. Ge et al. (2022) estimated that ET will detect about 600 FFP events, of which about 150 will have mass measurements. Our MF is about a factor of 1.4 higher normalization than that assumed in Ge et al. (2022) with similar slope. However, they assumed a flat MF for ≤1 M_⊕, while we continued the power-law slope down to the lower limit of 0.1 M_⊕. This renormalization will update the expected yield of ∼840 FFPs with masses down to that of Mars (including ∼210 with masses ≤M_⊕).

We thank the anonymous referee for the useful suggestions. We are grateful to S. Ida, S. A. Johnson, and K. Masuda for helpful comments. The MOA project is supported by JSPS KAKENHI grant Nos. JP24253004, JP26247023, JP23340064, JP16H06287, JP17H02871, and JP22H00153. N.K. was supported by the JSPS overseas research fellowship. D.P.B. acknowledges support from NASA grants 80NSSC20K0886 and 80NSSC18K0793.

Koshimoto et al. (2023) calculated the integrated detection efficiency for our FSPL event sample, ${\tilde{\epsilon }}_{\mathrm{FS}}({\theta }_{{\rm{E}}};{\rm{\Gamma }})$. This integrated detection efficiency is similar to the integrated detection efficiency, $\tilde{\epsilon }({t}_{{\rm{E}}};{\rm{\Gamma }})$, discussed in Sections 4.1.1 and 4.1.2, except that it has been integrated over t_E instead of θ_E for events with a significant FS signal, i.e., a measurement of ρ. The "relative detection efficiency" adopted for KMTNet by Gould et al. (2022) is actually a relative integrated detection efficiency in our nomenclature, which we think is more accurate. Their relative detection efficiency seems to be a ratio of the number of the events with the detection of the FS effect relative to the number of events with u₀ < ρ, while our ${\tilde{\epsilon }}_{\mathrm{FS}}({\theta }_{{\rm{E}}};{\rm{\Gamma }})$ in Koshimoto et al. (2023) is that relative to all events with u₀ ≤ 1. To compare these, we calculated the integrated detection efficiency with the FS effect relative to the events with u₀ < ρ, denoted as ${\tilde{\epsilon }}_{\mathrm{FS}}^{{\prime} }({\theta }_{{\rm{E}}};{{\rm{\Gamma }}}_{\mathrm{FS}})$ and shown in Figure 7. The integrated detection efficiency depends on the FFP MF, so we have used our best-fit MF to calculate these curves. This figure shows the MOA integrated detection efficiency as a function of θ_E for all sources (orange) and for giant sources with I_s,0 < 16 (blue). This is the same limit on I_s,0 as used by KMTNet (Gould et al. 2022), for their analysis of FSPL events. The green curve shows KMTNet's adopted relative integrated detection efficiency. Both the MOA I_s,0 < 16 curve and the KMTNet curves are normalized to match the MOA all-source integrated detection efficiency at ${\mathrm{log}}_{10}({\theta }_{{\rm{E}}})=-1.5$.

Zoom In Zoom Out Reset image size

The MOA sensitivity for giant sources is less than that for all sources at small θ_E because the large θ_* values for giant sources can significantly reduce the peak microlensing magnification. However, the sensitivity curve for KMTNet is very different from that of MOA, with a much sharper cutoff at small θ_E. This is partly because they directly cut off their integrated detection efficiency with a cut excluding events with θ_E < 3 μas. Gould et al. (2022) describe this cut by saying, "we complete this function linearly by imposing a threshold at θ_E = 3 μas, which is supported by the fact that all four FFPs are pressed up close to this limit." It is difficult to understand why they would need a cut like this given the sensitivity calculated for our analysis. Similarly, two of the four FFP events with FS effects found by OGLE (Mróz et al. 2018; Mróz et al. 2019, 2020, 2020) have θ_E < 3 μas (see Table 1) even though they have source stars with I_s,0 < 16. Perhaps the rationale for this cut that requires θ_E > 3 μas is to make their analysis consistent with their power-law prior assumption of 0.9 ≲ p ≲ 1.2. However, if this is the reason for this cut, it would raise the question as to why KMTNet has not been able to find events with θ_E < 3 μas in contrast to MOA and OGLE, which clearly have sensitivity well below this limit with bright sources with I_s,0 < 16. It would be helpful to see a full analysis for the KMTNet data set, including a complete detection efficiency analysis that includes both the t_E and θ_E dependence.

Note that our analysis does not use this integrated detection efficiency that depends only on θ_E. This integrated detection efficiency is included only for comparison with the Gould et al. (2022) analysis.