{'fraction_non_alphanumeric': 0.030033632527152375, 'fraction_numerical': 0.008611640710353249, 'mean_word_length': 4.65568209353041, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Maximal-Density Mass Function for Primordial Black Hole Dark Matter

The advent of gravitational wave astronomy has rekindled interest in primordial black holes (PBH) as a dark matter candidate. As there are many different observational probes of the PBH density across different masses, constraints on PBH models are dependent on the functional form of the PBH mass function. This complicates general statements about the mass functions allowed by current data, and, in particular, about the maximum total density of PBH. Numerical studies suggest that some forms of extended mass functions face tighter constraints than monochromatic mass functions, but they do not preclude the existence of a functional form for which constraints are relaxed. We use analytical arguments to show that the mass function which maximizes the fraction of the matter density in PBH subject to all constraints is a finite linear combination of monochromatic mass functions. We explicitly compute the maximum fraction of dark matter in PBH for different combinations of current constraints, allowing for total freedom of the mass function. Our framework elucidates the dependence of the maximum PBH density on the form of observational constraints, and we discuss the implications of current and future constraints for the viability of the PBH dark matter paradigm.


Introduction
The possibility that density fluctuations in the early universe collapsed into primordial black holes (PBH) has been studied for several decades [1]. Apart from their potential utility as a probe of the primordial universe, PBH are an excellent candidate for cosmological dark matter, as sufficiently large black holes are stable and dynamically cold. Further, with simple formation mechanisms, they can be produced with a cosmological density matching the observed density of dark matter.
If PBH account for a significant fraction of dark matter, it is possible that observed gravitational wave signals have a primordial origin. Direct observations of binary black hole mergers thus far all involve black holes with masses of several times 10M [2][3][4][5][6], in a range where microlensing constraints on the abundance of compact objects are ineffective. The observed merger rate is compatible with PBH as dark matter, and other constraints historically applied in the LIGO mass range are subject to large astrophysical uncertainties [7,8]. This has led to renewed interest in primordial production mechanisms, and it remains possible that PBH in this mass window account for much or all of dark matter [9].
However, depending on the formation mechanism, PBH may exist today with masses as small as 10 −16 M , or as large as those of supermassive black holes. Thus, constraining the total density contained in PBH requires the combination of constraints that span this vast range of mass scales. Such observables include microlensing surveys [10][11][12][13], CMB data [14], and the statistics of wide binaries [15]. In general, constraints from these observables have been computed under the assumption that all PBH have the same mass. The corresponding mass functions, comprising a single Dirac delta, are said to be monochromatic. However, as realistic production mechanisms necessarily result in an extended (non-monochromatic) mass function, it is essential to correctly combine constraints across all masses.
This problem has recently been studied by several authors [16][17][18]. In general, the constraints depend non-trivially on the functional form of the mass function, and statements about the implications of constraints for properties of the PBH population can be difficult to generalize. In particular, the total fraction f PBH of dark matter that may be accounted for by PBH varies with the form of the mass function, so f PBH = 1 is ruled out for some forms of the mass function, and allowed for others. This has led to confusion regarding the observational viability of the PBH dark matter scenario, and while prior work has established procedures for comparing specific extended mass functions with observables, general statements regarding the allowed total fraction of dark matter in PBH are lacking.
Depending on the set of constraints considered, observational data may or may not already rule out f PBH = 1 for monochromatic mass functions. Since the many constraints span a wide mass range, and since several do not overlap significantly, some authors have argued that broadening the mass function might relax constraints on PBH [19,20], possibly allowing for f PBH = 1 even if that possibility were excluded by constraints for monochromatic mass functions. However, [17,18] have evaluated the constraints numerically for several forms of extended mass functions, and found that extended mass functions are typically subject to stronger constraints than monochromatic mass functions.
These findings motivate the question we now pose: what is the theoretical maximum density of PBH permitted by constraints for a fully general mass function? Our goal is ultimately to clarify the observational status of PBH dark matter, and to understand the circumstances under which extending the mass function can relax constraints. We also seek a procedure which is flexible and simple enough to allow us to compare results for different sets of constraints, and to elucidate the dependence of the maximal density on the form of the constraints themselves. To that end, we derive the form of the mass function which optimizes the density subject to all observational constraints combined. This allows us to obtain a general bound on the density of PBH with minimal numerical computation, independently of the true form of the PBH mass function. Note that we do not propose a new prescription for the evaluation of constraints for a given extended mass function. Rather, we maximize the PBH density subject to constraints as evaluated using existing methods from the literature. The maximal-density mass functions we derive then provide insights into the overall impact of each individual observable. This paper is organized as follows. In section 2, we establish conventions and notations, and review the application of constraints from the monochromatic case to extended mass functions. In section 3, we present a pedagogical derivation of our main results regarding the maximum density of PBH, and we apply them to current data. We consider the impact of gravitational wave constraints separately in section 4. We discuss these results in section 5 and conclude in section 6. Finally, in appendix A, we validate our analytical results with direct numerical techniques.
2 Interpreting constraints for extended mass functions

The interpretation problem
Applying observational constraints to generic extended mass functions is non-trivial. It is not sufficient to check that the mass function does not intersect constraint curves, as experiments are typically sensitive to the integral of the mass function in each of a set of mass bins. Thus, most constraints are only trivial to interpret for monochromatic mass functions, i.e., mass functions of the form whose integrals are non-zero in only one bin. In this case, an observational constraint curve f max (M ) imposes the requirement that f 0 < f max (M 0 ). Transforming such constraints to the parameter space of a more general extended mass function involves summing contributions to observables from all mass bins. Multiple prescriptions for this procedure have been used in the literature.
The earliest systematic treatment of constraints for extended mass functions is due to [21]. They divide the mass range into N bins I 1 , . . . , I N , approximating the constraint functions as step functions on these bins. Within each bin, only the strongest constraint function f max (M ) is considered. A mass function ψ is excluded if for any k. This prescription is used by [16] to numerically transform observational constraints to the parameter space of a lognormal mass function. Their findings suggest that broadening the mass function does not generally relax constraints. However, as [16] treat the problem computationally, it is difficult to determine the relevance to their results of any particular constraint, or of the lognormal form of the trial mass functions. This provided partial motivation for the analysis of [17], who obtain similar numerical results for several additional constraints and forms of the mass function. Further, [17] derive a more rigorous prescription for transforming observational constraints to general extended mass functions. We review their derivation in section 2.2. Similar questions motivate the recent analysis of [18]. Rather than develop a prescription for translating constraints for monochromatic mass functions to suit a given extended mass function, the authors develop a prescription for converting the extended mass function into a set of monochromatic mass functions, each accounting for the contribution of the PBH population to one observable. The extended mass function is then subject to each constraint as it applies to the corresponding monochromatic mass function. This approach is used to constrain the parameter spaces of lognormal and power law mass functions, with results similar to those of [16] and [17].
Our methods bear some similarities to [18], in that we also find it sufficient to work with sets of monochromatic mass functions. However, the monochromatic mass functions we consider have a different interpretation, as discussed in section 3.3. Our goal is not to place constraints on any specific extended mass function, but rather to place bounds on PBH dark matter while allowing complete freedom in the mass function. Thus, our formalism is structured around the maximization problem, and we use our results to study both the current status of the PBH dark matter paradigm and the potential impact of future observables.

Constraint prescription
In this work, we seek a general result for the maximum allowed fraction of dark matter in PBH, independent of the form of the mass function, and in a form that elucidates the relevance of each observable. As such, it is necessary that we adopt a prescription for constraining a given mass function that allows for multiple simultaneous constraining observables, a requirement most naturally satisfied by that of [17]. Their prescription is thus the basis for our analytical work. We numerically confirm that similar results are obtained under the prescriptions of [18] and [21] (see appendix A.2).
We follow [17] to convert constraints for monochromatic mass functions to constraints for extended mass functions. We denote the mass function by ψ and adopt their normalization and conventions, such that where n is the number density of PBH at fixed mass. Most observables that can constrain primordial black holes are determined by the properties of single black holes, with no need to consider relationships between them. In such a case, an observable quantity A receives a linear combination of contributions from each mass bin, and the contribution from black holes of mass M is proportional to ψ(M ). As such, the observable can be written as a functional of ψ in the form We note in passing that there are some observables for which relationships between black holes are significant. For example, gravitational wave observations of mergers are dependent on the properties of pairs of black holes, and so one must combine contributions from pairs of mass bins. In the simplest case, where the contributions scale linearly with number in each mass bin, such an observable can clearly be written in the form and one can always express a generic observable by including higher-order terms of this form. Note that higher-order terms also account for non-linear dependence of A on ψ at fixed mass. For example, an observable which scales as ψ(M ) 2 can be expressed exactly at second order by setting We study the potential impact of gravitational wave observations in section 4. All of the other constraints that we consider in this work are of the simplest kind, and we will find eq. (2.4) sufficient. In this case, it is straightforward to relate constraints for a monochromatic mass function to constraints for a generic mass function, and we briefly review the argument given in [17].
is the largest coefficient allowed by constraints for a mass function of this form. If we take ψ(M ) = ψ mono (M 0 ; M ) in eq. (2.5), we obtain crosses a threshold value A obs . Then A[ψ mono ] = A obs by definition of f max , so eq. (2.6) gives K 1 (M ) independent of ψ. Substituting for K 1 (M ) in eq. (2.5) while leaving ψ generic gives the condition This expresses the constraint on a mass function ψ(M ) when the constraint for a monochromatic mass function is dM

Single-constraint case
For pedagogical purposes, we first consider the case of a single constraining observable. For such situations, when all observables can be expressed in the form of eq. (2.4), the constraint on the mass function has the form C[ψ] ≤ 1, with C[ψ] as defined in eq. (2.7). The problem is then to maximize dM ψ(M ) subject to this constraint. The optimal mass function saturates the constraint, so it suffices to require C[ψ] = 1.
Naively, this problem looks as though it can be solved using the method of Lagrange multipliers, by finding stationary points of the functional .
However, the Euler-Lagrange equation in ψ admits no non-trivial solutions. This is because dM ψ(M ) can be made arbitrarily large, even subject to C[ψ] = 1, unless ψ(M ) > 0 is imposed. Positivity can be imposed by setting ψ = φ * φ and performing an unconstrained optimization in φ, but the corresponding Euler-Lagrange equation leads to the condition that φ is, at every point, either zero or non-analytic.
The variational approach does not generalize to the case of multiple constraints, so we do not pursue it any further. Rather, we observe that since C[ψ] is linear, we have C C[ψ] −1 ψ = 1. Thus, we can impose C[ψ] = 1 by rescaling ψ → C[ψ] −1 ψ, and then the problem is to maximize the functional subject only to positivity. We call M[ψ] the normalized mass of ψ.
It is now simple to show that M[ψ] is maximized by taking ψ to be a monochromatic mass function. Let M max ≡ argmax f max (M ) and f mono ≡ f max (M max ), and define so that ψ 0 (M ) is the monochromatic mass function which maximizes the PBH density, and f mono is the maximum PBH density allowed for a monochromatic mass function. Choose any mass function ψ ≡ ψ 0 + δψ. Since ψ 0 vanishes everywhere except for M max , positivity of ψ requires that δψ(M ) ≥ 0 for all M = M max . Then we have .
, so no functional form allows a higher total PBH density than does the Dirac delta. In particular, for fixed PBH density, we conclude that an extended mass function is always more strongly constrained than the optimal monochromatic mass function. While this will not hold for the case of multiple constraints, it remains an excellent approximation if the constraints are weakest by far in a mass range where a single observable dominates.

Combining constraints
Realistically, the single-constraint case is too simplistic. In general, a mass function is ruled out on the basis of a χ 2 test statistic. If PBH are constrained by multiple observables A j , then the test statistic is found by adding the individual χ 2 statistics in quadrature. That is, To fail to reject ψ at some significance level requires that If we set N = 1, this reduces to is the analogue of f max (M ) for the jth constraint alone. For general N , [17] show that the constraint takes the form Since the individual constraints are added in quadrature, the argument applied to the singleconstraint case does not extend to the case of multiple constraints, and indeed, there are cases in which the density is not maximized by a monochromatic mass function. However, we will show that the maximizer is in general a linear combination of N monochromatic mass functions.

The general problem
For the case of several constraining observables, one has N constraint functions denoted by f max,1 , . . . , f max,N . For brevity, we define g j (M ) ≡ 1/f max,j (M ), and by analogy with eq. (2.7), we define Then the problem is to find ψ to maximize Since rescaling ψ does not change M[ψ], we can always set dM ψ(M ) = 1, and then the problem is equivalent to minimizing C[ψ] subject to this constraint. For convenience, we cast the integral in discrete form, writing where Q is not restricted to be finite. Thus, the problem is to minimize the norm of a sum of a k g(M k ) for some {M k } k=1,...,Q , subject to our normalization condition, which now takes the form Q k=1 a k = 1. Geometrically, this is the same as minimizing the norm over the convex hull of the g(M ), i.e., to compute where U is the mass range under consideration. We henceforth denote conv {g(M ) | M ∈ U } by conv(g). Since the minimizer is the projection of the origin onto a convex set, it is unique in the sense that any optimal mass function ψ must have the same C[ψ]. This does not require that the minimizing mass function is itself unique. Such a geometric formulation simplifies the interpretation of the problem. In particular, the result for the case of a single constraint is now immediate: the convex hull is 1-dimensional, so the point with minimum norm is simply the minimum value of g(M ). The corresponding mass function is monochromatic, with a peak at argmin g(M ). It is also clear that the monochromatic mass function is not generally the minimizer of the norm in the case of multiple constraints: we have no guarantee that g(M ) attains the minimum of the norm on conv(g) for any single M .
Still, minimizing the norm over the convex hull of a discretization of g(M ) is a simple computational problem, and it is easy to validate the result. We find an optimal mass function in three steps: We choose the M k using adaptive sampling to capture features of the constraint functions as precisely as possible. The convex hull of G is now a polytope A.
2. Find the point p min ∈ A with minimum norm. We implement the algorithm of [22], which requires only the extreme points of A as inputs. To avoid computing the convex hull in a high-dimensional space, we supply all of the points of G, of which the extreme points of A form a subset. The algorithm determines the facet S of A which contains p min , and gives the barycentric coordinates of p min in S as a vector w.

Define a mass function
Thus, ψ opt is a mass function which attains the maximum total dark matter fraction. In particular, for any mass function ψ, we have is an upper bound on the fraction of dark matter in PBH irrespective of the functional form of the mass function. We will refer to ψ opt as the semi-analytical optimum mass function.
We can now explain geometrically why the maximizing mass function is a linear combination of no more than N monochromatic mass functions. Observe that for any g(M ), the minimum of the norm must lie on the boundary of the convex hull conv(g), and since g(M k ) ∈ R N , this boundary has dimension at most N − 1. One can construct an arbitrarily refined triangulation of this boundary formed from (N − 1)-simplices, each with N points of G as vertices. The minimizer of the norm is a linear combination of these vertices, each of which is one of the original g(M k ), corresponding to a monochromatic mass function. We emphasize that at no step do we impose that the optimal mass function is a discrete linear combination of a finite number of monochromatic mass functions. This is a consequence of the fact that the optimum corresponds to a point in an (N − 1)-simplex, meaning that this mass function lies in a space which is spanned by at most N monochromatic mass functions.
Our method is deceptively similar to the procedure of [18], in that we also work with sets of monochromatic mass functions. However, the monochromatic mass functions considered in that work are used only to study the consequences of a given extended mass function for a given observable. The sum of the effective monochromatic mass functions corresponding to each observable does not generally give a single mass function with equivalent consequences for all observables combined. This is appropriate for the purposes of [18] because they investigate which constraints are most effective for mass functions of a fixed functional form.
In this work, the mass functions we derive maximize the density of PBH with respect to all constraints simultaneously. Since the constraints are statistical in nature, the combination of multiple independent constraints at a single mass is stronger than any one of them individually. We follow [17] in treating constraints simultaneously, and our resulting semi-analytical optima are indeed sums of monochromatic mass functions. This approach is necessary for our purposes because we investigate which constraints and mass ranges are most significant for overall constraints on PBH dark matter, irrespective of the functional form of the mass function.

Results
We perform the maximization explicitly for several sets of constraints. Set A includes robust constraints from evaporation [23]; GRB lensing [24]; microlensing from HSC [10], Kepler [11], EROS [12], and MACHO [13]; and CMB limits from Planck [14]. Set B includes dynamical constraints from Segue I [25], Eridanus II [26], and non-disruption of wide binaries [15]. Set C includes a constraint from white dwarf explosions [27], a constraint from neutron star capture [28] and a recently claimed constraint from SNe lensing in the LIGO window [29]. The constraints from evaporation and from Planck in A have been estimated differently in the literature, with important consequences for our analysis. Set A itself contains relatively non-restrictive estimates of these constraints. We incorporate more stringent versions (see section 5) of these constraints in a setĀ, which is otherwise identical to A.
We determine optimal mass functions for sets A,Ā, and all of their combinations with sets B and C. The results are summarized in table 1 and illustrated in fig. 1. We do not include cosmological constraints on the total matter density, so these values of f max,all may exceed 1. In particular, note that all combinations containing A have f max,all > 1, while all combinations containingĀ and C have f max,all < 1. The setĀ on its own has marginal status if only monochromatic mass functions are considered, but clearly f max,all > 1 in this case. With the constraints we consider in this work, f PBH = 1 is always allowed when using the less stringent set A, regardless of additional constraints.

Prospects for gravitational wave constraints
Gravitational wave observables are the major exception to the rule that measured quantities are linear in the PBH mass function in each mass bin. There are several methods by which gravitational waves might constrain the primordial black hole population. In principle, the simplest constraint arises from present-day measurements of the black hole binary (BHB)  Vertical lines denote the locations of Dirac deltas in the semi-analytical optimum mass function, with height indicating the weight given to each one. The labeled constraints are from BH evaporation (evap, [23]), GRB femtolensing observations (FL, [24]), white dwarf explosions (WD, [27]), Hyper Suprime-Cam (HSC, [10]), Kepler (K, [11]), EROS-II (EROS, [12]), supernova lensing (SNe, [29]), MACHO (MACHO, [13]), Segue I dynamics (SegI, [25]), Eridanus II dynamics (EriII, [26]), wide binary dynamics (WB, [15]), and CMB observables (CMB, [14,17]). merger rate, but this is weak for two reasons: first, it is difficult to distinguish primordial black holes from astrophysical black holes, for which a variety of additional physical mechanisms may affect the merger rate; and second, the observed merger rate is sufficiently uncertain as to be compatible with a wide range of PBH dark matter models [7,9].
An alternative method is to search for the stochastic gravitational wave background from primordial density fluctuations associated with inflationary production mechanisms [30]. However, such constraints are only effective within the context of this class of formation models, and within such a limited scope, our level of generality is excessive. Here it is sufficient to consider the mass functions that can be reasonably produced by such formation mechanisms, and constraints for such mass functions have been treated elsewhere in the literature.
A third technique is to search for the stochastic gravitational wave background due to BHB mergers throughout cosmic history [31]. While such an approach may ultimately produce strong constraints, there remains a great deal of uncertainty in modeling the merger rate, particularly for extended mass functions. This problem has only recently been treated in the literature [32,33], and the resulting constraints may not be robust. Still, it is useful to estimate these constraints, even imprecisely, in order to determine their relevance in the case of our semi-analytical optimum mass functions.
We now consider constraints from the non-detection of a stochastic background of gravitational waves from BHBs throughout cosmic history. This background is qualitatively different from all other observables considered in this work, since it has complicated non-linear dependence on the mass function. This means that determining constraints on a general mass function is non-trivial. In particular, one needs to include the higher order terms in the expansion of eq. (2.5), and the analogue of eq. (2.6) is then which does not constrain off-diagonal values of the kernels K n . Thus, gravitational wave constraints on the parameter space of monochromatic mass functions are insufficient to determine constraints on an extended mass function, even when the functional form is specified. This reflects the fact that gravitational wave constraints on extended mass functions are inherently model-dependent, in that one must determine the contribution to the background from binaries whose partners have unequal masses. The results of [32] provide a simple method for estimating the stochastic gravitational wave background given a particular mass function, which we now review. The observable characteristic strain amplitude is given by is a cosmology-dependent constant, ν GW is the gravitational wave frequency, τ merge is the mass-dependent binary merger rate per unit volume, and M c ≡ (m 1 m 2 ) 3/5 (m 1 + m 2 ) −1/5 is the chirp mass. The merger rate is determined by consideration of the capture rate for formation of PBH binaries, given in the Newtonian approximation by where n PBH is the local number density of PBH, v PBH is the characteristic velocity of a single black hole, and v rel ≡ √ 2v PBH is the characteristic relative velocity of two black holes. The local number density of PBH of mass m is parametrized as n PBH = δ PBH ρ PBH (m)/m, where ρ PBH (m) is the cosmological average density of PBH of mass m and δ PBH is the local density contrast of PBH. In particular, in our notation, this number density is given by Estimates for δ PBH range from 10 6 to 10 10 . In order to estimate conservative constraints, we henceforth take the relatively low value δ PBH = 10 7 .
In terms of the capture rate, the merger rate per unit volume is τ merge (m 1 , m 2 ) = (f PBH /δ PBH )τ capture (m 1 , m 2 )n PBH (m 2 ). Thus the strain amplitude is given by where C is a cosmology-dependent factor given by In particular, the dependence of h 2 c (ν GW ) on ψ admits a simple expansion of the form of eq. (2.5). If f PBH is fixed independently of ψ, then the only non-vanishing term has the form dm 1 dm 2 ψ(m 1 )ψ(m 2 )K 2 (m 1 , m 2 ), with Note that K 2 varies with ν GW , reflecting the fact that measurements of h c (ν GW ) at each frequency ν GW are independent. Thus, if the instrument used has sufficiently high frequency resolution, a large number of independent constraining observables can be measured, meaning that the maximizing mass function need not resemble a linear combination of a small number of monochromatic mass functions. In this case, direct numerical methods are necessary to determine the maximum density of PBH. It would be inappropriate to perform a full numerical optimization within our framework, since there are considerable theoretical uncertainties in the determination of merger rates. However, we can estimate the potential impact of gravitational wave constraints by checking compatibility of our optimal mass functions with existing gravitational wave observations. This serves to indicate the potential for future modeling work to constrain the PBH population: if the semi-analytical optimum for a given set of constraints is already consistent with gravitational wave observations as well, then we can predict that the detailed inclusion of this additional constraint will have a minimal impact on f max,all .
Current aLIGO bounds on h c (ν GW ) are strongest at ν GW 100 Hz, and we represent the current limit by h c (100 Hz) 10 −22 . Given a functional form for the mass function ψ, we can compute the maximum f PBH for which h c [ψ](100 Hz) satisfies this bound. We denote this maximum by f max,GW , and the values of f max,GW for our semi-analytical optima are shown in table 1. Of the sets of constraints we consider, only sets AC andĀC have f max,GW significantly less than one. This is to be expected, since the maximizing mass function in both of these cases has a large variance: [32] show that the gravitational wave background is strongly enhanced as the variance is increased. In the case of setĀC, we have f max,GW f max,all < 1, meaning that the overall maximum is minimally impacted by gravitational wave constraints. In particular, f PBH = 1 is ruled out regardless.
Only set AC has f max,GW < 1 < f max,all , which, in general, is difficult to interpret: in principle, there may be a different form for the mass function which relaxes gravitational wave constraints, retains f PBH ≥ 1, and remains consistent with the other constraints we consider in this work. A simple way to check this is to consider the maximizing monochromatic mass function, for which gravitational wave constraints should be relaxed compared with the high-variance semi-analytical optimum. Indeed, there is a monochromatic mass function which satisfies all non-gravitational constraints in set AC with f PBH = 1.371, and we compute f max,GW > 10 for this mass function. Thus, while it may not be possible to attain f max,all in this case without violating gravitational wave constraints, f PBH = 1 clearly remains allowed. As such, the addition of gravitational wave constraints does not change the overall status of the PBH dark matter paradigm for any of the constraint sets we consider. This reflects both the current status of observations and the large uncertainties in modeling the background. However, future experiments are expected to improve limits on h c (ν GW ) by 2-4 orders of magnitude, which would be sufficient to rule out all of the mass functions represented in table 1 even under fairly conservative assumptions.

Discussion
With the maximization procedure introduced in section 3.3, it is simple to determine the maximum PBH density consistent with constraints. We stress that this is a bound that applies for mass functions of all forms. Thus, given a set of observational constraints, we can determine a model-independent bound on the density of PBH.
Our results quantify, for the first time, the risks of using monochromatic mass functions to assess the overall status of the PBH dark matter paradigm. So long as one window in the constraint functions is much less constrained than all others, the difference between f max,all and f mono is generally very small. Set A is a clear example of such a case, and the correction is of order 0.1%. On the other hand, if PBH are constrained to a similar extent in multiple windows, the correction can be large. The most dramatic example is provided by setĀ, for which f max,all is larger than f mono by ∼ 50%. We conclude that, at worst, the bound on the total PBH density is related to the monochromatic bound by an O(1) factor.
The optimal mass functions themselves ( fig. 1) do not correspond to any well-motivated production scenario that we are aware of, and we certainly do not claim that the maximal density can be attained by producing PBH monochromatically at a discrete collection of masses spanning 15 orders of magnitude. Instead, the panels of fig. 1 should be interpreted as a tool to relate monochromatic constraint functions to their impact on the allowed total density of PBH. In particular, an immediate and non-trivial conclusion that can be drawn from the figures is that the addition of any new constraint which does not overlap the peaks of the optimal mass function will not reduce f max,all .
Further, the functional form of the optimal mass function clarifies the dependence of constraints on the variance of the mass function. In the single-constraint case, we showed that an extended mass function never outperforms the optimal monochromatic mass function. Indeed, in this case, increasing the variance of a narrow mass function will only relax constraints if f max is concave-up in the mass range of interest, i.e., if the monochromatic mass function under consideration is not the optimal one. When multiple constraints are considered, the relationship between the variance of the mass function and the allowed density is less obvious. Our semi-analytical optimum mass functions all exhibit some non-zero spread, and they definitively allow higher PBH densities than any zero-variance (i.e., monochromatic) mass function. However, extending a monochromatic mass function only slightly, without overlapping additional points of the semi-analytical optimum mass function, is not useful for relaxing constraints. In this respect, our findings are consistent with those of [17,18].
The most substantial differences in f max,all arise from differences between A andĀ. Set A contains more stringent forms of constraints from CMB anisotropy and PBH evaporation. The CMB constraint is strongly dependent on modeling poorly-understood accretion processes. Both versions of the constraint used in this work are drawn from [14]: the version in set A is obtained by considering only collisional ionization of the accreted gas, while the version in setĀ is obtained by including photoionization as well. The evaporation constraint is sensitive to uncertainties in the spectrum of extragalactic background radiation. We adopt the extreme cases considered by [17], with the relaxed form contained in set A and the more stringent form in setĀ.

Relative impact of constraints
The values of f max,all in table 1 demonstrate that the present observational status of PBH dark matter is strongly dependent on the constraints adopted. However, to rule out f PBH = 1, it is necessary to both take the more stringent constraintsĀ in place of A, and to include at least one of the constraints from set C: supernova microlensing [29], neutron star capture [28], and white dwarf explosions [27].
The supernova microlensing constraint is the most recent of those we consider, and its robustness is the subject of ongoing discussion in the literature [see e.g. 9]. We note that this constraint is dominant in the LIGO window only when dynamical constraints from set B are neglected, so the addition of this constraint alone to set AB orĀB will have a small impact on f max,all . The constraint from neutron star capture is also subject to astrophysical uncertainties, since it is dependent on the dark matter density in the cores of galactic clusters [28]. We consider the relatively restrictive constraint obtained by taking ρ DM = 10 4 GeVcm −3 . The strength of the constraint scales linearly with ρ DM , and more conservative estimates take ρ DM smaller by an order of magnitude or more. However, this constraint is most effective in a window shared with constraints from white dwarf explosions, so even if one of the two is subject to substantial uncertainties, the effect of set C on f max,all remains large.
The form of the optimal mass function allows us to rapidly identify the potential impacts of prospective constraints from future observations. For instance, constraints on intermediatemass black holes with M 10 2 M are already strong enough that our semi-analytical optimal mass functions are negligibly small throughout this region. Thus, the identification of additional dynamical systems that might tighten constraints in this region will not affect the overall bound on the PBH density at a level greater than one part in 10 4 . On the other hand, GRB femtolensing limits lie in a mass range where some of the semi-analytical optima have a large peak, and strengthening these constraints will have an immediate impact on the overall bound. In particular, upcoming Fermi GRB observations are expected to substantially strengthen constraints in this window, improving by a factor of five after 10 years of operation [34]. These results may ultimately rule out the PBH dark matter paradigm, with the exception of non-evaporating Planck-mass relics.
Constraints from gravitational wave observations are a special case, as they do not admit the linear interpretation that we take for the impact of extended mass functions on other observables. Further, the strain amplitude at each frequency is sensitive to PBH in a wide range of masses. Thus, it is not trivial to predict the effect of future gravitational wave constraints on our overall bound without direct numerical optimization. However, constraints from LISA and DECIGO [see e.g. 35] will eventually be capable of ruling out all of our semi-analytical optima, potentially lowering the upper bounds we set in this work.

Conclusions
We have found the form of the mass function which maximizes the PBH density subject to observational constraints, and we have used this to calculate an upper bound on the fraction of dark matter in PBH. Depending on the constraints adopted, we find f max,all as large as 27.25 (set A) or as small as 0.405 (setĀBC). The scenario in which all dark matter is composed of PBH is ruled out by stringent limits from evaporation and Planck if combined with the constraints from white dwarf explosions, neutron star capture and SNe lensing (set C). However, if relaxed constraints from evaporation and Planck are adopted, PBH dark matter is not ruled out by the addition of any other constraints we consider in this work. Estimated gravitational wave constraints do not affect these conclusions at the sensitivity of current instruments.
Our method provides a fast and robust technique to determine the total allowed density of PBH given a set of constraints (f max,all ), independent of the form of the PBH mass function. The optimal mass function itself allows an easy test of the impact of additional constraints on f max,all . While the optimal mass function is not exactly monochromatic, it is very nearly so for realistic constraints. The optimal mass function corresponding to each set of constraints we consider is approximately monochromatic, with additional components scaling the total allowed fraction by no more than an O(1) factor. Our results explain the findings of [17,18] that extended mass functions are generally more strongly constrained than monochromatic mass functions, and confirm that the monochromatic maximum density f mono is a good approximation of the allowed density across all mass functions.

A Numerical validation
Given a set of constraints, it is also possible to use numerical methods to find a mass function which maximizes the PBH density. There are significant caveats to such an approach. Most importantly, a maximization algorithm may converge to a local optimum rather than a global optimum. Additionally, computational costs may render numerical approaches impractical unless the functions involved are discretized sparsely. Even so, numerical optimization can be used to validate our analytical results: if the same set of masses is used for discretization, then the numerical result should never reach a greater normalized mass (cf. eq. (3.12)) than that of our corresponding semi-analytical result. Numerical methods can also be used to check that our semi-analytical optimum is a stationary point of the normalized mass functional.

A.1 Direct validation
We implement these validation steps using a simple Monte Carlo algorithm, as follows: we begin with an initial mass function of the form ψ 0 (M ) ∝ M −1 , which assigns equal PBH density to each log-spaced mass bin. We then perturb the value of ψ 0 in a random bin k by a value selected from a Gaussian distribution with mean 0 and variance σ 2 ψ 0 (M k ) 2 , where σ is a parameter of the maximization. We denote the resulting mass function by ψ 1 (M ). If ψ 1 (M k ) ≥ 0 and M[ψ 1 ] > M[ψ 0 ], we accept the step, replace ψ 0 by ψ 1 , and repeat. For simplicity, we do not accept any steps which reduce the normalized mass. This is not necessary in order to test whether our semi-analytical optimum mass function is a stationary point. We also reject steps which increase the normalized mass by less than 10 −10 to avoid exceeding the numerical precision of the semi-analytical result.
In order to make the problem numerically tractable, we use only 10 2 log-spaced mass bins. This discretization is different from the one used in table 1, and it does not capture sharp features of the constraints. Consequently, in order to compare the numerical results with semi-analytical results, we regenerate the semi-analytical mass function with the same discretization. Note that this affects both the form of the optimal mass function and the calculated f max,all .
We implement the numerical optimization with σ = 10 −2 . In what follows, we denote the numerical mass function by ψ N , and the semi-analytical optimum by ψ SA . The left-hand side of fig. 2 shows f PBH for the numerical mass function at each step as a fraction of the semi-analytical f max,all . The numerical f PBH converges to f max,all and immediately stabilizes, and in particular, in no step does f PBH exceed f max,all .
In principle, ψ N need not converge to ψ SA even given that f PBH converges to f max,all , since the mass function with maximal density is not necessarily unique. However, in the top-right panel of fig. 2, we show that ψ N tends to ψ SA in the L 2 norm. To compute this distance consistently, we treat the Dirac deltas of ψ SA as constant functions in their respective bins. As an additional test of convergence, we compute the acceptance rate, i.e., the fraction of steps which are accepted, during each window of 10 4 iterations. The acceptance rate vanishes as ψ N approaches ψ SA , which further demonstrates that ψ SA is a stationary point of the normalized mass.
The numerical and semi-analytical mass functions are shown in fig. 3. In order to compare Dirac deltas with the smooth mass function ψ N , the figure shows the integral of the mass function in each bin rather than ψ N and ψ SA themselves. It is clear that in this case, the numerical algorithm converges to the semi-analytical optimum. We have established via analytical arguments that this is not simply a local optimum, but indeed the global maximum of the normalized mass.

A.2 Sensitivity to the constraint prescription
Our analytical work is based on the prescription of [17] for evaluating constraints on extended mass functions. Since other prescriptions have been considered in the literature, it is important to determine the robustness of our results to variations on the constraint prescription. Assessing this analytically is intractable, as it requires the development of independent analytical frameworks for even slight modifications. However, numerical methods allow for a comparison of the bounds we obtain analytically with those that would result from any other specified prescription. We thus perform numerical optimization under the prescriptions of [18] and [21], and compare these with our semi-analytical results.
The situation is particularly simple for the constraint prescription of [18]: the major difference is that a mass function is allowed if it is allowed according to each individual constraint, rather than according to their statistical combination. Thus, the normalized   . For the case shown in fig. 3, these two numerical maxima agree to within 1%. We also implement the constraint procedure of [21], for which the constraints are treated as step functions on a set of mass bins. There is no universal prescription for the size of the bins across all constraints, but they should be chosen small enough that the minimum of each constraint function is not very different from its maximum within any single bin. Altogether, the procedure is as follows: 1. The mass range is divided into bins I 1 , . . . , I n .
2. In each bin, the dominant constraint is identified. If a bin captures the transition between two dominant constraints, the bin is subdivided at the transition point.
3. We evaluate the constraint of eq. (2.2), considering only the dominant constraint in each bin. We treat ψ as a smooth function, using a more refined set of bins for its numerical representation.
As in appendix A.1, we take a relatively coarse binning of the mass range for numerical testing purposes. We have evaluated the maxima attained for a range of different constraint bin counts, and we find that our determination of f max,all is robust to changes in binning at the 10% level. We note as well that prescriptions of this kind have been criticized in the literature for the fact that it is not trivial to determine the range of validity of the individual constraints, and hence to determine the limits of integration for each constraint curve in eq. (3.10). This is indeed a concern when individual constraints are considered, and it introduces significant potential uncertainty in cases where there are masses for which all constraints are at the edge of their range of validity. In these scenarios-for instance, with the constraints of set A-there are effectively "windows" in the constraints that make any upper bound on the density of PBH quite uncertain. However, in the other cases that we consider, the dominating constraints in a given mass range generally intersect well within their respective ranges of validity. Thus, modifying the limits of integration does not substantially impact our results.