Primordial black holes from scalar field evolution in the early universe

TL;DR

This paper explores how scalar field dynamics in the early universe, including condensate fragmentation and topological defects, can produce primordial black holes across a wide mass spectrum, potentially explaining LIGO observations.

Contribution

It introduces a mechanism for primordial black hole formation from scalar field evolution, extending the mass range beyond previous limits, relevant for gravitational wave detections.

Findings

Scalar condensates can fragment into Q-balls, leading to PBH formation.

The mechanism can produce black holes from 1 to 100 solar masses.

Topological defects can also generate PBHs through similar processes.

Abstract

Scalar condensates with large expectation values can form in the early universe, for example, in theories with supersymmetry. The condensate can undergo fragmentation into Q-balls before decaying. If the Q-balls dominate the energy density for some period of time, statistical fluctuations in their number density can lead to formation of primordial black holes (PBH). In the case of supersymmetry the mass range is limited from above by $10^{23}$g. For a general charged scalar field, this robust mechanism can generate black holes over a much broader mass range, including the black holes with masses of 1-100 solar masses, which is relevant for LIGO observations of gravitational waves. Topological defects can lead to formation of PBH in a similar fashion.

Equations

$$\omega \sim \Lambda Q^{\alpha -1},R\sim |Q{|}^{\beta }/\Lambda ,$$

$${\varphi }_0\sim \Lambda |Q{|}^{1-\alpha },M\sim \Lambda |Q{|}^{\alpha },$$

$$f_M(M)=\frac{M^{\frac{1-\alpha }{\alpha }}}{\alpha {\Lambda }^{1/\alpha }}\left[f_Q((M/\Lambda {)}^{1/\alpha })+f_Q(-(M/\Lambda {)}^{1/\alpha })\right].$$

$$f_M(M,Q_{\text{tot}},N)=\frac{\psi (M,Q_{\text{tot}},N)}{\int dM\psi (M,Q_{\text{tot}},N)},$$

$$\begin{array}{rl}\psi (M,Q_{\text{tot}},N)& =\int \left(\prod _{i=1}^{N}d{Q}_i{f}_Q(Q_i)\right)\\ & \times \delta \left(M-\Lambda \sum _{i=1}^N|Q_i{|}^{\alpha }\right)\delta \left(Q_{\text{tot}}-\sum _{i=1}^N{Q}_i\right),\end{array}$$

$$\tilde{\psi }({\xi }_M,{\xi }_{Q_{\text{tot}}},N)={\left[\int dQ{e}^{i({\xi }_M|Q{|}^{\alpha }+{\xi }_{Q_{\text{tot}}Q})}{f}_Q(Q)\right]}^N.$$

$$p(N)=e^{-N_f}\frac{N_f^N}{N!},$$

$$⟨M⟩=\sum _{N=1}^{\infty }{\int }_0^{\infty }dMM{F}_Q(M,Q_f,N)\approx M_f,$$

$$⟨N⟩=\sum _{N=1}^{\infty }{\int }_0^{\infty }dMN{F}_Q(M,Q_f,N)=N_f,$$

$$N_{\text{RMS}}={\left[\sum _{N=1}^{\infty }{\int }_0^{\infty }dM{N}^2{F}_Q(M,Q_f,N)\right]}^{1/2}=N_f,$$

$$p(N,V)=e^{-N_{f}V/V_f}\frac{(N_{f}V/V_f{)}^N}{N!}.$$

$$F_Q(M,V,N)=\delta \left(M-M_f{\left(\frac{N}{N_f}\right)}^{1-\alpha }{\left(\frac{V}{V_f}\right)}^{\alpha }\right)p(N,V).$$

$$\sum _{\{V\}}g(V)$$

$$\approx {\int }_1^{1+{\log }_{\chi }\frac{V_{\text{max}}}{V_{\text{min}}}}dig(V_{\text{max}}/{\chi }^{i-1})$$

$$=\frac{1}{\ln \chi }{\int }_{V_{\text{min}}}^{V_{\text{max}}}\frac{dV}{V}g(V),$$

$$⟨{\rho }_Q(t_f)⟩=\underset{V\to \infty }{\lim }\frac{⟨M⟩}{V}=\frac{M_f}{V_f},$$

$$\delta (t_f)=\frac{\delta {\rho }_Q}{⟨{\rho }_Q⟩}=\frac{M/V}{⟨{\rho }_Q⟩}-1={\left(\frac{N/N_f}{M/M_f}\right)}^{\frac{1-\alpha }{\alpha }}-1$$

$$⟨{\rho }_{\text{BH}}(t_R)⟩={\left(\frac{a_f}{a_R}\right)}^3\sum _{N=1}^{\infty }{\int }_{V_{\text{min}}}^{V_R}\frac{dV}{V}{\int }_0^{\infty }dM\left(\beta \frac{M}{V}\right)F_Q,$$

$$\frac{d⟨{\rho }_{\text{BH}}⟩}{dM}={\left(\frac{a_f}{a_R}\right)}^3\sum _{N=1}^{\infty }{\int }_{V_{\text{min}}}^{V_R}\frac{dV}{V}\left(\beta \frac{M}{V}\right)F_Q.$$

$${\rho }_R(t_{\text{RH}})=\frac{{\pi }^2}{30}{g}_{\ast }(T_{\text{RH}})T_{\text{RH}}^4\approx \frac{{\pi }^2}{327}{\Gamma }_I^2{M}_p^2,$$

$${\rho }_R(t)=\frac{{\pi }^2{M}_p^2}{327t^2};t_{\text{RH}}<t<t_Q$$

$$⟨{\rho }_Q(t)⟩$$

$$=\frac{3\Lambda Q_f^{\alpha }{N}_f^{1-\alpha }}{4\pi t_f^{3/2}{t}^{3/2}}{e}^{-(t-t_f)/{\tau }_Q};t_f<t<t_Q$$

$${\rho }_R(t)=\left[{\rho }_R(t_Q)+⟨{\rho }_Q(t_Q)⟩{\int }_{x_0}^{x}d{x}^{\prime }z(x^{\prime })e^{-x^{\prime }}\right]z^{-4},$$

$$⟨{\rho }_Q(t)⟩=\frac{3\Lambda Q_f^{\alpha }{N}_f^{1-\alpha }{t}_Q^{1/2}}{4\pi t_f^{3/2}{t}^2}{e}^{-(t-t_f)/{\tau }_Q};t_Q<t<t_R$$

$$1+{\left(\frac{t_R}{{\tau }_Q}\right)}^{-2/3}\Gamma \left(\frac{5}{3},\frac{t_Q}{{\tau }_Q},\frac{t_R}{{\tau }_Q}\right)={\left(\frac{t_R}{t_Q}\right)}^{2/3}{e}^{(t_Q-t_R)/{\tau }_Q},$$

$$t_Q=t_f{r}_f,t_R=t_f{r}_{f}r,{\tau }_Q=t_f{r}_f{r}_Q(r),$$

