Effect of clustering on primordial black hole microlensing constraints

Stellar microlensing observations tightly constrain compact object dark matter in the mass range (10-11–103)M ⊙. Primordial Black Holes (PBHs) form clusters, and it has been argued that these microlensing constraints are consequently weakened or evaded. For the most commonly studied PBH formation mechanism, the collapse of large gaussian curvature perturbations generated by inflation, the clusters are sufficiently extended that the PBHs within them act as individual lenses. We find that if the typical mass of the clusters is sufficiently large, ≳ 106 M ⊙, then the event duration distribution can deviate significantly from that produced by a smooth dark matter distribution, in particular at the shortest durations. As a consequence of this, the probability distribution of the number of observed events is non-Poissonian, peaking at a lower value, with an extended tail to large numbers of events. However, for PBHs formed from the collapse of large inflationary perturbations, the typical cluster is expected to contain ∼ 103 PBHs. In this case the effect of clustering is negligibly small, apart from for the most massive PBHs probed by decade-long stellar microlensing surveys (M PBH ∼ 103 M ⊙).


Introduction
The discovery of gravitational waves from mergers of tens of Solar mass black holes by LIGO-Virgo [1] has led to increased interest in Primordial Black Holes (PBHs) as a dark matter (DM) candidate [2][3][4][5].PBHs are black holes that may form in the early Universe [6,7].The most commonly studied formation mechanism is the collapse of large density perturbations produced by inflation (for reviews see e.g.Refs.[8,9]).
Stellar microlensing is the temporary amplification which occurs when a compact object passes close to the line of sight to a star [10].Various microlensing surveys have placed tight constraints on the abundance of compact objects in the Milky Way (MW) halo.The OGLE Galactic bulge survey [11] and observations of M31 using Subaru HSC [12,13] constrain planetary and sub-planetary masses, while the EROS [14], MACHO [15] and OGLE [16] surveys of the Large and Small Magellanic Clouds (LMC and SMC) constrain stellar and planetary masses.Following a proposal by Ref. [17], Ref. [18] has combined data from EROS-2 and MACHO to obtain sensitivity to long duration events, and hence constrain more massive compact objects.Taken at face value, the stellar microlensing constraints exclude PBHs with mass 10 −11 M M PBH 10 3 M making up all of the DM.However the calculation of these constraints involves various assumptions, for instance that the DM is smoothly distributed.
PBHs that form from the collapse of large gaussian perturbations generated by inflation do not form in gravitationally bound clusters [19].However, since PBHs are discrete objects, there are Poisson fluctuations in their initial distribution.As a consequence of these isocurvature fluctuations in the PBH density, PBH clusters form not long after radiation-matter equality [20].The abundance and properties of these clusters have been studied numerically [21] and analytically [22], using the spherical top-hat collapse model.Refs.[23][24][25] have argued that PBH clustering modifies the stellar microlensing constraints so that they are shifted to lower masses, and consequently multi-Solar mass PBHs can make up all of the DM.
In this manuscript we examine the effect of clustering on the stellar microlensing constraints for the clusters which form when PBHs are produced by the collapse of large gaussian inflationary density perturbations.In Sec. 2 we overview the properties of the PBH clusters.Next, in Sec. 3, we outline the calculation of the microlensing differential event rate, firstly for the standard case of a smooth MW halo (Sec.3.1) and then for clustered DM (Sec.3.2).We present our results in Sec. 4 before concluding with discussion in Sec. 5.

Cluster properties
Jedamzik used the spherical top-hat collapse model to calculate the properties of the PBH clusters that form when PBHs generated from the collapse of large inflationary density perturbations, with a single mass, make up all of the dark matter [22].The initial fluctuation in the density in a region containing N PBHs is δ(N ) = 1/ √ N , and these isocurvature fluctuations grow with time proportional to [20] where a eq is the scale factor at radiation-matter equality.A particular scale goes non-linear when the scale factor is equal to a coll , determined by D(a coll )δ(N ) ≈ δ c ≈ 1.68 and the cluster radius, R cl , can be estimated, from For initially Poisson distributed discrete objects, the number of clusters containing N cl objects, Ñ , is given (for N cl 1) by [21,26] where δ (a) = δ c /D(a) and (2.5) Ñ is always a monotonically decreasing function of N cl , and N grows with time.Therefore clusters containing a small number of PBHs are always the most abundant, however the number of clusters with large N cl increases with time.This behaviour has been confirmed numerically [21].Clusters containing small numbers of objects evaporate [27], and PBH clusters with N cl 10 3 will have evaporated by the present day [20,22].Therefore, for PBHs that form from large inflationary density perturbations, the most common cluster size today is expected to be N cl ∼ 10 3 , independent of the PBH mass.

