Constraints on Earth-mass primordial black holes from OGLE 5-year microlensing events

TL;DR

This paper uses OGLE 5-year microlensing data to constrain the abundance of Earth-mass primordial black holes, suggesting they could explain ultrashort-timescale events and providing new upper bounds on their dark matter contribution.

Contribution

It presents the first tight upper bounds on Earth-mass PBHs from OGLE data and proposes PBHs as an explanation for ultrashort-timescale microlensing events.

Findings

Upper bounds on PBH abundance in the Earth-Jupiter mass range.

PBHs can explain ultrashort-timescale microlensing events without free-floating planets.

Potential for future confirmation via M31 microlensing searches.

Abstract

We constrain the abundance of primordial black holes (PBH) using 2622 microlensing events obtained from 5-years observations of stars in the Galactic bulge by the Optical Gravitational Lensing Experiment (OGLE). The majority of microlensing events display a single or at least continuous population that has a peak around the light curve timescale $t_{\rm E}\simeq 20~{\rm days}$ and a wide distribution over the range $t_{\rm E}\simeq [1, 300]~{\rm days}$, while the data also indicates a second population of 6 ultrashort-timescale events in $t_{\rm E}\simeq [0.1,0.3]~{\rm days}$, which are advocated to be due to free-floating planets. We confirm that the main population of OGLE events can be well modeled by microlensing due to brown dwarfs, main sequence stars and stellar remnants (white dwarfs and neutron stars) in the standard Galactic bulge and disk models for their spatial and velocity distributions. Using the dark matter (DM) model for the Milky Way (MW) halo relative to the Galactic bulge/disk models, we obtain the tightest upper bound on the PBH abundance in the mass range $M_{\rm PBH}\simeq[10^{-6},10^{-3}]M_\odot$ (Earth-Jupiter mass range), if we employ null hypothesis that the OGLE data does not contain any PBH microlensing event. More interestingly, we also show that Earth-mass PBHs can well reproduce the 6 ultrashort-timescale events, without the need of free-floating planets, if the mass fraction of PBH to DM is at a per cent level, which is consistent with other constraints such as the microlensing search for Andromeda galaxy (M31) and the longer timescale OGLE events. Our result gives a hint of PBH existence, and can be confirmed or falsified by microlensing search for stars in M31, because M31 is towards the MW halo direction and should therefore contain a much less number of free-floating planets, even if exist, than the direction to the MW center.

Figures

Tables

Table 1. Table 1: Summary of the Galactic models for the mass and velocity distributions for stellar components and PBHs. α 𝛼 \alpha is the ratio of distances between lens and source, α ≡ d l / d s 𝛼 subscript 𝑑 l subscript 𝑑 s \alpha\equiv d_{\rm l}/d_{\rm s} . We employ the coordinate system as defined in Section  II.2.2 , and take the Galactic center as the coordinate origin. The optical depth ( τ 𝜏 \tau ) is calculated assuming an observation in the direction of ( b , l ) = ( − 2 ∘ ​ .389 , 1 ∘ ​ .0879 ) 𝑏 𝑙 superscript 2 .389 superscript 1 .0879 (b,l)=(-2^{\circ}\!\!.389,1^{\circ}\!\!.0879) which represents the OGLE Galctic bulge fields. For PBH case, we assume that PBHs constitute DM in the MW region by mass fraction, f PBH ≡ Ω PBH / Ω DM subscript 𝑓 PBH subscript Ω PBH subscript Ω DM f_{\rm PBH}\equiv\Omega_{\rm PBH}/\Omega_{\rm DM} , when computing the microlensing optical depth. We assume a Gaussian for the velocity profile, and the quantities, μ 𝜇 \mu and σ 𝜎 \sigma , denote the mean and dispersion for the Gaussian distribution. For PBH, we employ σ DM = 220 ​ km / s subscript 𝜎 DM 220 km s \sigma_{\rm DM}=220~{}{\rm km/s} for our fiducial model, which is taken from the rotational velocity of Galactic disk (see text).

lens	mass density profile: $\rho$ [$M_{\odot }{\mathrm{pc}}^{-3}$]	$\tau$ [${10}^{-6}$]	velocity profile: $(\mu ,\sigma )$ [km/s]
bulge		1.07	$f_y:\{-220\left(1-\alpha \right),\sqrt{1+{\alpha }^2}100\}$
	$3.53K_0\left(\frac{s}{667\mathrm{pc}}\right),\text{(}s\ge 938\text{pc)}$		$f_z:\{0,\sqrt{1+{\alpha }^2}100\}$
disk	$0.06\times \exp \left[-\left\{\frac{R-8000}{3500}+\frac{z}{325}\right\}\right]$	1.03	$f_y:\{220\alpha ,\sqrt{{\left(\kappa \delta +30\right)}^2+{\left(100\alpha \right)}^2}\}$
			$f_z:\{0,\sqrt{{\left(\lambda \delta +30\right)}^2+{\left(100\alpha \right)}^2}\}$
PBH	$4.88\times {10}^{-3}{f}_{\mathrm{PBH}}\frac{1}{\left(r/r_s\right){\left(1+r/r_s\right)}^2}$	0.18$f_{\mathrm{PBH}}$	$f_y:\{-220\left(1-\alpha \right),\sqrt{{\sigma }_{\mathrm{DM}}^2+{\left(\alpha 100\right)}^2}\}$
			$f_z:\{0,\sqrt{{\sigma }_{\mathrm{DM}}^2+{\left(\alpha 100\right)}^2}\}$

Table 2. Table 2: Summary of the mass spectrum for each of astrophysical objects: brown dwarfs, main-sequence stars, and stellar remnants (white dwarfs, neutron stars, and astrophysical black holes) in the standard Galactic bulge and disk models. We assume the Kroupa initial mass function as shown in Fig.  1 , and then assume that each massive star with initial masses M ≥ 1 ​ M ⊙ 𝑀 1 subscript 𝑀 direct-product M\geq 1M_{\odot} , as denoted in the column “initial mass range”, evolved into each stellar remnant. For white dwarf, we assume the relation between initial and end masses as given by M WD = 0.339 + 0.129 ​ M init subscript 𝑀 WD 0.339 0.129 subscript 𝑀 init M_{\rm WD}=0.339+0.129M_{\rm init} . The column “ d ​ n / d ​ M d 𝑛 d 𝑀 \mathrm{d}n/\mathrm{d}M ” denotes parameters of the mass spectrum for each object population, while we assume a Gaussian distribution with the mean and width values for the mass spectrum of neutron stars and black holes. The last column “ N 𝑁 N ” gives the number of each object population relative to that of main-sequence stars used in the calculation of microlensing event rate.

object	parameters in $dn/dM$	mass range $\left[M_{\odot }\right]$	initial mass range $\left[M_{\odot }\right]$	$N$
brown dwarf (BD)	Power-law ($M^{-0.8}$)	[0.01, 0.08]	$0.01\le M\le 0.08$	0.18
main-sequence star (MS)	Power-law ($M^{-2}$)	[0.5, 1.0]	$0.5\le M\le 1.0$	1
	Power-law ($M^{-1.3}$)	[0.08 0.5]	$0.08\le M\le 0.5$
white dwarf (WD)	Power-law (initially $M^{-2}$)	[0.34, 2.0]	$1.0\le M\le 8.0$	0.15
neutron star (NS)	Gaussian ($M_r=1.33,{\sigma }_r=0.12$)	[0.73, 1.93]	$8.00\le M\le 20.0$	0.013
black hole (BH)	Gaussian ($M_r=7.8,{\sigma }_r=1.2$)	[1.8, 13.8]	$20.0\le M\le 100.0$	0.0068

Equations

$$A(t)=\frac{u^2+2}{u\sqrt{u^2+4}},$$

$$u(t)=\sqrt{{\beta }^2+\frac{(t-t_0{)}^2}{t_{\mathrm{E}}^2}}.$$

$$t_{\mathrm{E}}\equiv \frac{R_{\mathrm{E}}}{v}=\frac{\sqrt{4GM{d}_{\mathrm{l}}{d}_{\mathrm{ls}}/d_{\mathrm{s}}}}{cv},$$

$$t_{\mathrm{E}}\simeq 44\mathrm{days}{\left(\frac{M}{M_{\odot }}\right)}^{1/2}{\left(\frac{d_{\mathrm{l}}{d}_{\mathrm{ls}}/d_{\mathrm{s}}}{4\mathrm{kpc}}\right)}^{1/2}{\left(\frac{v}{220km/s}\right)}^{-1}.$$

$$\tau \equiv {\tau }_{\mathrm{b}}+{\tau }_{\mathrm{d}}+{\tau }_{\mathrm{PBH}}.$$

$$\frac{\mathrm{d}{\Gamma }_a}{\mathrm{d}{t}_{\mathrm{E}}}\equiv \frac{{\mathrm{d}}^2{\tau }_a}{\mathrm{d}{t}_{\mathrm{obs}}\mathrm{d}{t}_{\mathrm{E}}},$$