$$\Lambda Q_f^{\alpha }=\frac{4\pi M_p^2{t}_f}{3\cdot 327r_f^{1/2}{N}_f^{1-\alpha }}{e}^{(1-1/r_f)/r_Q(r)},$$

$$f_i=\frac{1}{{\rho }_{\text{DM}}}{\int }_{M_i}^{M_{i+1}}dM\frac{d⟨{\rho }_{\text{BH}}(t_0)⟩}{dM},$$

$$\tau \approx \frac{1-Q_0^{1-a}}{(a-1){\Gamma }_0},$$

$$a=\frac{1}{4}(7+2(n+m-2)),$$

$${\Gamma }_0=112.7|g{|}^2{e}^{-0.236(n+m)}(n-m{)}^2{J}_{nm}\Lambda (\Lambda /{\Lambda }_{\ast }{)}^{2(n+m)-8},$$

$${\rho }_R(t)=\{\begin{array}{ll}\frac{{\pi }^2{M}_p^2}{327t^2}& t_{\text{RH}}<t<t_Q\\ \frac{{\pi }^2{M}_p^2}{327t_Q^2}\left[1+{\left(\frac{{\tau }_Q}{t}\right)}^{2/3}\Gamma \left(\frac{5}{3},\frac{t_Q}{{\tau }_Q},\frac{t}{{\tau }_Q}\right)\right]{\left(\frac{t_Q}{t}\right)}^{8/3}& t_Q<t<t_R\\ \frac{{\pi }^2{M}_p^2}{327t_Q^2}\left[1+{\left(\frac{{\tau }_Q}{t_R}\right)}^{2/3}\Gamma \left(\frac{5}{3},\frac{t_Q}{{\tau }_Q},\frac{t_R}{{\tau }_Q}\right)\right]{\left(\frac{t_Q}{t_R}\right)}^{8/3}{\left(\frac{t_R}{t}\right)}^2& t_R<t<t_{\text{eq}}\\ \frac{{\pi }^2{M}_p^2}{327t_Q^2}\left[1+{\left(\frac{{\tau }_Q}{t_R}\right)}^{2/3}\Gamma \left(\frac{5}{3},\frac{t_Q}{{\tau }_Q},\frac{t_R}{{\tau }_Q}\right)\right]{\left(\frac{t_Q}{t_R}\right)}^{8/3}{\left(\frac{t_R}{t_{\text{eq}}}\right)}^2{\left(\frac{t_{\text{eq}}}{t}\right)}^{8/3}& t_{\text{eq}}<t<t_0\end{array}$$

$$⟨{\rho }_Q(t)⟩=\{\begin{array}{ll}\frac{3M_f}{4\pi t_f^3}{\left(\frac{t_f}{t}\right)}^{3/2}{e}^{-(t-t_f)/{\tau }_Q}& t_f<t<t_Q\\ \frac{3M_f}{4\pi t_f^3}{\left(\frac{t_f}{t_Q}\right)}^{3/2}{\left(\frac{t_Q}{t}\right)}^2{e}^{-(t-t_f)/{\tau }_Q}& t_Q<t<t_R\\ \frac{3M_f}{4\pi t_f^3}{\left(\frac{t_f}{t_Q}\right)}^{3/2}{\left(\frac{t_Q}{t_R}\right)}^2{\left(\frac{t_R}{t}\right)}^{3/2}{e}^{-(t-t_f)/{\tau }_Q}& t_R<t<t_{\text{eq}}\\ \frac{3M_f}{4\pi t_f^3}{\left(\frac{t_f}{t_Q}\right)}^{3/2}{\left(\frac{t_Q}{t_R}\right)}^2{\left(\frac{t_R}{t_{\text{eq}}}\right)}^{3/2}{\left(\frac{t_{\text{eq}}}{t}\right)}^2{e}^{-(t-t_f)/{\tau }_Q}& t_{\text{eq}}<t<t_0\end{array}$$

$$⟨{\rho }_{\text{BH}}(t_R)⟩$$

Full text

††preprint: IPMU16-0192

I Introduction

It is well established that stellar core collapse can lead to formation of black holes. However, it remains an open question whether some processes in the early universe could produce primordial black holes (PBH) Zel’dovich and Novikov (1967); Hawking (1971); Carr and Hawking (1974); Garcia-Bellido et al. (1996); Khlopov (2010); Frampton et al. (2010); Kawasaki et al. (2012, 2016); Carr et al. (2016a); Inomata et al. (2016, 2017); Georg and Watson (2017); Carr et al. (2017); Domcke et al. (2017). PBHs can account for all or part of dark matter Zel’dovich and Novikov (1967); Hawking (1971); Carr and Hawking (1974); Garcia-Bellido et al. (1996); Khlopov (2010); Frampton et al. (2010); Kawasaki et al. (2016); Carr et al. (2016a); Inomata et al. (2016, 2017); Georg and Watson (2017); Carr et al. (2017); Domcke et al. (2017). Furthermore, they could be responsible for some of the gravitational wave signals observed by LIGO Bird et al. (2016), In addition, PBHs can invade and destroy neutron stars, ejecting neutron rich material in the process, which can account for for all or part of the r-process nucleosynthesis, as well as the 511-keV line in the galactic center Fuller et al. (2017). Finally, PBHs could provide seeds for supermassive black holes Kawasaki et al. (2012). A number of scenarios for black hole formation have been considered Khlopov (2010), and many of them rely on a spectrum of primordial density perturbations that has some additional power on certain length scales, which can be accomplished by means of tuning an inflaton potential.

It was recently pointed out that PBHs can form in a very generic scenario, which does not require any particular spectrum of density perturbations from inflation Cotner and Kusenko (2016a). Scalar fields with slowly growing potentials form a coherent condensate at the end of inflation Bunch and Davies (1978); Linde (1982); Affleck and Dine (1985); Starobinsky and Yokoyama (1994). In general, the condensate is not stable, and it breaks up in lumps, which evolve into Q-balls Kusenko and Shaposhnikov (1998). The gas of Q-balls contains a relatively low number of lumps per horizon, and the mass contained in these lumps fluctuates significantly from place to place. This creates relatively large fluctuations of mass density in Q-balls across both subhorizon and superhorizon distances. Since the energy density of a gas of Q-balls redshifts as mass, it can come to dominate the energy density temporarily, until the Q-balls decay, returning the universe to a radiation dominated era. The growth of structure during the Q-ball dominated phase can lead to copious production of primordial black holes. In this paper we will investigate this scenario in further detail.

Formation of Q-balls requires nothing more than some scalar field with a relatively flat potential at the end of inflation. For example, supersymmetric theories predict the existence of scalar fields with flat potentials. PBH formation in supersymmetric theories is, therefore, likely, even if the scale of supersymmetry breaking exceeds the reach of existing colliders.

A similar process can occur with topological defects, which can also lead to relatively large inhomogeneities. The discussion of topological defects is complicated by their non-trivial evolution. We will focus primarily on Q-balls, and will briefly comment on topological defects.

The format of this paper is as follows: in Section II, we describe the fragmentation of the condensate and the production of Q-balls, then in Section III we derive the formalism for calculating the statistical moments of collections of Q-balls. In Section IV we use the results of the previous Section to calculate the expected PBH density and mass spectrum, and in Section V we account for the effects on cosmological thermal history and evolve the PBH distribution to the present day. In Section VI, we then compare our results with current observational constraints, and in Section VII explore the available parameter space. In Section VIII, we comment on the applicability of this mechanism to topological defects.

II Formation of Q-balls

Formation of a scalar condensate after inflation and its fragmentation Kusenko and Shaposhnikov (1998) is a fairly generic phenomenon. While supersymmetry is a well-motivated theory for scalar fields carrying global charges and having flat potentials Affleck and Dine (1985); Dine and Kusenko (2003), our discussion can be easily generalized to an arbitrary scalar field with a global U(1) symmetry in the potential. Supersymmetric potentials generically contain flat directions that are lifted only by supersymmetry breaking terms. Some of the scalar fields that parameterize the flat directions carry a conserved U(1) quantum number, such as the baryon or lepton number. During inflation, these field develop a large vacuum expectation value (VEV) Bunch and Davies (1978); Linde (1982); Affleck and Dine (1985); Starobinsky and Yokoyama (1994), leading to a large, nonzero global charge density. When inflation is over, the scalar condensate $\phi(t)=\phi_{0}(t)\exp\{i\theta(t)\}$ relaxes to the minimum of the potential by a coherent classical motion with $\dot{\theta}\neq 0$ due to the initial conditions and possible CP violation at a high scale.

The initially homogeneous condensate is unstable with respect to fragmentation into non-topological solitons, Q-balls Coleman (1985). Q-balls exist in the spectrum of every supersymmetric generalization of the Standard Model Kusenko (1997a, b), and they can be stable or long-lived along a flat direction Kusenko and Shaposhnikov (1998); Dvali et al. (1998). In the case of a relatively large charge density (which is necessary for Affleck-Dine baryogenesis Affleck and Dine (1985); Dine and Kusenko (2003)), the stability of Q-balls can be analyzed analytically Kusenko and Shaposhnikov (1998); Enqvist and McDonald (1998a, 1999); these results agree well with numerical simulations Kasuya and Kawasaki (2000a). One finds that the almost homogeneous condensate develops an instability with wavenumbers in the range $0<k<k_{\text{max}}$, where $k_{\text{max}}=\sqrt{\omega^{2}-V^{\prime\prime}(\phi_{0})}$, and $\omega=\dot{\theta}$. The fastest growing modes of instability have a wavelength $\sim 10^{-2\pm 1}$ of the horizon size at the time of fragmentation, and they create isolated lumps of condensate which evolve into Q-balls. Numerical simulations Kasuya and Kawasaki (2000a, b) indicate that most of the condensate ends up in lumps. However, since the mass of Q-balls is a non-linear function of the Q-ball size, Q-ball formation, in general, leads to a non-uniform distribution of energy density in the matter component represented by the scalar condensate. Q-balls can also form when the charge density is small or zero, in which case both positively and negatively charged Q-balls are produced Kasuya and Kawasaki (2000a); here we do not consider this possibility.

Depending on the potential, the Q-balls with a global charge $Q$ have the following properties Dvali et al. (1998); Kusenko and Shaposhnikov (1998); Kusenko et al. (2005); Cotner and Kusenko (2016b):

[TABLE]

where $\Lambda$ is the energy scale associated with the scalar potential, $Q$ is the global U(1) charge and $0<\alpha<1$, $0<\beta<1$ denotes which type of $Q$-ball is under consideration (and also depends on the form of the scalar potential). For “flat direction" (FD) $Q$-balls, $\alpha=3/4$ ($\beta=1/4$), and for “curved direction" (CD) $Q$-balls, $\alpha=1$ ($\beta=1/3$) Dvali et al. (1998); Kusenko and Shaposhnikov (1998).

III Q-ball charge/mass distributions

Numerical simulations of condensate fragmentation and $Q$-ball formation have been performed in the past, from which we are able to determine the resulting charge and mass distributions Multamaki and Vilja (2002); Kasuya and Kawasaki (2000b). These distributions appear to be very sensitive to initial conditions in the condensate, such as the ratio of energy to charge density ($x=\rho/mq$), and to the details of the scalar potential. In addition, the resultant charge distribution can be very non-Gaussian due to the high degree of nonlinearity and chaos in the fragmentation process.

It should be understood that the results of these simulations are statistical in nature: a large number of $Q$-balls are created within the simulation volume so that the charge distribution tends towards a statistical average. In reality, if one were to perform a large simulation and look at the charge distributions in a number of small sub-volumes, you will find a large degree of variation, with more variation on smaller scales due to small sample sizes (this is not to say that the variance will be larger, just that the differences between distributions are large), as can be seen in Fig. 1.

It is these large fluctuations relative to the mean that will be the source of density perturbations.

Once the resultant charge distribution of Q-balls $f_{Q}(Q)\,dQ$ has been calculated from these numerical simulations, we can use this to calculate the mass distribution for single Q-balls using $M=\Lambda|Q|^{\alpha}$ (we will absorb all numerical factors into the definition of $\Lambda$ without loss of generality):

[TABLE]

It is important to note that a distribution well-localized in charge is also well-localized in mass. We can also use probability theory to calculate the mass of a collection of Q-balls. Under the assumption that a charge $Q_{\text{tot}}$ is distributed amongst $N$ Q-balls whose distribution is described by $f_{Q}(Q)$, the probability distribution function (PDF) for the total mass of this collection of Q-balls is given by

[TABLE]

where $\psi$ admits a simple-looking Fourier transform in $M$ and $Q_{\text{tot}}$:

[TABLE]

The power of $\alpha$ prevents analytic calculation of this PDF for all but the simplest charge distributions. Specifically, if we take the charge distribution to be a delta function: $f_{Q}(Q)=\delta(Q-Q_{0})$, then the mass distribution is also a delta function: $f_{M}(M,Q_{\text{tot}},N)=\delta(M-NM_{0})$, where $M_{0}=\Lambda Q_{0}^{\alpha}$, and $Q_{0}=Q_{\text{tot}}/N$ to satisfy charge conservation (this constraint comes from a mathematical issue that arises due to the canceling of a delta function of the form $\delta(Q_{\text{tot}}-NQ_{0})/\delta(Q_{\text{tot}}-NQ_{0})$; we can see that if we consider $\delta(x)$ as the limit of a smooth function that approaches this distribution, then this ratio is unity provided $Q_{0}=Q_{\text{tot}}/N$).

For ease of computation, we will assume the delta function charge/mass distribution for the rest of this paper. This also has good theoretical motivation, as the Affleck-Dine baryogenesis scenario requires a large nonzero charge density, which tends to result in a highly-localized charge distribution.

III.1 Single length scale

One should notice that the mass distribution function calculated earlier is also a function of both the total charge $Q_{\text{tot}}$ and the number of Q-balls $N$. During the chaotic fragmentation procedure, the number of Q-balls will fluctuate between horizons. So in order to get a full description of the fluctuations, we must supplement the mass distribution with a number distribution $p(N)$. This can be calculated from a simulation by simply counting the number of Q-balls within the simulation volume. Here, we will assume that the number of Q-balls per horizon $N$ is described by a Poisson distribution, as is typical for a random process such as fragmentation:

[TABLE]

where $N_{f}$ is the average number of Q-balls per horizon at fragmentation. We then combine Equations 4 and 6 to create a joint PDF which describes the distribution of mass $M$ within a horizon composed of $N$ Q-balls (we also set $Q_{\text{tot}}=Q_{f}$, the total charge on the horizon at $t_{f}$): $F_{Q}(M,N)=f_{M}(M,Q_{f},N)p(N)$. This is manifestly normalized since $\sum_{N}\int dM\,F_{Q}=\sum_{N}(\int dM\,f_{M})p=\sum_{N}p=1$. We can then use this to calculate statistical moments such as the average horizon mass, average horizon Q-ball number, RMS fluctuations, etc:

[TABLE]

where $M_{f}=\Lambda Q_{f}^{\alpha}N_{f}^{1-\alpha}$ is the horizon mass of Q-balls at $t_{f}$ (the first relation is only approximate because $\braket{N^{1-\alpha}}\approx N_{f}^{1-\alpha}$, though the relative error scales as $|\braket{N^{1-\alpha}}-N_{f}^{1-\alpha}|/\braket{N^{1-\alpha}}\approx 1/10N_{f}$, so totally negligible for large $N_{f}$).

III.2 Multiple length scales

The previous treatment has the shortcoming that it can only describe Q-ball distributions with spatial extent the size of the horizon at the time of fragmentation. We now generalize this to handle distributions on an arbitrary scale. First, when considering a physical volume $V$ at the time of fragmentation, the charge contained within this volume (assuming initial uniformity of the condensate) is given by $Q_{\text{tot}}=Q_{V}=Q_{f}(V/V_{f})$, where $V_{f}=\frac{4\pi}{3}t_{f}^{3}$ is the horizon volume at $t_{f}$. Second, the number distribution is altered so that the number of Q-balls within volume $V$ (assuming the same average number density $n_{f}=N_{f}/V_{f}$ across all scales) is described by

[TABLE]

The joint PDF for a mass $M$ composed of $N$ Q-balls contained within a volume $V$ at the time of fragmentation $t_{f}$ is then given by

[TABLE]

Note that if $V=V_{f}$, this reduces to the single-scale, horizon-size treatment.

In addition to being able to calculate quantities on each scale $V$ individually, we will also want to sum the contributions from each scale in some cases (such as contributions to the PBH density from both subhorizon and superhorizon modes). To do so, we will sum over all volume scales from $V_{\text{min}}$ to $V_{\text{max}}$ using a coarse-graining method. We consider an arbitrary function of volume $g(V)$. The sum of the contributions from each scale $V_{i}=V_{\text{max}}/\chi^{i-1}$ is then given by

[TABLE]

where we have used Euler-Maclaurin to approximate the sum, and $\chi\sim\text{few}$ is a parameter of the spacing between intervals of the coarse-graining procedure. We will take $\chi=e$ from now on for simplicity; another choice does not significantly affect the outcome provided it is not too close to unity. $V_{\text{min}}$ is the smallest volume under consideration; there will be a natural cutoff due to the fact that Q-balls have a finite size, and so this scale is generally defined as the volume which contains some number $N_{\text{min}}\sim 10$ Q-balls on average: $V_{\text{min}}/N_{\text{min}}=N_{f}/V_{f}$.

IV Q-ball and PBH densities

Using the framework of Section III, we are now in a position to begin calculating the energy densities associated with the Q-balls, fluctuations in that energy density, and the resulting density of black holes.

IV.1 Q-ball density at fragmentation

Using the formalism of Section III.2, we can calculate the background energy density (over the largest scales) of Q-balls at $t_{f}$:

[TABLE]

which will be important in the discussion of density perturbations in Section IV. Since Q-balls are formed at rest, the evolution of the Q-ball density after fragmentation is simply that of decaying nonrelativistic matter $\braket{\rho_{Q}(t)}=\braket{\rho_{Q}(t_{f})}(a_{f}/a)^{3}e^{(t_{f}-t)/\tau_{Q}}$, where $\tau_{Q}=1/\Gamma_{Q}$ is the lifetime of the Q-balls. Q-balls are generally considered stable with respect to decay into the quanta of the scalar field, but it is possible to decay through other processes. For example, if a coupling of the scalar field to a light fermion with mass $m<\omega$ exists, Q-balls can decay to these fermions through an evaporation process Cohen et al. (1986); Kawasaki and Yamada (2013); Enqvist and McDonald (1998b). Q-balls can also decay if the U(1) symmetry is broken by some higher-dimension operators Kusenko et al. (2005); Kawasaki et al. (2005); Kasuya and Kawasaki (2014); Cotner and Kusenko (2016b). We define $\Gamma_{Q}$ to include all such decay channels.

IV.2 Q-ball density perturbations due to fluctuations

Due to the stochastic nature of the fragmentation process, volumes of space can arise within which the number density of Q-balls exceeds the average number density. Due to the nonlinear relationship between Q-ball mass and charge $M=\Lambda|Q|^{\alpha}$, this also gives rise to fluctuations in the energy density within that volume. The density contrast in Q-balls at fragmentation $\delta(t_{f})$ for a volume $V$ containing mass $M$ is defined as

[TABLE]

where in the last line we have used the argument of the delta function in Equation 11 to eliminate $V$ (this will be justified by an integral over $V$ later). Note that if the Q-ball mass-charge relationship were linear ($\alpha=1$), the perturbations would vanish identically.

The subhorizon density perturbations ($V<V_{f}$) are frozen during the initial radiation dominated era, but they grow linearly in the scale factor during the Q-ball dominated epoch: $\delta(t)=\delta(t_{f})(a/a_{Q})=\delta(t_{f})=(t/t_{Q})^{2/3}$, where $t_{Q}$ is the beginning of the era of Q-ball domination. The structure growth generally goes nonlinear and decouples from the expansion around $\delta>\delta_{c}\sim 1.7$, at which point the overdense regions collapse and become gravitationally bound. However, some structures with $\delta<\delta_{c}$ can still collapse, and not all structures with $\delta>\delta_{c}$ are guaranteed to collapse into black holes. Due to nonsphericity of the gravitationally-bound structures, only a fraction $\beta=\gamma\delta^{13/2}(t_{R})(M/M_{Q})^{13/3}$ (where $\gamma\approx 0.02$ is a factor due to the nonsphericity, $M_{Q}=M_{f}(t_{Q}/t_{f})^{3/2}$ is the horizon mass at the beginning of the Q-ball dominated era, and $t_{R}$ is the end of the Q-ball dominated era, when the radiation comes to dominate again) will actually collapse to black holes Carr et al. (2016a); Polnarev and Khlopov (1985); Georg et al. (2016) by the end of the Q-ball dominated era. We assume that structures with $\delta\geq\delta_{c}$ do not continue to grow past the point of nonlinearity, as they have already collapsed and had their chance to form a PBH; for these perturbations we set $\beta=\gamma\delta_{c}^{13/2}(M/M_{Q})^{13/3}$ for $\delta(t_{R})>\delta_{c}$. This refinement may not be necessary, as the average density perturbations are generally so small they never reach $\delta_{c}$, and indeed, changing the value of $\delta_{c}$ does not seem to significantly alter the outcome.

Additional care must be taken to extend this to scales which enter the horizon at later times, and thus may not subject to the same amount of growth as subhorizon modes. Those that enter the horizon between $t_{f}<t<t_{Q}$ can be treated as effectively subhorizon since they enter the horizon before the Q-ball dominated epoch begins, and thus fluctuations are subject to the same amount of amplification as initially subhorizon modes. This includes all volumes $V<V_{Q}$, where $V_{Q}=\frac{4\pi}{3}t_{Q}^{3}(t_{f}/t_{Q})^{3/2}$ is defined as the (initially superhorizon) physical volume at $t_{f}$ which enters the horizon at $t_{Q}$: $(a_{Q}/a_{f})^{3}V_{Q}=\frac{4\pi}{3}t_{Q}^{3}$. Fluctuations which enter the horizon during the Q-ball dominated epoch are treated slightly differently, as they are only subjected to amplification from the time they enter the horizon $t_{h}$ until the radiation comes to dominate again at $t_{R}$. We can account for this by calculating $t_{h}$ for a given scale $V$ via $(a(t_{h})/a_{f})^{3}V=\frac{4\pi}{3}t_{h}^{3}$ (which gives us $t_{h}=(3/4\pi)(V/t_{f}^{3/2}t_{Q}^{1/2})$), and then replacing the scale factor $a_{R}/a_{Q}$ with $a_{R}/a(t_{h})$ in the definition of $\beta$ above. This treatment is valid for all scales between $V_{Q}<V<V_{R}$, where $V_{R}=(4\pi/3)t_{R}^{3}(t_{Q}/t_{R})^{2}(t_{f}/t_{Q})^{3/2}$ is the physical volume at $t_{f}$ which enters the horizon at $t_{R}$.

In addition to these details, we also enforce the constraint $\beta\leq 1$ in order to prevent PBH production probabilities over unity, though this does not become relevant unless the Q-ball dominated era is extremely long.

IV.3 Primordial black hole density

We are now in a position to calculate the average energy density in PBH created during the Q-ball dominated era. We do this by calculating the energy density of Q-balls at $t_{f}$ that will eventually form black holes by $t_{R}$ by weighting the Q-ball energy density $M/V$ by the collapse fraction/probability $\beta$ evaluated at $t_{R}$, summing over all scales $V$, and then redshifting this value appropriately. The expression for this procedure is given by

[TABLE]

where it should be understood that the integral over $V$ is broken up into two separate domains, $[V_{\text{min}},V_{Q}]$ and $[V_{Q},V_{R}]$, where separate definitions of $\beta$ apply, as described in Section IV.2. Due to the complicated piecewise nature of the function $\beta$, the authors are unaware of any analytic solution, and further progress must be made numerically.

We find that Equation 17 can be rewritten in such a way that it only depends on the dimensionless numbers $N_{f}$, $r_{f}=t_{Q}/t_{f}$, and $r=t_{R}/t_{Q}$. $r_{f}$ and $r$ can be interpreted as measures of the duration of the era between the fragmentation and the beginning of Q-ball domination, and the length of the Q-ball dominated era, respectively. The effect of these parameters on the black hole density can be seen in Figure 2.

Larger $r_{f}$ will reduce the fraction of Q-ball energy that goes into making black holes due to the dilution of the number density and increased horizon mass at $t_{Q}$ due to the delay of the Q-ball dominated era. Larger $r$ leads to an increased fraction of Q-ball energy that goes into black holes due to more amplification of the density perturbations, leading to a higher probability of PBH formation. Larger $N_{f}$ reduces the fraction because of higher suppression of fluctuations due to the Poisson statistics. The form of the contours in this plot suggest that this ratio roughly scales as $\braket{\rho_{\text{BH}}}/\braket{\rho_{Q}}\sim(4.6\times 10^{-4})r_{f}^{-5.9}r^{4.4}N_{f}^{-3.7}$.

IV.4 Black hole mass spectrum

One can derive the mass spectrum of the black holes by not integrating over $M$ in Equation 17:

[TABLE]

This yields the differential black hole energy density $d\braket{\rho_{\text{BH}}}/dM$. We find that the spectrum can be rewritten in terms of the parameter $\eta=M/M_{f}$ (fraction of horizon mass at $t_{f}$), along with the previously mentioned parameters $r_{f}=t_{Q}/t_{f}$, $r=t_{R}/t_{Q}$, and $N_{f}$. Calculation of this function can be done by evaluating Equation 18 at multiple values of $\eta$ and then interpolating. An example is given in Figure 3.

First, it’s obvious from the normalization of each curve that the lower the number of Q-balls per horizon, the more black holes that are created. This is expected, as the Poisson statistics suppress the density fluctuations for large Q-ball number. The normalization also increases with $r$, as explained in Figure 2. Second, there is a hard lower cutoff in the PBH mass, which occurs at $\eta=N_{\text{min}}/N_{f}$, which is due to the lower cutoff in the volume mentioned earlier. Above that, the BH number sharply increases with a power law $\propto\eta^{2.85\pm 0.15}$; the extent of this region depends on the magnitude of $r$, with larger values leading to a larger range. We suspect that this is due to the fact that the small-scale density fluctuations have already reached their critical value $\delta_{c}$ and can no longer continue growing, whereas the large-scale fluctuations (which started out smaller) still have room to do so. Above that, the spectrum becomes approximately flat ($\propto\eta^{-0.15}$), meaning that the number of black holes in each decade of mass are comparable. Of course, the upper end of this range dominates the energy density of the distribution. Then, at around $M=M_{Q}$, there is a sharp transition and the slope becomes strongly negative ($\propto\eta^{-4.5}$) due to the reduced growth the superhorizon modes are subject to. Then, there is an upper exponential cutoff at $\eta\sim 10^{8}/N_{f}$ due once again to the Poisson statistics (the cutoff appears to take precedence over previously mentioned transitions).

V Cosmological history

We now give a detailed account of how the Q-balls, radiation, and black holes evolve throughout the history of the universe up until the present day, as seen in Figure 4.

In summary: we assume an initial period of inflation and reheating in order to create a radiation dominated era with a uniform charged scalar field as a subdominant component of the energy density. The scalar field fragments into Q-balls at $t_{f}$, which then come to dominate the energy density at $t_{Q}$. During the Q-ball dominated epoch, primordial black holes are produced, and at $t_{R}$, their density is frozen in and evolves as nonrelativistic matter. After this initial matter dominated epoch, the Standard Model of cosmology resumes, and evolves through all the eras we are familiar with (BBN, matter-radiation equality, etc.) up to the present day.

The functions used to model the energy density evolution for each species is summarized in Appendix A.

V.1 Initial radiation dominated era

After the end of inflation, the Universe enters a brief matter dominated era due to the coherent oscillations of the inflaton field. The decay of the quanta of this field at time $t_{\text{RH}}=\Gamma_{I}^{-1}$ reheats the Universe, which enters a radiation dominated epoch with temperature $T_{\text{RH}}=0.55g_{*}^{-1/4}(\Gamma_{I}M_{p})^{1/2}$ and radiation energy density

[TABLE]

where $\Gamma_{I}\sim 1/t_{\text{RH}}$ is the decay width of the inflaton oscillations. The radiation density redshifts as $\rho_{R}(t)=\rho_{R}(t_{\text{RH}})(a_{\text{RH}}/a)^{4}=\rho_{R}(t_{\text{RH}})(t_{\text{RH}}/t)^{2}$ during this epoch, which ends up canceling the factor of $\Gamma_{I}$ to give us

[TABLE]

At some point $t_{f}$, the scalar condensate fragments into Q-balls, resulting in an energy density given by Equation 15. The Q-balls then redshift as decaying nonrelativistic matter:

[TABLE]

V.2 Q-ball dominated era

At some point $t_{Q}$, the Q-balls come to dominate the energy density. This time is defined by $\rho_{R}(t_{Q})=\braket{\rho_{Q}(t_{Q})}$, using the equations of the previous section. During this era, Q-ball decays begin to affect the radiation density, causing the radiation temperature to decrease less slowly than it normally would due to the expansion. Following the analysis of Scherrer and Turner Scherrer and Turner (1985), the radiation density in this epoch due to the decay of the Q-balls can be modeled as

[TABLE]

where $x\equiv\Gamma_{Q}t$, $x_{0}=\Gamma_{Q}t_{Q}$, and $z=(x/x_{0})^{2/3}$. The Q-balls continue to redshift and decay, leading to

[TABLE]

As the Q-balls decay, eventually the radiation comes to dominate again at $t_{R}$, defined by $\rho_{R}(t_{R})=\braket{\rho_{Q}(t_{R})}$. Using Equations 20, 21, 23 and 24, this gives us the relation

[TABLE]

where $\Gamma$ is the generalized incomplete gamma function. This allows us to solve (numerically) for $r_{Q}\equiv\tau_{Q}/t_{Q}$ as a function of $r=t_{R}/t_{Q}$. At this point, if we specify $t_{f}$, $r_{f}$ and $r$, we can calculate the other parameters $t_{Q}$, $t_{R}$, $\tau_{Q}$, and $\Lambda Q_{f}^{\alpha}$ via

[TABLE]

from which we can calculate all other quantities of interest ($M_{f}=\Lambda Q_{f}^{\alpha}N_{f}^{1-\alpha}$, $M_{Q}=M_{f}r_{f}^{3/2}$, etc.).

V.3 Standard cosmological era

After the Q-balls have decayed sufficiently, the universe returns to a radiation dominated era, and the standard cosmology begins. In order to evolve the radiation, Q-ball, and black hole densities to the present day, one would naïvely use $a_{1}/a_{2}=(t_{1}/t_{2})^{n}$, where $n=\frac{1}{2}\,(\frac{2}{3})$ in a radiation (matter) dominated era, keeping in mind that the universe transitions between the two at $z_{\text{eq}}\approx 3360$, or $t_{\text{eq}}\approx 4.7\times 10^{4}\text{ yr}$. However, due to the extended era of matter domination, the time at which cosmological events (such as BBN, matter-radiation equality, or recombination) occur are not the same as in the standard cosmology. Instead, one must evolve according to the universe’s thermal history, where cosmological events occur at specific temperatures. In this case, one must use $a_{1}/a_{2}=g_{*S}^{1/3}(T_{2})T_{2}/g_{*S}^{1/3}(T_{1})T_{1}$ and evolve from $T_{R}$ (defined by $\rho_{R}(t_{R})=(\pi^{2}/30)g_{*}(T_{R})T_{R}^{4}$) to $T_{0}=2.7\text{ K}=2.3\text{ meV}$. This has the advantage of accurately accounting for any deviation from cosmological history. We can then find the time at which some event $X$ occurs by solving $\rho_{R}(t_{X})=(\pi^{2}/30)g_{*}(T_{X})T_{X}^{4}=\rho_{R}(T_{R})(a(t_{R})/a(t))^{4}$. In order to ensure that this early matter dominated era does not spoil the canonical cosmological thermal history, we enforce an additional constraint $T_{R}>T_{\text{BBN}}\sim\text{MeV}$, so that the entropy injection from Q-ball decays does not interfere with nucleosynthesis.

VI Observational constraints

We now examine the observational constraints on primordial black holes and where our results fit in. The constraints come from a wide variety of sources, such as extragalactic gamma rays from evaporation Johnson et al. (2017); Carr et al. (2016b), femtolensing of gamma ray bursts (GRB) Barnacka et al. (2012), capture by white dwarfs Graham et al. (2015), microlensing observations from HSC Niikura et al. (2017), Kepler Griest et al. (2013, 2014), and EROS/MACHO/OGLE Tisserand et al. (2007), measurements of distortion of the CMB Tashiro and Sugiyama (2008); Ricotti et al. (2008), and bounds on the number density of compact X-ray objects Inoue and Kusenko (2017) (constraints summarized in Carr et al. (2010); Frampton et al. (2010); Kuehnel and Freese (2017)). The constraints are typically expressed in a form that assumes a monochromatic distribution of PBH masses. However, in the case of an extended mass distribution (such as we have in this scenario), care must be taken to apply the limits correctly. To do this, we follow the procedure outlined in Carr et al. (2016a); Kuehnel and Freese (2017), which amounts to dividing the mass spectrum into a number of bins (labeled by the index $i$), then integrating the dark matter fraction over the interval contained in the bin:

[TABLE]

which is then compared with the constraints on a bin-by-bin basis. We find that for sufficient choices of the parameters $N_{f}$, $t_{f}$, $r_{f}$, and $r$, our model can produce black holes over practically the entire parameter space allowed by the constraints (see Fig. 5).

Notably, this Figure illustrates two interesting points: 1) that this mechanism is capable of generating black holes which can account for both 100% of the dark matter in the region $M\sim 10^{20}\text{ g}$ and production or r-process elements Fuller et al. (2017), and 2) it is also capable of generating black holes with sufficient mass to explain the recent LIGO observation GW150914 Abbott (2016). Some studies have even argued that PBH can account for 100% of the DM in this range by contesting the CMB constraints Bird et al. (2016). The three contours in Figure 5 are, however, simply chosen by hand for illustrative purposes, and are therefore not representative of the entire parameter space allowed to this mechanism, which is much wider than suggested by the given parameters.