Microlensing event rate
In this section we outline the calculation of the microlensing differential event rate, first for the standard case of a smooth halo (Sec.3.1) and then for clustered DM (Sec.3.2).

Smooth halo
The microlensing differential event rate, dΓ/d t, towards the LMC for a smooth halo composed entirely of compact objects with mass M PBH and a Maxwellian velocity distribution is given by [28,29]: where t is the time taken to cross the Einstein diameter, R E (x) is the Einstein radius G is the Gravitational constant, L ≈ 50 kpc is the distance to the LMC, x is the distance of the lens from the observer in units of L and ), where v c = 220 km s −1 is the circular speed.
The standard halo model usually assumed in microlensing studies ('Model S')2 is a cored isothermal sphere with density profile and local dark matter density ρ 0 = 0.0079M pc −3 , core radius R c = 5 kpc and Solar radius R 0 = 8.5 kpc.The differential event rate, Eq.(3.1), is then given by [29] where The expected number of events, N exp , is given by where E is the exposure in star years and ( t) is the detection efficiency i.e. the probability that a microlensing event that occurs with duration t is detected.

Clustered halo
The typical separation of PBHs in a cluster is much larger than the Einstein Radius (for Therefore the individual PBHs act as lenses, and not (as argued in Ref. [23,24]) the cluster as a whole.Appendix A2 of Ref. [25] argues that even when clusters are sufficiently diffuse that the PBHs act as individual lenses, lensing by the cluster as a whole renders the magnification from lensing by a single PBH unobservable.Their argument, however, relies on a significant underestimate of the Einstein radius (see Appendix A for further details).A fraction of the PBHs may be in binaries [22,32,33].Ref. [34] has however shown that the time separation of the lensing events caused by PBHs in a binary would be of order 100 yr, and hence the PBHs act as separate, individual lenses.
Our method for calculating the microlensing event rate from clusters is similar to Refs.[35,36].We assume that the surface area of the LMC is circular, so that microlensing events can be caused by compact objects within a cone with apex at the Earth and base at the LMC.We take the cone half angle to be θ = 5.2 • , to match the 84 deg 2 of the LMC monitored by EROS-2.Using Eq. (3.4) for the density profile of the MW, the total mass of DM in this cone is M cone ≈ 9 × 10 8 M .We assume that a fraction f of the DM is in the form of PBHs3 , and all PBHs are in clusters containing N cl PBHs with mass M PBH .We saw in Sec. 2 that clusters with N cl 10 3 will have evaporated by the present day, and therefore some (probably quite large) fraction of PBHs will not be clustered today.Therefore assuming that all PBHs are in clusters with N cl 10 3 provides an upper limit on the actual effect of clustering on the EROS-2 microlensing constraints.As noted by Petač et al. [34], for large N cl and/or M PBH the cluster radii, R cl , given by Eq. ( 2.3) from the spherical top-hat collapse model are unphysically large.In particular they are larger than the typical separation between clusters.Therefore we follow Ref. [34] and set R cl = 10 pc.
In order to take into account clusters that lie only partly within the microlensing cone, we simulate clusters within a larger region which is centered on the microlensing cone and has radius at each line of sight distance, x, equal to the radius of the microlensing cone plus the cluster radius: r tcone (x) = xL tan θ + R cl , i.e. a truncated cone with the narrow end at the Earth.For each combination of M PBH and N cl we first calculate the average number of clusters within the truncated cone, N tcone , where M tcone is the mass within the truncated cone.For each realisation we first draw the actual number of clusters from a Poisson distribution.The line-of-sight position, x cl , and transverse velocity, v ⊥,cl , of each cluster are generated assuming the cored isothermal sphere density profile, Eq. (3.4), and a Maxwellian velocity distribution with v c = 220 km s −1 .We also generate a value for the distance of the centre of the cluster from the axis of the microlensing cone such that, at each x, the clusters are uniformly distributed within the circular crosssection of the truncated cone.This distance is then used to calculate f , the fraction of the cluster within the microlensing cone.
The velocity dispersion of PBHs within a cluster is of order [22] σ cl ≈ 0.6 This is negligible compared with the cluster transverse velocity, and therefore all PBHs within a given cluster will cause microlensing events with the same duration ≈ 300 Next we need to calculate the rate at which lensing events occur for each cluster.The optical depth is the probability that a star lies within the Einstein radius of a lens (e.g.Ref. [37]).For a cluster which lies entirely within the microlensing cone the optical depth, τ cl , is the product of the lensing cross section (πR 2 E ), the surface number density of lenses and the fraction of the solid angle to the LMC, Ω LMC , covered by the cluster [36]: where is the solid angle subtended by the cluster.In a time dt the lensing area swept out by a lens is dA = 2R E v ⊥,cl dt and hence the probability of a new microlensing event occurring is The rate at which microlensing occurs, Γ cl = dτ cl /dt, is therefore where f is the fraction of the cluster which lies within the microlensing cone.
For each realisation, we calculate the total differential event rate, dΓ/d t, from all clusters by summing the binned values of the event durations, tcl , for each cluster, weighted by their rates, Γ cl .The mean number of events produced by each cluster is given by Ncl = E ( tcl )Γ cl . (3.15) For each cluster we draw the observed number of events, N cl,obs , from a Poisson distribution with mean Ncl .The total number of observed events, N obs , is the sum of N cl,obs over all clusters.

Results
We use the method described in Sec.3.2 to calculate the differential event rate for 10 4 realisations of each combination of the number of PBHs in a cluster, N cl , and the PBH mass, M PBH .For cases where the number of clusters in the cone to the LMC is large, 10 4 , then the differential event rate for a single realisation only has the expected small stochastic deviations from the differential event rate produced by smoothly distributed DM.However when the number of clusters in the cone is smaller than this, there are systematic deviations in the differential event rate for short-duration events.For most realisations there is a deficit of short-duration events, however for a small fraction of realisations there is a large excess of short events.Fig. 1 shows the differential event rate, dΓ/d t, for three different realisations for N cl = 10 6 and both M PBH = 1 and 10 M , compared with the case of smoothly distributed DM.For both values of M PBH we show two 'typical' realisations, which have a deficit of short events, and one 'rare' realisation with an excess.This behaviour can be understood by considering the dependence of the Einstein radius, R E , and the cross-sectional area of the cone to the LMC on x, the fractional distance along the line of sight.The Einstein radius is proportional to [x(1 − x)] 1/2 while (for the standard halo model) the lens transverse velocity distribution is independent of x.Therefore short-duration events are typically produced by lenses (in the case of clustered DM, clusters) at small or large x.The cross-sectional area of the cone to the LMC is proportional to x 2 , therefore the probability of there being a cluster within (or partly within) the cone at small x is small.Most realisations don't have clusters at very small x, and hence have a deficit of short-duration events.For the small fraction of realisations which do have a cluster at very small x, that cluster subtends a large fraction of the solid angle to the LMC and hence produces a high rate of short-duration events.More quantitatively, see Eq. (3.14), the total lensing rate by the cluster, Γ cl , is proportional to R E /x 2 cl .The 'rare' realisations in Fig. 1 both have a cluster with a small x cl value.
These variations in the rate of short events emerge when the number of clusters in the cone to the LMC is smaller than ∼ 10 3 (which corresponds to a number of PBHs per cluster N cl 10 6 (M /M PBH )) and become larger if the number of clusters is decreased.Since the Einstein radius increases with increasing M PBH , so does the value of t at which the variations in the differential event rate appear.We note that for the standard PBH formation mechanism, the collapse of large inflationary density perturbations, most clusters are expected to have N cl ∼ 10 3 , and not all PBHs are in clusters.Therefore, for this formation mechanism, we expect this effect to be negligible apart from for the most massive PBHs probed by stellar microlensing, M PBH ∼ 10 3 M .
Next we study the effect of these variations in the differential event rate on the number of events predicted in LMC microlensing surveys.We consider two different microlensing survey configurations: • An EROS-2-like survey, with exposure E = 3.77 × 10 7 star years and detection efficiency, ( t), given in Fig. 11 of Ref. [14], which observes no microlensing events.
• A 'toy' long-duration event survey, with E = 2.5 × 10 9 star years and ( t) = 0.4 for 400 day < t < 15 years and zero otherwise, which observes no microlensing events.
For the later survey we have chosen the exposure and maximum event duration to, roughly, mimic catalogues 2 and 3 in Ref. [17].The minimum event duration matches the cut-off imposed in Ref. [18] to remove backgrounds from lensing by stars in the LMC or MW disk, and the efficiency roughly matches that obtained in their analysis.Fig. 2 shows the probability distribution of the observed number of events, P [N obs (f = 1)], if all of the MW halo is in clusters containing a fixed number of PBHs for i) M PBH = 1M and the EROS-2 like survey and ii) M PBH = 10 3 M and the 'toy' long-duration event survey.For the former we consider both the EROS-2 detection efficiency, and also perfect detection efficiency, ( t) = 1 for all t.We see that if the number of clusters in the cone to the LMC is less than of order a thousand, the probability distribution deviates from the Poisson distribution expected for smoothly distributed DM; the peak of the distribution is shifted to a smaller value of the number of events, and there is an extended tail to large numbers of events.This behaviour is a direct consequence of the variations in the differential event rate for the shortest events discussed above.When the number of clusters is not large most realisations have a deficit of short events and hence a lower observed number of events than for smoothly distributed DM, while a small fraction of realisations have an excess of short events and hence a high observed number of events.For M PBH = 1M the deviation of the probability distribution from Poissonian only emerges for N cl 10 6 , much larger than the typical size of clusters for the standard PBH formation mechanism (N cl ∼ 10 3 ).The deviations from the Poisson distribution are smaller for the EROS-2 detection efficiency than for perfect efficiency, because the EROS-2 efficiency is largest for t ≈ 200 days (and for M PBH = 1M the variations in the event duration distribution manifest at smaller values of t where the efficiency is smaller).For M PBH = 10 3 M the deviations are visible, but relatively small, for N cl = 10 3 .
Finally we study the effect of clustering on the constraints on the fraction of the MW halo in PBHs, f .For a survey which observes zero events, a 95% confidence limit on the PBH halo fraction can be calculated, as in Ref. [14], by finding (for each value of M PBH ) the value of f for which P [N obs (f ) = 0] = 0.05.For smoothly distributed DM, P [N obs (f )] is Poissonian and hence The differential event rate is directly proportional to the local dark matter density, ρ 0 , and therefore the expected number of events for smoothly distributed DM is directly proportional to f : N exp (f ) = f N exp (f = 1), where N exp (f = 1) is the expected number of events for f = 1, calculated using Eq.(3.6).Setting Eq. (4.1) equal to 0.05 gives N exp (f ) = 3.0 and therefore f = 3.0/N exp (f = 1).The constraints on f for smoothly distributed DM obtained for the 'EROS-2-like' survey match those found by the EROS collaboration [14] to within ∼ 10% (e.g.Ref. [38]).For clustered DM (if the number of clusters within the cone is small) the probability distribution of the observed number of events is non-Poissonian, and the 95% exclusion limit on f has to be found by explicitly calculating P [N obs (f )] for a range of f values, to find the value of f for which P [N obs (f ) = 0] = 0.05.
For the EROS-2 like survey the change in the constraint on the halo fraction is negligible for all values of M PBH for which f < 1, unless N cl is many orders of magnitude larger than expected for the standard PBH formation mechanism.For the 'toy' long-duration survey the change in the constraint is only non-negligible (for physically relevant values of N cl ) for large values of M PBH .The 95% confidence limit on the halo fraction in PBHs with M PBH = 10 3 M is f < 0.076 for smoothly distributed DM.For N cl = 10 3 the increased probability of small values of N obs leads to a weakening of the constraint to f < 0.096.