$${\tau }_{\mathrm{b}}$$

$$=\frac{4\pi G}{c^2{N}_{\mathrm{s}}}{\int }_{d_{\mathrm{s},\min }}^{d_{\mathrm{s},\max }}\mathrm{d}{d}_{\mathrm{s}}{n}_{\mathrm{s}}(d_{\mathrm{s}})\sum _i{\int }_{d_{\mathrm{s},\min }}^{d_{\mathrm{s}}}\mathrm{d}{d}_{\mathrm{l}}{\rho }_{\mathrm{b},i}(d_{\mathrm{l}})D,$$

$$N_{\mathrm{s}}\equiv {\int }_{d_{\mathrm{s},\min }}^{d_{\mathrm{s},\max }}\mathrm{d}{d}_{\mathrm{s}}{n}_{\mathrm{s}}(d_{\mathrm{s}}).$$

$$\mathrm{d}{\Gamma }_{\mathrm{b}}$$

$$f_{\mathrm{b}}({\mathbf{v}}_{\perp },v_{\parallel })=f_{\mathrm{b}}({\mathbf{v}}_{\perp })f_{\mathrm{b}}(v_{\parallel }).$$

$$\frac{\mathrm{d}{\Gamma }_{\mathrm{b}}}{\mathrm{d}{t}_{\mathrm{E}}}$$

$$\times {\int }_0^{\infty }\mathrm{d}{v}_{\perp }{\int }_{-\pi /2}^{\pi /2}\mathrm{d}\theta v_{\perp }^2\cos \theta f_{\mathrm{b},i}(v_{\perp },\theta )$$

$$\times {\delta }_D\left(t_{\mathrm{E}}-\frac{2R_{\mathrm{E}}\cos \theta }{v_{\perp }}\right).$$

$${\delta }_D\left(t_{\mathrm{E}}-\frac{2R_{\mathrm{E}}\cos \theta }{v_{\perp }}\right)={\delta }_D\left(v_{\perp }-\frac{2R_{\mathrm{E}}\cos \theta }{t_{\mathrm{E}}}\right)\frac{v_{\perp }^2}{2R_{\mathrm{E}}\cos \theta },$$

$$\frac{\mathrm{d}{\Gamma }_{\mathrm{b}}}{\mathrm{d}{t}_{\mathrm{E}}}$$

$$\times \sum _i{\int }_{d_{\mathrm{s},\min }}^{d_{\mathrm{s}}}\mathrm{d}{d}_{\mathrm{l}}\frac{{\rho }_{\mathrm{b},i}(d_{\mathrm{l}})}{M_i}{\int }_{-\pi /2}^{\pi /2}\mathrm{d}\theta v_{\perp }^4{f}_{\mathrm{b},i}(v_{\perp },\theta ),$$

$${\tau }_{\mathrm{d}}=\frac{4\pi G}{c^2}{\int }_0^{{\bar{d}}_{\mathrm{s}}}\mathrm{d}{d}_{\mathrm{l}}\sum _i{\rho }_{\mathrm{d},i}(d_{\mathrm{l}})D,$$

$$\frac{\mathrm{d}{\Gamma }_{\mathrm{d}}}{\mathrm{d}{t}_{\mathrm{E}}}=\pi \sum _i{\int }_0^{{\bar{d}}_{\mathrm{s}}}\mathrm{d}{d}_{\mathrm{l}}\frac{{\rho }_{\mathrm{d},i}(d_{\mathrm{l}})}{M_i}{\int }_{-\pi /2}^{\pi /2}\mathrm{d}\theta v_{\perp }^4{f}_{\mathrm{d},i}(v_{\perp },\theta ),$$

$${\tau }_{\mathrm{PBH}}=\frac{4\pi G}{c^2}{\int }_0^{{\bar{d}}_{\mathrm{s}}}\mathrm{d}{d}_{\mathrm{l}}{\rho }_{\mathrm{DM}}(d_{\mathrm{l}})D,$$

$$\frac{\mathrm{d}{\Gamma }_{\mathrm{PBH}}}{\mathrm{d}{t}_{\mathrm{E}}}=\pi {\int }_0^{{\bar{d}}_{\mathrm{s}}}\mathrm{d}{d}_{\mathrm{l}}\frac{{\rho }_{\mathrm{DM}}(d_{\mathrm{l}})}{M_{\mathrm{PBH}}}{\int }_{-\pi /2}^{\pi /2}\mathrm{d}\theta v_{\perp }^4{f}_{\mathrm{DM}}(v_{\perp },\theta ),$$

$${\rho }_{\mathrm{b}}(x,y,z)$$