VII Parameter space

We now explore the parameter space available to this mechanism in which it is possible to account for a considerable fraction of the dark matter while avoiding observational constraints. To do so, we develop an algorithm to accomplish this task in the following manner: since the function $((a_{f}/a_{R})^{3}/\braket{\rho_{Q}(t_{f})})\,(d\braket{\rho_{\text{BH}}(t_{R})}/d\eta)$ is determined solely by the parameters $r_{f}$ and $r$, we generate a list of such functions by sampling the $r-r_{f}$ plane at various points. Then, for each $(r,r_{f})$ pair, we vary $t_{f}$ using a weighed bisection method until $\max(|1-f_{i}/f_{\text{con}}(M_{i})|)<\epsilon=10^{-1}$, with $f_{i}$ given by Equation 28. This determines the value of $t_{f}$ which gives the maximum dark matter fraction allowed by the constraints for the given values of $r$ and $r_{f}$. Once that has been determined, we can calculate all other relevant quantities of interest, such as $M_{\text{BH,peak}}$, $T_{R}$, and the dark matter fraction $f=\Omega_{\text{PBH}}/\Omega_{\text{DM}}$. The results for $N_{f}=10^{6}$ are shown in Figure 6.

We can see from this Figure that the contours of constant $f$ are highly correlated with the values of $M_{\text{BH,peak}}$. This is due to the fact that the observational constraints are solely a function of $M$, and the calculated spectrum $f(M)$ is quite sharply peaked at $M_{\text{BH,peak}}$. The region of $f\geq 1$ roughly follows the contour of $M_{\text{BH,peak}}=10^{20}\text{ g}$, where the constraints are weakest, due to this fact. The general tendency appears to be that for increasing $r$, the spectrum favors heavier black holes, while for increasing $r_{f}$, it favors lighter black holes. This is because the longer the period of structure formation, the more the large mass perturbations (which initially have small $\delta$) grow, increasing the probability that they collapse to black holes, whereas the longer the period between fragmentation and Q-ball domination, the more the energy density of Q-balls is diluted, so that the value of $t_{f}$ has to be lower in order to achieve the same value of $f$ that one would with smaller $r_{f}$. The smaller $t_{f}$, the smaller the horizon mass, and the smaller the PBH masses. This can be quantified, as the contours of constant $f$ very roughly follow $r^{0.63}/r_{f}\sim\text{const}$.