Discussion
We have revisited the constraints on PBH DM from stellar microlensing towards the LMC, taking into account the clustering of PBHs expected when PBHs form from the collapse of large gaussian perturbations generated by inflation.In this case the PBH clusters are sufficiently diffuse that the PBHs act as individual lenses, and clusters containing N cl ∼ 10 3 are expected to be most abundant, with smaller clusters having evaporated.For simplicity we assume that all PBHs have the same mass, M PBH , and are in clusters containing a fixed number of PBHs, N cl .In fact some fraction of the PBHs, including those that were previously in clusters with N cl 10 3 , will be unclustered today, and therefore our results provide an upper limit on the effect of clustering on the LMC stellar microlensing constraints.
We find that if the number of clusters in the cone to the LMC is sufficiently small, 10 3 , or equivalently the number of PBHs in each cluster, N cl , is greater than 10 6 (M /M PBH ), then the differential event rate deviates significantly from that produced by a smooth halo for short event durations.This is because the probability of there being a cluster close to the observer is small, however if there is such a cluster it produces short-duration events at a high rate.Consequently most realisations, which don't have a cluster close to the observer, have a deficit of short events (see top two rows of Fig. 1).However the rare realisations where there is a cluster close to the observer have a high rate of short events (see bottom row of Fig. 1).
Consequently, as shown in Fig. 2, the probability distribution of the observed number of events deviates from the Poisson distribution produced by a smooth DM distribution.It peaks at a smaller value, since most realisations have a deficit of short events, and has a long tail to large values, from the rare realisations with a cluster close to the observer which produces a high rate of short events.However the number of clusters is only small enough for these effects to occur if either the number of PBHs per cluster, N cl , is larger than expected, and/or the PBH mass is large.Even for the most massive PBHs probed by decade long microlensing surveys (M PBH ∼ 10 3 M ), the change in the constraint on the halo fraction in PBHs, f , is only of order ten-percent if all of the PBHs are in clusters with N cl ∼ 10 3 (in fact not all of the PBHs are expected to be in clusters).
In summary, PBH clustering could have a significant effect on stellar microlensing constaints if the clusters are sufficiently compact (so that the cluster as a whole acts as a lens) or have a sufficiently large mass (so that the number of clusters in the cone to the LMC is small, 10 3 ).However for the most commonly studied PBH formation mechanism, the collapse of large gaussian perturbations generated by inflation, the clusters are expected to be diffuse enough that the PBHs act as lenses individually, and the number of PBHs in a typical cluster sufficiently small (N cl ∼ 10 3 ), that the change in the constraints on the PBH abundance is small, even for the most massive PBHs probed by decade-long microlensing surveys.
While we were completing this work similar work by Petač et al. [34] appeared on the arXiv.They use a different method and take into account some effects that we neglect (e.g. the variation of the surface density of stars in the LMC and the density profile of the PBH clusters).Nonetheless our results for the probability distribution of the observed number of events (see Fig. 2) are in very good agreement with theirs.In addition we have shown that the variations in this probability distribution arise from the effect of rare clusters close to the observer on the rate of the shortest duration events.Also, following the appearance of Ref. [18], we have looked explicitly at the more massive PBHs, M PBH up to 10 3 M , probed by long-duration microlensing surveys.
cos b cos l and b = −32.8• and l = 281 • are the galactic latitude and longitude, respectively, of the LMC.