$$=\{\begin{array}{ll}1.04\times 10^6{\left(\frac{s}{0.482\mathrm{pc}}\right)}^{-1.85}{M}_{\odot }\mathrm{p}{\mathrm{c}}^{-3},& \text{(s<938pc)},\\ 3.53K_0\left(\frac{s}{667\mathrm{pc}}\right)M_{\odot }\mathrm{p}{\mathrm{c}}^{-3},& \text{(s≥938pc)},\end{array}$$

$$r(d)=\sqrt{R_{\oplus }^2-2R_{\oplus }d\cos l\cos b+r^2},$$

$${\rho }_{\mathrm{d}}(R,z)=0.06\times \exp \left[-\left\{\frac{R-8000}{3500}+\frac{z}{325}\right\}\right]M_{\odot }{\mathrm{pc}}^{-3}.$$

$${\rho }_{\mathrm{NFW}}(r)=\frac{{\rho }_c}{(r/r_s)(1+r/r_s{)}^2},$$

$${\mathbf{v}}_{\perp }$$

$$={\mathbf{v}}_{\mathrm{l}}-\left[\alpha {\mathbf{v}}_{\mathrm{s}}+(1-\alpha ){\mathbf{v}}_{\mathrm{o}}\right],$$

$${\bar{v}}_{\mathrm{b}y}\equiv ⟨v_{\mathrm{l}y}-\left[\alpha v_{\mathrm{s}y}-(1-\alpha )v_{\mathrm{o}y}\right]⟩=-220(1-\alpha )\text{mbox}km/s,$$

$${\bar{v}}_{\mathrm{b}z}=0$$

$${\sigma }_{\mathrm{b}y}^2$$

$$=⟨v_{\mathrm{l}y}^2⟩+{\alpha }^2⟨v_{\mathrm{s}y}^2⟩$$

$$=(1+{\alpha }^2)(100\text{mbox}km/s{)}^2,$$

$${\sigma }_{\mathrm{b}z}^2=(1+{\alpha }^2)(100km/s{)}^2.$$

$$f_{\mathrm{b}}(\mathbf{v})=f_{\mathrm{b}}(v_y)f_{\mathrm{b}}(v_z),$$

Full text

dark matter-black holes-Galactic halos-Spiral arms and bars; galactic disks

I Introduction

Dark matter (DM) is one of the most essential ingredients in the standard model of cosmic hierarchical structure formation scenario (e.g., Springel et al., 2005). However, the nature of DM is not yet known and one of the most important, unresolved problems in physics and astronomy. Unknown stable elementary particle(s) beyond the Standard Model of Particle Physics, so-called Weakly Interacting Massive Particle(s) (WIMP), have been thought of as a viable candidate of DM, but has yet to be detected either in direct experiments, collider experiments or indirect searches (e.g., Jungman et al., 1996; Arina, 2018; Hooper, 2018, for a review). Primordial black holes (PBH), first proposed in pioneer works Hawking (1971); Carr and Hawking (1974); Carr (1975), are alternative viable candidate of DM (also see Carr et al., 2016, for a review). PBHs can be formed by gravitational collapse of Hubble patch in the early universe if the patch has a large primordial overdensity of $\delta\sim O(0.1)$. Various works have proposed a mechanism to generate such a large overdensity such as an inflation model (e.g., Inomata et al., 2017a; Cotner and Kusenko, 2017; Cotner et al., 2018), producing PBHs with a target mass scale and a target abundance. Furthermore, there is a renewed interest in PBH scenario because of recent claims (Sasaki et al., 2016; Bird et al., 2016; Sasaki et al., 2018; Ali-Haïmoud et al., 2017; Nakamura et al., 2016) (also see Ioka et al., 1998) that PBHs can be progenitors of binary black holes whose gravitational wave events have been detected by the LIGO/Virgo experiments (Abbott et al., 2016; Abbott and the others, 2017).

There are various attempts at constraining PBHs over almost twenty orders of magnitudes in their mass scales; gamma-ray background from PBH evaporation (Carr et al., 2010), femtolensing of gamma-ray bursts (Gould, 1992) (although Ref. Katz et al. (2018) recently pointed out that a finite-source size effect of the gamma-ray burst progenitor significantly relaxes or even removes the constraint), supernovae of white dwarfs triggered by PBH (Graham et al., 2015), PBH capture by neutron stars (Capela et al., 2013; Fuller et al., 2017), microlensing constraints (Alcock et al., 2000; Tisserand et al., 2007; Griest et al., 2014; Niikura et al., 2017), caustics-network lensing in the galaxy cluster region (Kelly et al., 2018; Diego et al., 2018; Oguri et al., 2018), $X$-ray background from gas accretion on PBH (Inoue and Kusenko, 2017), and the effect of PBH gas accretion on the cosmic microwave background optical depth (Ricotti et al., 2008) (see Refs. (Ali-Haïmoud and Kamionkowski, 2017; Aloni et al., 2017) for the revisited calculations), the effect of PBH on pulsar timing array observation Schutz and Liu (2017) and the effect on type-Ia supernova observation (Zumalacárregui and Seljak, 2018). Except for a mass window of $M_{\rm PBH}\simeq[10^{-16},10^{-11}]M_{\odot}$, these constraints rule out a scenario that PBHs constitute the dominant fraction of DM if PBH has a narrow mass spectrum.

Microlensing is the most robust, powerful tool among various methods to probe a compact, macroscopic DM in the Milky Way (MW) halo region (Paczynski, 1986; Griest et al., 1991), because lensing is a gravitational effect and can directly probe mass (gravity strength) of a lensing object irrespective of whether a lensing object is visible or not. The Optical Gravitational Lensing Experiment (OGLE 111http://ogle.astrouw.edu.pl) collaboration (Udalski et al., 1994, 2015) has been making invaluable long-term efforts, more than a decade, to make monitoring observations of million stars in the Galactic bulge fields. The OGLE team has been finding more than two thousands of microlensing events and obtained various constraints on exoplanetary systems, brown dwarfs, low-mass stars as well as presented even an indication of free-floating planets in inter-stellar space (Mróz et al., 2017) (also see Sumi et al., 2003, 2011, for the similar constraints from the MOA microlensing experiments).

In this paper we use the 5-years OGLE data containing 2622 microlensing events in Mróz et al. (2017) to constrain the PBH abundance. Interestingly the OGLE data indicates 6 ultrashort-timescale microlensing events that have their light curve timescales of $t_{\rm E}\simeq[0.1,0.3]~{}{\rm days}$ (also see Mroz et al., 2018, for the new candidates), which is a distinct population from the majority of OGLE events. The ultrashort-timescale events indicate Earth-mass “unbounded” (wide-orbit or free-floating) planets (Sumi et al., 2011; Mróz et al., 2017). However, the origin of such free-floating planets is poorly understood because it involves complicated physics of star formation, planetary system formation and interaction of planetary system with other stars/planets. Hence we pay a particular attention to a possibility of whether PBHs can give an alternative explanation of the ultrashort-timescale events. For this purpose, we first study the standard Galactic bulge and disk models to estimate event rates of microlensing events due to astrophysical objects including brown dwarfs, main sequence stars and stellar remnants (white dwarfs, neutron stars and astrophysical black holes) following the pioneer work in Han and Gould (1995) (also see Han and Gould, 1996). After comparing the model predictions of astrophysical objects with the OGLE data including the calibration factor that takes into account observational effects, we use the OGLE data to constrain the PBH abundance using the standard MW halo model for the spatial and velocity distributions of DM (therefore PBHs). In doing this we employ two working hypotheses. First, we employ “null hypothesis” that there is no PBH microlensing in the OGLE data and then derive an upper bound on the PBH abundance. Second, we employ the assumption that the 6 ultrashort-timescale events are due to PBHs and derive an allowed region of PBHs in the mass and abundance parameter space. To obtain the results, we properly use the likelihood function of OGLE events assuming the Poisson uncertainty in the counts of microlensing events.

The structure of this paper is as follows. In Section II, we briefly review basics of microlensing and give equations relevant for the event rate calculations. In Section III we review the standard models for the Galactic disk and bulge describing the spatial and velocity distributions for brown dwarfs, stars and stellar remnants as the constitutions, and also describe the MW halo model for the distributions of DM, i.e. PBH in our study. In Section IV, we give the main results of this paper, after reviewing the OGLE data; the upper bound on the PBH abundance and a possible implication of Earth-mass PBHs. We will then give conclusion and discussion in Section V.

II Microlensing for bulge stars

II.1 Microlensing basics

When a source star and a lensing PBH are almost perfectly aligned along the line-of-sight direction of an observer, the star is multiply imaged due to strong lensing (Paczynski, 1986) (also see Meylan et al., 2006; Dodelson, 2017, for a review). In case these multiple images are unresolved, the flux from the star appears magnified. The light curve of such a microlensing magnification event is given by

[TABLE]

where $u(t)$ is the separation between the source star and the lens at an observation epoch $t$, and is given by

[TABLE]

Here $\beta$ is the dimensional-less impact parameter, and $t_{0}$ is the time where the source and the lens become closest in separation. Throughout this paper we use a crossing time of the Einstein radius, denoted as $t_{{\rm E}}$, to characterize a timescale of microlensing light curve:

[TABLE]

where $R_{\rm E}$ is the Einstein radius, $M$ is the mass of lensing object (we assume a point mass for lens throughout this paper), $v$ is the (total) relative velocity on the two-dimensional plane perpendicular to the line-of-sight direction (see below), $d_{\rm l}$ and $d_{\rm s}$ are distances to the lens and source, respectively, and $d_{\rm ls}$ is the distance between lens and source ($d_{\rm ls}\equiv d_{\rm s}-d_{\rm l}$).

If we plug typical values of the physical quantities into Eq. (3), we can find a typical timescale of the microlensing light curve as

[TABLE]

This equation shows that a lighter-mass lens causes a shorter timescale light curve. As we will show below, the 5-year OGLE data gives a sufficient sampling of microlensing light curves over the range of timescales, $[10^{-1},300]~{}$days, roughly corresponding to lenses of mass scales, $[10^{-6},10]~{}M_{\odot}$. In reality, since lensing objects have a velocity distribution, it causes a distribution of the light curve timescales even if lensing objects are in a narrow mass bin, which we need to properly take into account.

II.2 Microlensing optical depth and event rate for a star in the Galactic bulge

II.2.1 Definition of microlensing optical depth and event rate

Here we define the optical depth and event rate of microlensing for a single star in the Galactic bulge region. The optical depth is defined as the probability for a source star to be inside the Einstein radius of a foreground lensing object on the sky at a certain moment. This corresponds to the probability for the lensing magnification to be greater than $A\geq 1.34$. The total optical depth due to lensing objects in the bulge and disk regions as well as due to PBHs in the MW halo region is formally expressed as

[TABLE]

Hereafter we employ abbreviations: “b” for “bulge” and “d” for “disk”, respectively, and we ignore a multiple lensing case for a single star (this is a good approximation given the low optical depth as we show below). For a lensing object in the bulge and disk regions, we consider brown dwarfs and stellar components, where the latter includes main sequence stars and stellar remnants (white dwarfs, neutron stars and astrophysical black holes), as we will explain in detail later.

The differential event rate of microlensing is defined as the frequency of microlensing events of a given light curve timescale (denoted as $t_{\rm E}$) for a single source star per unit observational time ($t_{\rm obs}$):

[TABLE]

where the subscript $a={\rm bulge}$, disk or PBH, respectively.

II.2.2 Coordinate system

It would be useful to explicitly define the coordinate system we employ in the following calculations. For the rectangular coordinate system, denoted as $(x,y,z)$, we choose the Galactic center as the coordinate origin. Without loss of generality, we can take the $x$-direction to be along the direction connecting the Galactic center and the Earth position (an observer’s position). We assume that the Earth is located at the position, $(x,y,z)_{\oplus}=(8~{}{\rm kpc},0,0)$, i.e. 8 kpc in distance from the Galactic center. Furthermore, we take the $y$-direction to be along the Earth’s rotation direction in the Galactic disk plane, and the $z$-direction to be in the direction perpendicular to the disk plane.

In this paper we consider the microlensing datasets obtained from the 5-years OGLE survey (Udalski et al., 2015)222The OGLE-IV fields can be found from http://ogle.astrouw.edu.pl/sky/ogle4-BLG/.. In the Galactic coordinates, the OGLE fields are located in the range of -15^{\circ}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}b\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}15^{\circ} and -20^{\circ}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}l\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}20^{\circ}. Throughout this paper we simply assume that the OGLE field is in the direction to the field BLG505 with $(b,l)=(-2^{\circ}\!\!.389,1^{\circ}\!\!.0879)$, which has the largest number of background stars among the OGLE fields. We believe that this approximation is valid because our results are based on relative contributions of microlensing due to stellar components in the bulge and disk regions compared to microlensing due to PBHs in the MW halo region.

II.2.3 Bulge lens

First we consider a microlensing that both lens and source are in the Galactic bulge region. The average optical depth of microlensing due to the $i$-th stellar component for a single source star is given by

[TABLE]

where the integration is along the line-of-sight direction (see below), $D\equiv d_{\rm l}d_{\rm ls}/d_{\rm s}$, and the index $i$ stands for the $i$-th stellar component as a lensing object for which we will consider brown dwarfs, main sequence stars, white dwarfs, neutron stars and astrophysical black holes (see below). $N_{\rm s}$ is the surface number density of source stars defined by a line-of-sight integration of the three-dimensional number density distribution of source stars, $n_{{\rm s}}$, as

[TABLE]

The function $\rho_{{\rm b},i}(d_{\rm l})$ is the mass density profile for the $i$-th stellar component. $d_{\rm s,min}$ and $d_{\rm s,max}$ are the maximum and minimum distances to the boundary of the bulge region from an observer’s position. Throughout this paper we employ $d_{\rm s,min}=4~{}{\rm kpc}$ and $d_{\rm s,max}=12~{}{\rm kpc}$; that is, we assume that the bulge has a size of $4~{}{\rm kpc}$ around the center ($d_{\rm s}=8~{}{\rm kpc}$) in depth from an observer. The integration over $d_{{\rm s}}$ or $d_{{\rm l}}$ in Eq. (7) is along the line-of-sight direction of an observer towards the source star in the direction ($y,z$). As long as the lens distribution $\rho_{{\rm b},i}$ and the source star distribution are given, the line-of-sight integration is straightforward to perform as we will show later.

Now we consider the event rate of microlensing. To do this, we start from the geometry and variables defined in Fig. 4 of Griest et al. (1991) (also see Figure 7 of Niikura et al., 2017), which defines the differential event rate of a lensing object entering a volume element along the line-of-sight where the lens causes a microlensing with magnification above a certain threshold value:

[TABLE]

where $f_{{\rm b},i}({\bf v}_{\perp},v_{\parallel})$ is the velocity distribution of the $i$-th stellar component, defined so as to satisfy the normalization condition $\int\!\!\mathrm{d}^{2}{\bf v}{\perp}\int\!\!\mathrm{d}v_{\parallel}~{}f({\bf v}_{\perp},v_{\parallel})=1$; ${\bf v}{\perp}$ is the perpendicular components of relative velocity between an observer, lens and source star (see below), defined as ${\bf v}_{\perp}=v_{\perp}(\cos\theta,\sin\theta)$; $\alpha$ is the azimuthal angle in the two-dimensional ($y,z$)-plane, defined as $(y,z)=\sqrt{y^{2}+z^{2}}(\cos\alpha,\sin\alpha)$; $\theta$ is the angle between the line connecting the source and the lens center and the direction of the transverse velocity ${\bf v}_{\perp}$. In this paper we define the microlensing “event” if the lensing magnification is greater than a threshold magnification, $A>A(R_{\rm E})=1.34$, which is satisfied if the separation between lens and source is closer than the threshold separation, $b\leq R_{{\rm E}}$. The parameters vary in the range of $\theta\in[-\pi/2,\pi/2]$, $\alpha\in[0,2\pi]$, and $v_{\perp}\in[0,\infty)$.

We assume that the velocity distribution can be simplified as

[TABLE]

As we discussed, for a characteristic timescale of microlensing light curve, we employ a crossing time scale of the Einstein ring, defined as $t_{\rm E}=2R_{\rm E}\cos\theta/v_{\perp}$. This simplification is not critical for the following discussion because we study the PBH microlensing contribution relative to those due to the stellar components in the disk and bulge regions. In this case, the microlensing event rate due to the $i$-th stellar components is given as

[TABLE]

Using the Dirac delta function identity

[TABLE]

the above equation is simplified as

[TABLE]

where $v_{\perp}=2R_{{\rm E}}\cos\theta/t_{{\rm E}}$. With this condition, the tangential velocity $v_{\perp}$ depends on integration variables, $d_{\rm s},d_{\rm l}$, and $\theta$ via $R_{\rm E}=R_{\rm E}(d_{\rm l},d_{\rm s})$.

II.2.4 Disk lens

Next we consider an event rate for microlesning due to stellar components in the disk region for a single source star in the bulge region. The calculation is very similar to the case for bulge lens in the preceding section. In this case we employ a single source plane approximation for simplicity; that is, we assume that all source stars are at distance of $8~{}$kpc, the Galactic center. Under this assumption, the optical depth is given by

[TABLE]

where $\rho_{{\rm d},i}(d_{\rm l})$ is the mass density distribution of the $i$-th stellar component, $D=d_{{\rm l}}d_{{\rm ls}}/\bar{d}_{{\rm s}}$, $\bar{d}_{{\rm s}}$ is the mean distance to source stars, i.e. $\bar{d}_{\rm s}=8~{}{\rm kpc}$, and $d_{{\rm ls}}=\bar{d}_{{\rm s}}-d_{{\rm l}}$.

Similarly, the event rate is

[TABLE]

where $v_{\perp}=2R_{{\rm E}}\cos\theta/t_{{\rm E}}$, $f_{{\rm d},i}({\bf v}_{\perp})=f_{{{\rm d}},i}(v_{\perp},\theta)$ is the velocity distribution for velocity components perpendicular to the ling-of-sight direction for the $i$-th stellar component in the disk region.

II.2.5 PBH lens

Now we consider a scenario that PBHs constitutes some mass fraction of DM in the MW halo region. We call “PBHs in the halo region” because PBHs are distributed from the Galactic center through the outer halo region due to the large velocity dispersion. When a lensing PBH happens to pass across a source star in the bulge on the sky, it could cause microlensing effect on the source star. Throughout this paper we consider a monochromatic mass distribution for PBHs. Similarly to the disk microlensing, the optical depth of microlensing due to PBHs is

[TABLE]

where $\rho_{{\rm DM}}(d_{\rm l})$ is the dark matter distribution. If PBHs consist only some partial mass fraction of DM in the MW region, denoted as $f_{\rm PBH}\equiv\Omega_{\rm PBH}/\Omega_{\rm DM}$, we replace $\rho_{\rm DM}$ in the above and following equations with $f_{\rm PBH}\rho_{\rm DM}$.

Similarly, the event rate of microlensing due to PBHs for a single source star in the bulge is

[TABLE]

where $f_{{\rm DM}}({\bf v}_{\perp})$ is the velocity distribution of PBHs.

III Models of Galactic disk and bulge and Milky Way dark matter

As we described, once we give the density and velocity distributions for stellar components in the MW bulge and disk regions as well as those for PBHs (equivalently DM) in the halo region, we can compute the event rates of microlensing for a star in the bulge region. In this subsection, we briefly review the standard model for the MW bulge and disk following Han and Gould (1995). Then we describe our model for the density and velocity distributions for PBHs in the MW region.

III.1 The mass density distribution

For the mass density distribution of stellar population in the bulge region, we adopt the model in Kent (1992) that describes the following bar-structured model:

[TABLE]

where $K_{0}(x)$ is the modified Bessel function, $s$ is the radius from the Galactic center in the elliptical coordinates, defined as $s^{4}\equiv R^{4}+(z/0.61)^{4}$ with $R\equiv(x^{2}+y^{2})^{1/2}$, and all the coordinate components ($x,y,z,s,R$) are in units of pc. As defined in Section II.2.2, the coordinate origin is the Galactic center. Note that the above profile is continuous at $s=938~{}{\rm pc}$, and we consider the above profile is for the total contribution of visible objects, i.e main sequence stars, as we will describe below. Using the Galactic celestial coordinate variables $(l,b)$, a star at the distance $d$ from an observer (the Earth’s position) is at the distance from the Galactic center, $r$, given as

[TABLE]

where $r=\sqrt{x^{2}+y^{2}+z^{2}}$, $x=d\cos b\cos l$, $y=d\cos b\sin l$, and $z=d\sin b$. This variable transformation between $d$ ($d_{\rm l}$ or $d_{\rm s}$) and ${x,y,z}$ enters into the above equations such as Eq. (13).

For the mass density distribution in the disk region, we employ the model in Bahcall (1986):

[TABLE]

Note that, as we defined, $R(=\sqrt{x^{2}+y^{2}})$ denotes the radial distance in the cylindrical coordinates and $z$ is in the direction perpendicular to the Galactic disk (variables are in units of pc). This model assumes that the disk has an exponential distribution with vertical and radial scale lengths of 325 pc and 3500 pc, respectively. Although the mass-to-light ratio of disk stellar population is not well understood, we normalize the above density profile to $\rho_{\rm d0}=0.06~{}M_{\odot}\mathrm{pc}^{-3}$ at the solar neighborhood ($R=8000~{}{\rm pc}$).

For the spatial distribution of DM (therefore PBHs) between the Galactic center and an observer (the Earth), we assume the Navarro-Frenk-White (NFW) model (Navarro et al., 1997):

[TABLE]

where $r_{s}$ is the scale radius and $\rho_{c}$ is the central density parameter. For this model we assume spherical symmetry for the DM distribution for simplicity. In this paper we adopt the halo model in Klypin et al. (2002): $M_{\rm vir}=10^{12}M_{\odot}$, $\rho_{c}=4.88\times 10^{6}~{}{M_{\odot}/{\rm kpc}^{3}}$, and $r_{s}=21.5~{}{\rm kpc}$, taken from Table 2 in the paper. The DM profile with these parameters has been shown to fairly well reproduce the observed rotation curve in the MW. However, there might still be a residual uncertainty in the total mass (mostly DM) of MW within a factor of 2 (Callingham et al., 2018).

Table 1 summarizes models of the mass density profiles for stellar components in the bulge and disk regions and for dark matter in the MW halo region. The table also gives the optical depth for a single source star for lenses in the bulge or disk region and for PBHs, respectively. Note that the optical depth does not depend on a lens mass as indicated from Eq. (7). The table shows that the optical depth due to PBHs is smaller than that of astrophysical objects in the disk or bulge region, by a factor of 5, reflecting that astrophysical objects are more centrally concentrated due to the dissipation processes.

III.2 The velocity distribution

A timescale of the microlensing light curve (see Eqs. 13, 15, and 17) is determined by a transverse component of the relative velocity for source-lens-observer system on the sky (Griest et al., 1991; Han and Gould, 1995):

[TABLE]

where ${\bf v}_{{\rm l}}$, ${\bf v}_{\rm s}$ and ${\bf v}_{\rm o}$ are the transverse velocities for lens, source star and an observer, respectively, and we have introduced the notation $\alpha\equiv d_{{\rm l}}/d_{{\rm s}}$. As we described in Section II.2.2, the $x$-direction is along the direction from the observer to the Galactic center (i.e. a source star), which is equivalent to the line-of-sight direction, the $y$-direction is along the direction of disk rotation, and the $z$-direction is perpendicular to the line-of-sight direction. Hence we need to model the mean and distribution of the transverse velocity components, ${\bf v}_{\perp}=(v_{y},v_{z})$. Hereafter we often omit the subscript “$\perp$” in ${\bf v}_{\perp}$ for notational simplicity.

III.2.1 Bulge lens

First we consider the velocity distribution for the bulge microlensing where both lens and source star are in the bulge region. For the velocity distribution, we assume that the stellar components are supported by an isotropic velocity dispersion, and do not have any rotational velocity component. Under these assumptions, the mean of the transverse velocities is

[TABLE]

where we have assumed that an observer is in the rest frame of the rigid body rotation of Galactic disk, has the rotational velocity of 220 km/s with respect to the Galactic center, and has no mean velocity in the disk height direction.

The velocity dispersion for the $y$-component of relative velocity can be computed as

[TABLE]

Here we assumed that the velocity dispersion per component $\sigma_{y}=100$ km/s, and assumed that the source and lens have independent random motions; $\left\langle v_{{\rm s}y}v_{{\rm l}y}\right\rangle=0$. The velocity dispersion for the velocity $z$-component is

[TABLE]

Following Han and Gould (1995), we assume that the velocity distribution is given by a Gaussian and that the velocity distribution for the bulge microlensing is given by

[TABLE]

where

[TABLE]

III.2.2 Disk lens

Next we consider the velocity distribution for stellar components in the disk region. We assume that the stellar components have a rigid rotation on average:

[TABLE]

For the velocity dispersion, we assume the linear disk velocity dispersion model in Table 1 of Han and Gould (1995):

[TABLE]

with

[TABLE]

where $\kappa$ and $\lambda$ are the velocity dispersion gradient coefficients. Hence the velocity distribution functions are given by the similar equations to Eq. (27).

III.2.3 PBH lens

Now we consider a microlensing due to PBHs, acting as DM, in the MW halo region. PBHs are tracers of the MW halo that is much more extended than the bulge size ($\sim 200$ kpc vs. a few kpc in radius). The large extent of DM halo reflects the fact that PBHs have a larger velocity dispersion than that of bulge stars (100 km/s). First we assume that PBHs have isotropic velocity distribution with respect to the halo center for which we assume the Galactic center. Hence the mean relative velocity for a PBH lens is

[TABLE]

For a PBH causing microlesning effect on a bulge star, it should be located somewhere between the Galactic center and the Earth, which is a very inner region compared to the halo size. Such a PBH (more generally dark matter) tends to have a large velocity when passing through the central region of DM halo; DM tends to have a larger velocity at the closest point to the halo center (i.e. the Galactic center), while a bounded DM should stop at an apocenter point of its orbit, which tends to be around the outer boundary of the halo. Hence we assume that PBHs causing the microlensing have a large velocity dispersion whose amplitude is similar to the rotation velocity. To keep generality of our discussion, we introduce a parameter to model the velocity dispersion of DM per one direction, $\sigma_{\rm DM}$:

[TABLE]

where we have again assumed the isotropic velocity dispersion. For our fiducial model, we assume $\sigma_{\rm DM}=220~{}{\rm km/s}$. We checked that a change in $\sigma_{\rm DM}$, say by $\pm 10\%$, gives only a small change in the following PBH constraints. Such a large velocity dispersion in the central region within the halo is supported by $N$-body simulation studies that simulate MW-scale halos, for example, Figure 2 of Ref. (Vogelsberger et al., 2009).

Table 1 gives a summary of the velocity distributions for stellar components in the bulge and disk regions and for PBHs (DM) in the halo region, respectively.

III.3 Mass spectrum of astrophysical lensing objects

As implied by Eq. (4), a timescale of microlensing light curve varies with a lens mass. To make a quantitative modeling of microlensing events as a function of the light curve timescale, we need to take into account the mass distribution of stellar lenses in the bulge and disk regions. Following Mróz et al. (2017) (also see Gould, 2000; Sumi et al., 2011), we consider brown dwarfs (BD), main sequence stars (MS), white dwarfs (WD), neutron stars (NS), and astrophysical black holes (BH) as the constituents in the bulge and disk regions.

To do this, we first assume the Kroupa-like broken power-law initial mass function (IMF) for the stellar components (Kroupa, 2001):

[TABLE]

for BD, low-mass MS stars ($0.08<M/M_{\odot}<0.5$), and high-mass MS stars ($M\geq 0.5M_{\odot}$), respectively; $A_{\rm BD}$ and $A_{\rm MS}$ are normalization parameters for which we will discuss below, and $\alpha_{\rm BD},\alpha_{\rm MS1}$ and $\alpha_{\rm MS2}$ are the power-law index parameters for these components, respectively. Following Mróz et al. (2017) we assume $\alpha_{\rm BD}=0.8$, $\alpha_{\rm MS1}=1.3$ and $\alpha_{\rm MS2}=2$, respectively, where the slope for low-mass stars of \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}1~{}M_{\odot} is taken from the study of the Galactic bulge IMF in Ref. (Zoccali et al., 2000). Throughout this paper we assume the same population composition of stellar components in the disk and bulge regions. The formation of BDs is still poorly understood (see Burrows et al., 2001, for a review). Some fraction of BDs can be found in the planetary disk around a primary main-sequence star. Other population of BDs can form at the center of protoplanetary disk as a primary gravitating object of the system. Moreover, there might be some population of BDs ejected from the host system due to three-body scattering, which would be observed as “free-floating planets” in the interstellar space. A wide-orbit BD or a free-floating BD causes microlensing event characterized by the BD mass. If a lens system has both primary star and BD in the close orbit, the microlensing event is characterized by the total mass (mainly the host star). Thus, as discussed in Mróz et al. (2017) (Sumi et al., 2011), there is still a lot of discussion for the origin of microlensing events in a short timescale corresponding to BD masses or even shorter (smaller-mass) events. Hence, the amplitude of BD mass function is uncertain, and needs to be further study. In the following results, we will treat the BD normalization parameter $A_{\rm BD}$ as a free parameter, and determine it so that the model prediction matches the OGLE data in the corresponding short timescales.

Massive stars with masses $M_{\rm init}\geq 1~{}M_{\odot}$ have a rapid time evolution during the age of MW, and evolved into stellar remnants. Following Mróz et al. (2017), we assume that all stars with initial masses $1\leq M/M_{\odot}\leq 8$ evolved into WDs following the empirical initial-final mass relation, $M_{\rm WD}=0.339+0.129M_{\rm init}$, after the mass loss; stars with $8\leq M/M_{\odot}\leq 20$ evolved into NSs for which we assume a Gaussian distribution with peak mass $M_{\rm final}=1.33M_{\odot}$ and width $\sigma=0.12M_{\odot}$ for the end masses; stars with $M\geq 20M_{\odot}$ evolved into astrophysical BHs for which we assume a Gaussian distribution with peak mass $M=7.8M_{\odot}$ and width $\sigma=1.2M_{\odot}$. We adopt the number conservation between initial stars and stellar remnants; each massive star evolved into each stellar remnant. Under this assumption, we found the ratio of the number of each stellar remnant relative to that of main sequence stars as

[TABLE]

Throughout this paper we refer to stars with masses $0.08\leq M/M_{\odot}\leq 1$ as “main sequence stars” (MS).

Fig. 1 displays the model for the initial or final mass spectrum of BDs, MSs and stellar remnants, which we use in this paper. For the BD mass function, we will determine the normalization parameter so that the model prediction matches the OGLE microlensing events as we show below (we here adopt a normalization that is continuous with the stellar IMF at $M=0.08M_{\odot}$). The number of each lensing population determines the frequency of microlensing. Then if we focus on the event rates for a particular light curve timescale, the events arise mainly from lensing objects of the corresponding mass scales (Eq. 4). Thus, by studying the event rate as a function of the light curve timescales, one can distinguish contributions from different populations of lensing objects. Table 2 also gives the summary of our model for the mass spectrum of BDs, MS stars or stellar remnants.

Furthermore, we assume the binary fraction $f_{\rm bin}=0.4$; the fraction of MS stars or stellar remnants are in binary systems. For simplicity we consider equal-mass binary systems: we treat a microlensing of of a binary system by that of a lens with mass $M_{\rm binary}=2M$. We do not consider binary systems that contain two objects of different masses and contain two objects of different populations (e.g., MS-WD system) for simplicity. Consequently we decrease the number of lens systems from the above numbers in Fig. 1 by the binary fraction. Including the binary systems give a slightly improved agreement between the model predictions and the OGLE data, but it is not an important assumption for our main results.

To perform a calculation of microlensing event rates, we need to specify the normalization parameter of MS IMF, $A_{\rm MS}$ (Eq. 36). Recalling that the mass of Galactic bulge and disk regions is dominated by the total mass of low-mass MS stars, we determine $A_{\rm MS}$ by the condition

[TABLE]

Here $\rho_{\ast}$ is the normalization coefficient of mass density profile in the bulge and disk regions as given in Table 1. With this normalization, $A_{\rm MS}$ has a dimension of $[{\rm pc}^{-3}]$. We assume the same composition of stars and stellar remnants everywhere in the disk and bulge regions; that is, we ignore a possible dependence of the composition on a position in the Galactic region. Details of our model are different from the model in Mróz et al. (2017), so we will introduce a fudge normalization parameter later to model a possible uncertainty in the normalization: $A_{\rm MS}\rightarrow f_{\rm A}A_{\rm MS}$. However, we find $f_{\rm A}\simeq 1$, implying that our model is sufficiently close to the best-fit model in Mróz et al. (2017) or equivalently that the standard Galactic bulge/disk models are fairly accurate to reproduce the observed timescale (mass) distribution of microlensing events as we will show below.

IV Results

IV.1 OGLE data

The results shown in this paper are all based on the microlensing data taken in the OGLE-IV sky survey during the 5 years, 2011–2015 (Udalski et al., 2015; Mróz et al., 2017). The OGLE survey uses the dedicated 1.3m Warsaw Telescope, located at Las Campanas Observatory, Chile. The OGLE survey carried out a long-term monitoring observation of the nine fields towards the Galactic bulge region with a cadence of either 20 min or 60 min, covering 12.6 square degrees in total. After the careful cuts in a selection of microlensing events, the team created a catalog of 2622 microlensing events for each of which a timescale of the microlensing light curve ($t_{\rm E}$) is measured. Thus the OGLE datasets are quite rich and allow us to constrain the abundance and mass distribution of each lensing object population. As carefully studied in Mróz et al. (2017), a majority of the OGLE events can be fairly well explained by superposition of microlensing events due to BDs, MS stars, and stellar remnants. Among these contributions, the origin of BDs is not well understood as we discussed (see Burrows et al., 2001, for a review). Thanks to the unique power of microlensing that can probe a gravitational mass of a lensing object regardless of whether it is visible, the OGLE data in timescales less than $\sim 10~{}$days can be used to identify a microlensing contribution of BDs. However, the abundance and mass spectrum has not been fully understood yet. In addition to the BD contribution, those papers discussed a possible contribution of unbounded planets (wide-orbit planets or free-floating planets). Mróz et al. (2017) discussed that unbounded Jupiter-mass planets are about 0.05 – 0.25 planets per main sequence star, which is smaller than previously advocated in Ref. Sumi et al. (2011). The large OGLE dataset indicates even shorter timescale events that correspond to unbounded Earth-mass planets. Interestingly, the OGLE data also indicates a ”gap” (no microlensing event) at timescales between the BD or unbounded Jupiter-mass microlensing and the Earth-mass planets. Since planetary formation theory would predict a continuous mass spectrum, the gap, if real, seems very challenging to explain. For example, a mechanism preferentially scattering Earth-mass planets from the planetary system is needed. All these results are very interesting, and worth to further explore. All the datasets we use in this paper are taken from Extended Data Table in Mróz et al. (2017).

In this paper, to derive PBH constraints, we employ the following two working hypotheses:

• (1)

Null hypothesis of PBH microlensing: we assume that all the OGLE microlensing events are due to astrophysical objects, i.e. BDs, stars and remnants, so do not contain any PBH microlensing event. Then we use all the OGLE events to obtain an upper limit on the abundance of PBHs assuming the monochromatic mass spectrum.

• (2)

PBH hypothesis of the 6 ultrashort-timescale OGLE events, $t_{\rm E}\simeq[0.1,0.3]~{}{\rm days}$. The OGLE data found the 6 ultrashort-timescale mcirolensing events (Fig. 5), which seem a different population from the majority of events in the longer timescales due to BDs, stars and remnants. The timescale implies an Earth-mass lens. Although the origin might be unbounded Earth-mass planets, we here assume that the ultrashort-timescale events are due to PBHs, but other longer timescale events are due to astrophysical objects as in case (1). Under this hypothesis, we derive an allowed region of PBHs in two parameter space of the abundance and mass scale, assuming the monochromatic mass spectrum.

IV.2 Event rate of microlensing

We are now in a position to compute event rates of microlensing by plugging the model ingredients, which we have discussed up to the preceding subsection, into the equations such as Eq. (13).

Fig. 2 shows the expected differential number of microlensing events per logarithmic interval of the light curve timescale $t_{\rm E}$, for a single source star in the bulge region, assuming the 5-years observation as in the OGLE data. For PBH microlensing, we adopt the model ingredients in Sections III for the mass density profile and velocity distribution, assuming the monochromatic mass scale. We assumed that all DM is made of PBHs of each mass scale: $f_{\rm PBH}=1$. If we consider lighter-mass PBHs, the number density of PBHs increases and such PBHs yield a higher frequency of microlensing events with shorter timescales. In particular, for microlensing events with timescales shorter than a few days, PBHs with M_{\rm PBH}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{-1}M_{\odot} could produce a larger number of microlensing events than MS stars of $\sim 1~{}M_{\odot}$ do, if such PBHs constitute a significant fraction of DM.

In Fig. 3 we study relative contributions of MS stars in the bulge and disk regions to the total of MS microlensing events. It can be found that stars in the disk region gives a dominant contribution, while the bulge star contribution is significant for shorter timescale events.

IV.3 Comparison with the 5-years OGLE data

We now compare the model predictions of microlensing events with the 5-years OGLE data. The OGLE data contains 2622 events over the range of light curve timescales, $t_{\rm E}=[10^{-1},300]~{}{\rm days}$ (see Extended Data Table 4 in Mróz et al. (2017)). The expected number of microlensing events per a given timescale interval of $[t_{\rm E}-\Delta t_{\rm E}/2,t_{\rm E}+\Delta t_{\rm E}/2]$ is computed as

[TABLE]

where $t_{\rm obs}$ is the total observation time, $N_{\rm s}$ is the total number of source stars in the OGLE bulge fields, and $\epsilon(t_{\rm E})$ is the “detection efficiency” quantifying the probability that a microlensing event of timescale $t_{\rm E}$ is detected by the OGLE data. For the OGLE data, $t_{\rm obs}=5~{}{\rm years}$ and $N_{\rm s}=4.88\times 10^{7}$ (see Extended Data Table 2 in Mróz et al., 2017). We employ the detection efficiency, $\epsilon(t_{\rm E})$, that is taken from Extended Data Figure 2 in Mróz et al. (2017), which is explicitly shown in Fig. 4. We do not include variations of the detection efficiency in the different OGLE fields for simplicity. The coefficient $f_{\rm A}$ is a fudge normalization factor that takes into account a possible difference in details of our model calculations and the model of Mróz et al. (2017).

In Fig. 5 we compare the model prediction of microlensing event rates with the 5-years OGLE data. First of all, a majority of the OGLE microlensing events has a single peak around the timescale, $t_{\rm E}\sim 20~{}{\rm days}$, and has a gradual decrease at the shorter and longer timescales than the peak timescale. Thus the OGLE data suggests only a single population of the underlying lensing objects, except for the 6 ultrashort-timescale events, $t_{\rm E}=[0.1,0.3]~{}{\rm days}$, which we will discuss later. Interestingly, the model assuming the standard Galactic bulge and disk models (see Section III) can fairly well reproduce event rates for the main population of OGLE microlensing events. Furthermore, by employing the mass distribution of BD, stars and stellar remnants, the model can reproduce the distribution of light curve timescales (see Section III.3). Although we introduced a fudge factor in Eq. (39) to model a possible difference between our model and the model in Mróz et al. (2017), we found $f_{A}=0.99$ to have a nice agreement of our model prediction with the OGLE data at timescales greater than the peak timescale. The best-fit $f_{\rm A}$ value is close to unity, and reflects the fact that the Galactic bulge and disk model, constructed based on observations and the previous knowledges, is fairly accurate. As we discussed, the origin and nature of unbounded BDs, which cause shorter timescale events, is poorly understood. We found that the model matches the microlensing events at timescales shorter than the peak timescale, if we assume 0.18 BDs per main-sequence star (see Table 2). The figure clearly shows that MS stars with $0.08\leq M/M_{\odot}\leq 1$ give a dominant contribution to the OGLE events at timescales, t_{\rm E}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10~{}{\rm days}, while stellar remnants give secondary contributions. BDs give a dominant contribution at the shorter timescales. However, the figure shows that, as long as we assume a smooth model for the density and velocity distributions of BDs, the model cannot reproduce the ultrashort-timescale events of $t_{\rm E}\sim 0.1~{}{\rm days}$. This clearly indicates a distinct, second population of small-mass lensing objects. On the other hand, PBHs do not necessarily follow the similar timescale distribution of microlensing events to that of BD, stars or remnants, because DM has different spatial and velocity distributions from the stellar populations. As an example, Fig. 5 shows the result for case that PBHs with mass $M_{\rm PBH}=10^{-3}M_{\odot}$ are DM. A sum of the PBH and stellar population contributions give too many microlensing events compared to the OGLE data. In other words, such a PBH population is not allowed by the OGLE data. Thus we can use the OGLE data to obtain an upper bound on the abundance of PBHs with varying PBH mass scales.

IV.4 Upper bound on the PBH abundance under null hypothesis

As we showed in the preceding subsection, the Galactic bulge and disk models including the stellar components fairly well reproduce the OGLE data except for the 6 ultrashort-timescale events. In other words, the OGLE data does not necessarily imply the existence of PBH microlensing in the data. As the first working hypothesis (see Section IV.1), we here employ “null hypothesis” to obtain an upper limit on the abundance of PBHs. That is, we assume that all the observed OGLE microlensing events, including the 6 ultrashort-timescale events, are due to the stellar components, or equivalently there is no PBH lensing in the OGLE data. This would give us a most stringent upper bound on the PBH abundance. If we allow a possible PBH contribution to the OGLE data in addition to the stellar events, it would give us a more relaxed upper bound or could even allow for a detection of PBH. However, this requires a perfect knowledge of the Galactic stellar components, which is not straightforward.

We assume that the OGLE counts of microlensing events at each timescale bin follows the Poisson distribution. This is a good assumption because the same lensing object very unlikely produces multiple lensing events (lensing for multiple source stars) because of smallness of the lensing optical depth, $\tau\sim 10^{-6}$ (see Table 1). Hence we can safely assume that different microlensing events are independent and uncorrelated with each other. Under these assumptions, we assume that the log likelihood of OGLE microlensing events is given by

[TABLE]

where $N_{\rm obs}(t_{{\rm E},i})$ is the observed number of events at the $i$-th timescale bin ($t_{{\rm E},i}$); d is the data vector, ${\bf d}\equiv\{N_{\rm obs}(t_{{\rm E},1}),N_{\rm obs}(t_{{\rm E},2}),\dots,N_{\rm obs}(t_{{\rm E},n_{\rm bin}})\}$ in our case, $n_{\rm bin}$ is the number of timescale bins ($n_{\rm bin}=25$ as can be found from Fig. 5); ${\bm{\theta}}$ is the model vector; $\lambda(t_{{\rm E},i})$ is the expectation number of events at the bin. When we include PBH microlensing contributions, we model the expectation number by

[TABLE]

Here $N_{\rm exp}^{\rm PBH}(t_{{\rm E},i})$ is the expected number of microlensing events due to PBHs at the $i$-th timescale bin, which is computed from Eq. (39) once the PBH mass fraction to DM, $f_{\rm PBH}$, is specified for an assumed PBH mass scale ($M_{\rm PBH})$; $N_{\rm exp}^{\rm PBH}({{\rm E},i})\propto f_{\rm PBH}$. As a conservative approach, we use the observed counts, $N_{\rm obs}(t_{{\rm E},i})$ for the expectation value of microlensing events due to stellar components. In the following, we assume that the MW DM model for the spatial and velocity distributions for PBH in Sections III, and we treat the PBH mass fraction parameter, $f_{\rm PBH}$, as a free parameter for an assumed PBH mass scale ($M_{\rm PBH}$). Namely we consider a single model parameter for an assumed PBH mass scale (we will discuss later for a possible extension of this assumption). When $f_{\rm PBH}=0$, i.e. $N_{\rm exp}^{\rm PBH}=0$, the maximum likelihood is realized because of $N^{\rm PBH}_{\rm exp}\geq 0$. The last term in the above log likelihood is irrelevant for parameter inference, because it is a fixed number irrespectively of model parameter $(f_{\rm PBH})$.

Given the likelihood function and the PBH model (denoted as ${\cal M}$), the posterior distribution of model parameter, $f_{\rm PBH}$, is computed based on the Bayes’s theorem as

[TABLE]

where $\Pi(f_{\rm PBH})$ is a prior of $f_{\rm PBH}$ and $P({\bf d}|{\cal M})\equiv{\cal E}$ is the evidence. In this paper, we assume a flat prior, $f_{\rm PBH}\leq 1$; the total PBH mass in the MW region cannot exceed the DM mass. By computing the above equation with varying the model parameter $f_{\rm PBH}$, we can obtain the posterior distribution for an assumed mass scale of PBH. Fig. 6 shows some examples for the posterior distribution for a given PBH mass scale, obtained from the above method.

Fig. 7 shows 95% CL upper bound on the PBH abundance at different mass scales. The OGLE data improves the constraints on the abundance for PBHs in the mass range $M_{\rm PBH}\simeq[10^{-3},10^{-6}]$. The dashed curve shows the upper bound if we do not include the 6 ultrashort-timescale events in the 4 shortest timescale bins for the null hypothesis. Thus the upper bound at small mass scales is sensitive to the assumption of whether we include the short timescale bins in the analysis. The results can be compared with other constraints such as those from the Subaru Hyper Suprime-Cam (HSC 333https://hsc.mtk.nao.ac.jp/ssp/) observation of Andromeda galaxy (M31) (Niikura et al., 2017) and the earlier MACHO/EROS experiments (Alcock et al., 2000; Tisserand et al., 2007). The OGLE bound is stronger than that of MACHO/EROS due to the larger sample of microlensing events. The OGLE constraint is complementary to the Subaru HSC result that uses even denser (2 min) cadence data of M31 to search for microlensing for a larger number of source stars, but from only a single night observation. If we want to extend the constraint to PBHs at larger mass scales, we need to use the microlensing data extending to longer timescales such as year timescales.

IV.5 A possible detection of Earth-mass scale PBHs from short-timescale OGLE data

Now we employ the second working hypothesis in Section IV.1. That is we consider a case that the 6 ultrashort-timescale OGLE events, in $t_{\rm E}\simeq[0.1,0.3]~{}{\rm days}$, are due to PBH microlensing. In this case, we assume that the expectation number of microlensing events at each of the first 4 timescale bins (in the ultrashort-timescale bins) is given as

[TABLE]

for the Poisson distribution of micrlensing event counts (Eq. 40). Here we should again note that we assume the monochromatic mass scale for PBHs. We also assume that the OGLE events in the longer timescales, i.e. the majority of OGLE events, are not due to PBHs (i.e no PBH microlensing as in the preceding section).

Fig. 8 shows a 95% CL allowed region of PBHs in two parameter space of its mass and abundance. The thin-blue shaded region corresponds to the allowed region obtained from the OGLE data alone, while the thick shaded region is the allowed region when combining the OGLE results with the null hypothesis of PBH lensing for the longer timescale OGLE events (t_{\rm E}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}0.5~{}{\rm days}) and the the HSC constraints for M31. The best-fit model, which has a maximum likelihood, is a model with $M_{\rm PBH}=9.5\times 10^{-6}M_{\odot}$ and $f_{\rm PBH}=0.026$. The figure shows that PBHs of Earth-mass scales ($3\times 10^{-6}M_{\odot}$) can well reproduce the 6 ultrashort-timescale OGLE events if the abundance is at a per cent level. This allowed region is also consistent with null PBH results in the HSC data and the longer timescale OGLE data (if assuming the monochromatic mass spectrum).

Fig. 9 compares the best-fit PBH model with the timescale distribution of OGLE events in the shortest timescales. Interestingly, the width of OGLE timescale distribution is nicely reproduced by the velocity distribution of PBHs in the MW DM model. For comparison, the two dashed curves denote the model predictions for two models that are close to the boundary of 95% C.L. intervals in the allowed region, which are specified by model parameters $(M_{\rm PBH}/M_{\odot},f_{\rm PBH})=(1.6\times 10^{-6},0.062)$ or $(6.9\times 10^{-5},0.014)$, respectively. These failed models under- or over-predict the microlensing event rats over the range of lightcurve timescale bins, and also become inconsistent with the upper bounds of the HSC M31 data and the longer timescale OGLE data. This result implies that, if PBHs have a wide mass spectrum extending to larger masses than the best-fit mass, such a model generally fills the gap around $t_{\rm E}\simeq 0.4~{}{\rm days}$ between the ultrashort-timescale events and the main population. Note that PBH models with smaller masses than the best-fit mass generally predict too many microlensing events at even shorter timescales $t_{\rm E}\leq 0.1~{}{\rm days}$, however, the OGLE data does not have a sensitivity due to the limitation of the cadence data (20 min cadence), and therefore the OGLE data has no sensitivity to events at $t_{\rm E}\leq 0.1~{}{\rm days}$, which is taken into account by the detection efficiency in Fig. 4. The original event rates can have an increasing function at t_{\rm E}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.1~{}{\rm days}. These mass-scale PBHs are well constrained by the HSC results.

As we have shown, the PBH model assuming the MW DM model can give an alternative explanation of the ultrashort-timescale OGLE events. Since there is an uncertainty in the MW DM model, especially the DM distribution around the halo center, one might worry whether such an uncertainty in the PBH distribution around the halo center is sensitive to our results. Fig. 10 shows which distant PBHs for the best-fit model contribute the microlensing events at timescale $t_{\rm E}=0.2~{}{\rm days}$. Due to lensing efficiency, PBHs over a wide range of distances between the Earth and the Galactic center equally contribute the lensing events. The PBHs near the Galactic center, where the DM distribution is most uncertain, is not particularly sensitive to the final result. This lens distance dependence is contrary to DM annihilation that is quite sensitive to details of the DM density in the Galactic center.

V Discussion and Conclusion

In this paper we have used the largest sample of microlensing events for stars in the Galactic bulge, obtained from the 5-years OGLE observation, to constrain the abundance of PBH that is a viable candidate of DM in the MW region. The 2622 microlensing data contains rich information on event rates and light curve timescales that correspond to the abundance and mass scale of lensing “compact” objects, because kinematical or velocity structures of lensing objects, even for DM in the MW halo region, are well constrained by various observations; the relative velocity of lens-source-observer determines the crossing time of the lensing Einstein radius, or equivalently the microlensing light curve timescale, once mass of a lensing object is assumed because the mass determines the Einstein radius. Thus we can use the invaluable OGLE data, which covers the wide range of timescales $t_{\rm E}\simeq[0.1,300]~{}{\rm days}$, allows us to explore the abundance of unknown “compact” objects over the wide range of mass scales.

To do this, we first revisited the Galactic bulge/disk models to estimate event rates of microlensing due to astrophysical objects such as brown dwarfs, MS stars, and stellar remnants (white dwarfs, neutron stars, and astrophysical black holes), following Ref. Mróz et al. (2017). Since the mass of Galactic disk/bulge regions is dominated by low-mass main-sequence (MS) stars around $M\sim 0.5M_{\odot}$, astronomers have a quite good knowledge of the abundance of stars from various observations of star number counts. In addition, stellar remnants are from massive stars, so we can infer their abundances from the initial mass function of low-mass MS stars that is well constrained by various observations. We showed that, even if details of our models would be different from that of Ref. Mróz et al. (2017), the standard Galactic bulge/disk models including MS stars and stellar remnants nicely reproduce the OGLE events at timescales t_{\rm E}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}20~{}{\rm days}, corresponding to objects with M\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}1M_{\odot}. For the shorter-timescale events, we need to add contributions from brown dwarfs (BD) that are invisible or difficult to directly observe. The origin and nature of BDs (0.01\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}M/M_{\odot}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.08) are not well understood. Some of BDs should form around a primary MS star, while some of BDs would form, as a primary gravitating object, in a protoplanetary disk. The BDs contributing the short timescale microlensing are “unbounded” BDs, because the timescale becomes too long if the host primary star contributes microlensing (because of much larger mass compared to that of BD). Nevertheless, the timescale distribution of OGLE data continuously extends to shorter timescales than the MS peak timescale, and it suggests a population of BD-mass objects with a continuous abundance to that of MS stars, except for the 6 ultrashort-timescale events in $t_{\rm E}\simeq[0.1,0.3]~{}{\rm days}$. We showed that our model can reproduce the entire timescale distribution for the main population of OGLE events if we adjust the abundance of BDs as done in Mróz et al. (2017) (about 0.18 BDs per MS star in our model).