However, at around $r\gtrsim 3.2\times 10^{5}$, the character of the plot changes so that the spectrum appears to no longer depend on $r$, only on $r_{f}$. This is likely due to the definition of $\beta$ discussed in Section IV.2, in which density perturbations are not allowed to grow past $\delta_{c}$, which makes it so that further increases in $r$ have no effect on the spectrum. Further increases in $r_{f}$, however, serve to further dilute the Q-ball density before structure formation can occur, lowering the necessary $t_{f}$ as mentioned above.

The pink region where $T_{R}>1\text{M}eV$ is ruled out since in this region Q-ball decays begin to interfere with nucleosynthesis. We can see that it is correlated with high mass, as the later the fragmentation (meaning larger horizon mass), the less time there is for the Q-balls to dominate before the universe cools sufficiently enough that nucleosynthesis begins. One might notice that the mass of the black holes seen at LIGO in the GW150914 event lies beyond the range indicated in this Figure, but in Figure 5 we have plotted a contour that has $M_{\text{BH,peak}}\sim 30\text{ M}_{\odot}$. This is because the contour in Fig. 5 is for $N_{f}=10^{5}$, whereas Fig. 6 has $N_{f}$ fixed at $10^{6}$. We can see in Fig. 3 that the density spectrum has a peak at higher masses for lower values of $N_{f}$, so that the equivalent plot to Fig. 6 for $N_{f}<10^{6}$ would see the $T_{R}<1\text{ MeV}$ constraint pushed to higher values of $M_{\text{BH,peak}}$; enough so that they could explain the $30\text{ M}_{\odot}$ black holes while avoiding the nucleosynthesis constraint.