2 Figure 1 .
Figure1.Example realisations of the differential event rate, dΓ/d t, for clustered DM (blue lines) compared with the standard smooth DM halo (black).In all six cases all of the DM is in clusters containing cl = 10 6 PBHs and the PBHs have mass M PBH = 1 and 10M in the left and right hand columns respectively.The top two rows show 'typical' realisations, where the absence of any clusters close to the observer leads to a deficit of short-duration events.The bottom row shows examples of rare realisations where there is a cluster close to the observer which produces short-duration events at a high rate (note the different range of the y-axis in this case).

Figure 2 .
Figure 2. The probability distribution, P [N obs (f = 1)], of the observed number of microlensing events, N obs (f = 1), if all of the MW halo is in compact objects with M PBH = 1M (left panel) and M PBH = 10 3 M (right panel).Left panel: The orange, red and brown lines assume all of the DM is in clusters with M cl = 10 6 , 10 7 and 10 8 M respectively.The solid lines use the EROS-2 detection efficiency (see Sec. 3 for further details) while the dotted lines assume perfect efficiency, i.e. ( t) = 1 for all t.The black dashed and dotted lines show the Poisson distribution, which arises in the standard case of a smooth DM halo, with N exp (f = 1) = 25 and 60 for the EROS-2 and perfect detection efficiencies respectively.Right panel: The green, blue and purple lines assume all of the DM is in clusters with M cl = 10 3 , 10 4 and 10 5 M respectively for a 'toy' long-duration event survey (see Sec. 3).The black dashed line shows the Poisson distribution, which arises in the standard case of a smooth DM halo, which has N exp (f = 1) = 40.