Given a justification of the standard Galactic bulge/disk models, we employ the “null hypothesis”, i.e. no PBH microlensing event in the OGLE data, to obtain the stringent upper bound on the abundance of PBHs in the mass range $M_{\rm PBH}\simeq[10^{-6},10^{-3}]M_{\odot}$ (from Earth to Jupiter mass scales), assuming the monochromatic mass spectrum (Fig. 7). The upper bounds are tighter than the previous bound from the MACHO/EROS experiments (Tisserand et al., 2007) and the Subaru/HSC microlensing search for M31 (Niikura et al., 2017). This result shows the power of microlensing for exploring the PBH abundance.

Even more interestingly, we showed that the 6 ultrashort-timescale events can be well explained by PBHs of Earth-mass scales if such PBHs constitute about 1% of DM in the MW region (Fig. 8). Even if we employ the monochromatic mass spectrum for simplicity, the timescale distribution naturally arises from the velocity distribution of PBHs expected for the DM kinematical structures in the MW region. There is a mechanism in inflation model to produce PBHs in such a narrow mass range or with right abundance (e.g. Inomata et al., 2018). If this is a real PBH microlensing, it would be a big discover. Such a small mass black hole cannot be made by any astrophysical process, so this would also give an evidence of the large primordial perturbations at the corresponding Hubble horizon in the early universe. If we include a possible distribution of PBH masses, the results would be changed. Nevertheless it is rather straightforward to translate our results into a specific PBH model with a given mass spectrum, following the methods developed in Carr et al. (2016) (also see Inomata et al., 2017b).