VIII Topological defects

Topological defect formation can also lead to the production of PBHs if the topological defects come to dominate the energy density. The analysis is sufficiently different from that of Q-balls, primarily because typically only one defect per horizon is produced at the time of formation due to the Kibble mechanism Kibble (1976). However, the general mechanism remains the same: small number densities of defects lead to large fluctuations relative to the background density, these fluctuations become gravitationally bound and collapse to form black holes once the relic density has come to dominate, and the relics decay due to some instability (such as gravitational waves or decay to Nambu-Goldstone bosons in the case of cosmic strings). In order to accurately model production of PBHs from these defects, one should calculate the expected density perturbations on initially superhorizon scales, which only begin to grow once these scales pass back within the horizon and the defects come to dominate the universe’s energy density.

Cosmic strings are probably the most likely candidate for primordial relics due to the fact that they are typically cosmologically safe, as the energy density in string loops is diluted during expansion at the same rate as radiation, $a^{-4}$ Vachaspati and Vilenkin (1984); Turok (1983). In contrast, the string “network" (i.e. infinitely long strings) energy density redshifts as $a^{-2}$ so that they quickly come to dominate the universe’s energy density. However, once long strings start intercommuting to produce loops, and these loops subsequently self-intersect to fragment into smaller loops, the string network approaches a scaling solution which leads to the $a^{-4}$ dilution of string loops Bennett (1986).

However, this scaling solution critically relies on the probability of string intercommutation being very close to unity so that the long strings can efficiently break into small loops. If this probability was sufficiently low, then the string density could redshift as $a^{-2}$ or $a^{-3}$, survive until a matter/string dominated era, initiate structure formation, collapse to form PBHs, and then subsequently decay. As long as these conditions are satisfied, cosmic strings could act as a source of PBH. In addition, there exists a large class of solutions to the string equations of motion which never self-intersect Kibble and Turok (1982), making this scenario plausible.

IX Discussion

The mechanism we have discussed has a number of advantages over some other models. It is extremely effective in creating primordial black holes across a broad range of masses, and it does not require the tuning the inflaton potential Inomata et al. (2017); Garcia-Bellido et al. (1996); Inomata et al. (2016); Kawasaki et al. (2016). We did not have to make any ad-hoc assumptions regarding density fluctuations; the fluctuations are calculable from first principles.

The mechanism is also generally applicable to practically any complex scalar field with a conserved global charge and flat potential, so that the formation of PBH is now a general prediction of any theory containing such charged scalars. In particular, supersymmetric extensions to the Standard Model typically have such fields, making the production of primordial black holes a general prediction of such theories. For the case of supersymmetric Q-balls with the SUSY-breaking scale $\Lambda_{\text{SUSY}}>10\text{ TeV}$, the fragmentation time cannot be much longer than the Hubble time $H^{-1}\sim M_{p}/g_{*}^{1/2}\Lambda_{\text{SUSY}}\lesssim 8\times 10^{-15}\text{ s}$, which corresponds to peak PBH masses of about $10^{23}\text{g}$ (assuming $N_{f}\sim 10^{6}$). The solid black curve illustrated in Figure 5 satisfies this bound, thus primordial black holes from supersymmetric Q-balls can account for 100% of the dark matter.