A usual explanation of the ultrashort-timescale OGLE events is due to “unbounded” Earth-like planets, where “unbounded” is needed to have the right microlensing timescale. These unbounded, more exactly wide-orbit or free-floating planets, could be formed by the formation of planetary system or scattering of planetary systems; for example, if a planetary system encounters a massive planet or star, such an Earth-mass planet might be scattered. Even if this happens, why there is a gap around $t_{\rm E}\simeq 0.4~{}{\rm days}$ between the ultrashort-timescale events and the majority of events (i.e. the main population). Of course the gap might be an apparent statistical fluctuation due to the low number statistics, but it would be completely a mystery if if this is genuine, because there would be a continuous mass spectrum expected for unbounded planets from Earth (or even smaller-masses) to Jupiter masses. These involve complicated, nonlinear astrophysics in planetary or star formation, so a further observational study would be a more direct path to resolving the origin and nature of these short timescale events. For example, a more detailed study of lightcurve for each ultrashort-timescale event would be very useful (see Mroz et al., 2018, for such an attempt). This requires a denser cadence data of microlensing search. Such a denser cadence data would be also useful to distinguish genuine microlensing events from other contaminating events such as stellar flare.

To confirm or falsify the PBH hypothesis of short-timescale events against free-floating planets, there is a very promising, robust way. It is a microlensing search for stars in M31 using the Subaru HSC data (or eventually LSST data towards the Magellanic Clouds). The angular direction of M31 is in a high latitude in the Galactic coordinates, i.e. far from the Galactic disk. If the short timescale microlensing events are due to free-floating planets, we should expect a much smaller number of events towards M31, because there is a much less number of stars in the high latitude direction compared to the direction to the Galactic center. On the other hand, if the PBH scenario is true, we should find microlensing events with a frequency predicted by the standard DM model of the MW halo region that is supported by the disk rotation measurements. As shown in Niikura et al. (2017), the Subaru and HSC combination is ideal because its FoV can cover the entire disk region of M31 and the large aperture allows us to use main-sequence stars in M31 for the microlensing search even with short exposure (e.g., 90 sec). We are now carrying out a monitoring observation of M31 with Subaru HSC, and we envision that we can address these important questions in the near future.

Acknowledgments

We would like to thank the OGLE collaboration for making their microlensing data publicly available to us via Mróz et al. (2017). We also thank to Chiaki Kobayashi, Toshiki Kurita, Hitoshi Murayama, Toyokazu Sekiguchi, Sunao Sugiyama and Misao Sasaki for useful discussion. This research was supported by World Premier International Research Center Initiative (WPI), MEXT, Japan. This work was in part supported by MEXT Grant-in-Aid for Scientific Research on Innovative Areas (No. JP15H05887, JP15H05893, JP15K21733, JP18H04356), and by JSPS KAKENHI Grant Numbers JP17J03653 (HN), JP15H03654 (MT), and 15K17659 (SY).