Supersymmetric Q-balls themselves have been suggested as the source of dark matter in models where they are entirely stable Kusenko et al. (2005); Cotner and Kusenko (2016b); Kusenko (1997a); Enqvist and McDonald (1998b); Kusenko and Shaposhnikov (1998). However, the stability is model-dependent, and it only applies to Q-balls carrying baryonic charge (so-called B-balls), since those carrying leptonic charge (L-balls) would quickly evaporate to neutrinos Cohen et al. (1986). In our scenario, a short evaporation timescale is precisely what is needed to end the early era of Q-ball domination before nucleosynthesis begins. The L-balls would then be composed of slepton fields which subsequently decay to neutrinos at an early time. Since neutrinos do not decouple from the plasma until just before nucleosynthesis ($T\gtrsim\text{few MeV}$), they thermalize quickly. If they decay early enough, it may even be possible to generate the baryon asymmetry through conversion of lepton number to baryon number via sphaleron processes during the electroweak phase transition Fujii et al. (2001).

There are some remaining open questions, such as how well the assumption that all Q-balls within the volume $V$ are the same charge models this scenario. Simulations show that for scalar condensates with a high ratio of charge density to energy density, this is a good prediction, as all the Q-balls formed from this initial condition typically have similar sizes. This is also theoretically understood for a scalar condensate with a sufficiently large charge density Kusenko and Shaposhnikov (1998); Enqvist and McDonald (1998a); Dine and Kusenko (2003). However, condensates with a large energy density and small charge density generally produce broad charge distributions, nearly symmetric about $Q=0$, since the excess energy cannot be contained in Q-balls with the same sign of $Q$ while also conserving charge. We suspect that in this scenario the production of PBH will be reduced, since the charge conservation does not play as significant a role. Loss of energy due to scalar radiation in the fragmentation process may still be able to produce energy density inhomogeneities, but this will require further study. In the same manner, the production of oscillons from the fragmentation of a real scalar field may be able to produce significant numbers of PBH as well.

One may also wonder what sort of mechanism is needed in order to ensure the Q-balls decay at the correct time. As an example, following the work of Cotner and Kusenko (2016b), the lifetime of a Q-ball with initial charge $Q_{0}$ decaying to pseudo-Goldstone bosons through the effects of a charge-violating operator of the form $V_{\not}{Q}(\phi)=g\phi^{n}(\phi^{*})^{m}/\Lambda_{*}^{n+m-4}+\text{c.c.}$ is given by

[TABLE]

where

[TABLE]

and $J_{nm}\sim O(10^{-7}-10^{-6})$. For $g\sim 0.1$, $\Lambda_{*}\sim 10^{16}\text{ GeV}$, and $\Lambda\sim 10^{9}\text{ GeV}$, the lifetime of a Q-ball decaying through an operator with $(n,m)=(2,3)$ is about $\tau\sim 10^{-13}\text{ s}$, which is sufficient to explain the curve of Figure 4 (and satisfies the SUSY bound). Decay through these higher-dimension operators isn’t the only way to induce the decay of Q-balls though; many other scenarios have been explored in the literature Kusenko et al. (2005); Kawasaki et al. (2005); Kasuya and Kawasaki (2014).

This work also begs the question of what possible observables exist that could show the Q-ball clusters collapse to black holes. We assume that the collapse will produce a stochastic gravitational wave background Nakama et al. (2017), which could be detected by future observatories (or put constraints on the model). Further evolution of the PBH population could see successive mergers, which in addition to creating another stochastic GW background García-Bellido (2017), could also alter the distribution of black hole masses (in addition to evaporation/accretion effects Bettwieser and Glatzel (1981); Custodio and Horvath (1998)). We propose to calculate the gravitational wave spectrum in a future publication.

X Conclusion

In summary, we have shown that the number density fluctuations of a Q-ball population in the early universe can lead to production of primordial black holes with sufficient abundance to explain the dark matter. Scalar fields and Q-ball formation are general features of supersymmetric extensions to the Standard Model, which provides a good motivation for this mechanism. A similar mechanism using solitons, topological defects, or other compact objects associated with scalar fields in the early universe can also lead to a copious production of primordial black holes.

XI Acknowledgements

This work was supported by the U.S. Department of Energy Grant No. DE-SC0009937. A.K. was also supported by the World Premier International Research Center Initiative (WPI), MEXT, Japan.

Appendix A Evolution of energy densities

Here we tabulate the functional form of the energy density for each species (radiation, Q-ball, black hole) up until the present day. The values for $t_{f}$, $t_{Q}$, $t_{R}$, and $\tau_{Q}$ are taken as input parameters (subject to some self-consistency conditions), while the values of $t_{\text{eq}}$ and $t_{0}$ are calculated from the procedure described in Section V.3.

A.1 Radiation

The radiation density begins after reheating and is given by Equation 20. From this point we evolve it through time to the present day, taking into account the contribution due to Q-ball decays during the period $t_{Q}<t<t_{R}$.

[TABLE]

A.2 Q-balls

The Q-balls are created at the time of fragmentation $t_{f}$, and evolve as decaying nonrelativistic matter. The magnitude of the energy density becomes insignificant shortly after $t_{R}$. $M_{f}=\Lambda|Q_{f}|^{\alpha}N_{f}^{1-\alpha}$ can be determined from specifying $t_{f}$, $r_{f}$, $r$, and $N_{f}$, as given in Section V.

[TABLE]

A.3 Black holes

The black holes are created towards the end of the initial Q-ball dominated era, and their density at $t_{R}$ is given by Equation 17:

[TABLE]

where

[TABLE]

and

[TABLE]

where $x_{\text{min}}=N_{f}V_{\text{min}}/V_{f}=N_{\text{min}}$, $x_{Q}=N_{f}V_{Q}/V_{f}$, and $x_{R}=N_{f}V_{R}/V_{f}$. After this has been evaluated, the evolution of the black hole density is fairly straightforward:

[TABLE]