Analytic Description of Primordial Black Hole Formation from Scalar Field Fragmentation

TL;DR

This paper presents an analytic framework for primordial black hole formation from scalar field fragmentation during the early universe, highlighting its differences from other models and its potential role as dark matter.

Contribution

It introduces a novel analytic approach to describe PBH formation via scalar field fragmentation, emphasizing the role of scalar self-interactions and matter fluctuations.

Findings

Scalar field fragmentation can produce PBHs during a matter-dominated era.

The process involves scalar lumps like Q-balls or oscillons forming from instabilities.

PBHs formed in this scenario are a plausible dark matter candidate.

Abstract

Primordial black hole (PBH) formation is a more generic phenomenon than was once thought. The dynamics of a scalar field in inflationary universe can produce PBHs under mild assumptions regarding the scalar potential. In the early universe, light scalar fields develop large expectation values during inflation and subsequently relax to the minimum of the effective potential at a later time. During the relaxation process, an initially homogeneous scalar condensate can fragment into lumps via an instability similar to the gravitational (Jeans) instability, where the scalar self-interactions, rather than gravity, play the leading role. The fragmentation of the scalar field into lumps (e.g. Q-balls or oscillons) creates matter composed of relatively few heavy "particles", whose distribution is subject to significant fluctuations unconstrained by comic microwave background (CMB) observations and unrelated to the large-scale structure. If this matter component comes to temporarily dominate the energy density before the scalar lumps decay, PBHs can be efficiently produced during the temporary matter-dominated era. We develop a general analytic framework for description of PBH formation in this class of models. We highlight the differences between the scalar fragmentation scenario and other commonly considered PBH formation models. Given the existence of the Higgs field and the preponderance of scalar fields within supersymmetric and other models of new physics, PBHs constitute an appealing and plausible candidate for dark matter.

Figures

Tables

Table 1. Table 1 : Comparative summary of standard PBH formation mechanism from inflationary perturbations and PBH production from scalar field fragmentation mechanism.

	PBH Production Scenario
	Inflationary Perturbations	Field Fragmentation
	(common mechanism)	(our mechanism)
Source and type of large	inflaton fluctuations,	inflaton fluctuations,
(CMB-scale) perturbations	curvature	curvature
Source and type of small	inflaton fluctuations,	stochastic field fragmentation,
(PBH-scale) perturbations	curvature	isocurvature (fragment-lumps)
PBH source field	inflaton	inflaton or spectator field
Required potential condition	inflaton potential fine tuning	no new restrictions on inflaton
		potential, scalar field potential
		shallower than quadratic
		(attractive self-interactions)
PBH formation era $\left(t_{\mathrm{PBH}}\right)$	$t_{\mathrm{BBN}}\gtrsim t_{\mathrm{PBH}}\gtrsim t_{\mathrm{reh}}$,	$t_{\mathrm{BBN}}\gtrsim t_{\mathrm{PBH}}\gtrsim t_{\inf }$,
and type	after reheating,	before or after reheating,
	radiation-dominated era	temporary matter-dominated era
PBH size $\left(r_{\mathrm{BH}}\right)$ vs. horizon $\left(r_H\right)$	$r_{\mathrm{BH}}\sim r_H\sim H^{-1}$	$r_{\mathrm{BH}}\ll r_H\sim H^{-1}$
at formation
PBH spin $\left(a_s\right)$	$a_s\sim 0$	$a_s\sim 𝒪\left(1\right)$ possible

Table 2. Table 2 : Model input parameters for production of PBHs from clustering of Q-balls. As explained in the text, Model A is consistent with electroweak-scale supersymmetry.

Model	$t_f$ [s]	$r_f$	$r$	${\overline{N}}_H$	$\alpha$	$ϵ$
A	$3\times {10}^{-16}$	2	$100$	$100$	$3/4$	1
B	$7\times {10}^{-9}$	2	$40$	$30$	$3/4$	1
C	$4\times {10}^{-3}$	1.1	$30$	$30$	$3/4$	1

Table 3. Table 3 : Derived quantities for models A, B, C described in Table 2 .

		$T_R$	$M_{\text{Q-ball}}$	$M_{\text{PBH,peak}}$	$M_{\text{PBH,peak}}$
Model	$f_{\text{DM}}$	[MeV]	[GeV]	[g]	[M⊙]
A	1.0	$9.0\times {10}^5$	$6.1\times {10}^{42}$	$7.8\times {10}^{20}$	$3.9\times {10}^{-13}$
B	$1.2\times {10}^{-2}$	$3.4\times {10}^2$	$4.8\times {10}^{50}$	$1.7\times {10}^{28}$	$8.4\times {10}^{-6}$
C	$8.5\times {10}^{-4}$	$1.1$	$3.6\times {10}^{56}$	$6.3\times {10}^{33}$	$3.1$

Equations

$$U_{\mathrm{total}}=U_I({\Phi }_I)+U(\varphi )\approx U_I({\Phi }_I).$$

$$U(⟨|\varphi |⟩)\sim H_I^4,$$

$$H_I=\sqrt{\frac{U_{\mathrm{total}}}{3M_{pl}^2}}\approx \sqrt{\frac{U_I({\Phi }_I)}{3M_{pl}^2}}$$

$$0<k<k_{\max }=k_{\mathrm{J}}=2(\pi G\rho {)}^{1/2}/c_s,$$

$$\mathcal{L}=|{\partial }_{\mu }\varphi {|}^2-U(\varphi ).$$

$$\rho =|{\partial }_{\mu }\varphi {|}^2+U(\varphi )$$

$$p=|{\partial }_{\mu }\varphi {|}^2-U(\varphi )$$

$$p=\Big{(}\frac{n-2}{n+2}\Big{)}\rho$$

$$\varphi (x,t)=\varphi (t)\equiv \Phi (t)e^{i\Theta (t)}.$$

$$\displaystyle\left\{\begin{array}[]{lll}\ddot{\Theta}+3H\dot{\Theta}-\dfrac{1}{a^{2}(t)}\Delta\Theta+\Big{(}\dfrac{2\dot{\Phi}}{\Phi}\Big{)}\dot{\Theta}-\dfrac{2}{a^{2}(t)\Phi}(\partial_{i}\Theta)(\partial^{i}\Phi)&=&0~{},\\ \ddot{\Phi}+3H\dot{\Phi}-\dfrac{1}{a^{2}(t)}\Delta\Phi-\dot{\Theta}^{2}\Phi+\dfrac{1}{a^{2}(t)}(\partial_{i}\Theta)^{2}\Phi+\Big{(}\dfrac{\partial U}{\partial\Phi}\Big{)}&=&0~{},\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ll}&\ddot{\delta\Theta}+3H\dot{(\delta\Theta)}-\dfrac{1}{a^{2}(t)}\Delta(\delta\Theta)+\Big{(}\dfrac{2\dot{\Phi}}{\Phi}\Big{)}\dot{(\delta\Theta)}+\Big{(}\dfrac{2\dot{\Theta}}{\Phi}\Big{)}\dot{(\delta\Phi)}-\Big{(}\dfrac{2\dot{\Phi}\dot{\Theta}}{\Phi^{2}}\Big{)}\delta\Phi~{}=~{}0~{},\\ &\ddot{\delta\Phi}+3H\dot{(\delta\Phi)}-\dfrac{1}{a^{2}(t)}\Delta(\delta\Phi)-2\Phi\dot{\Theta}\dot{(\delta\Theta)}+U^{\prime\prime}\delta\Phi-\dot{\Theta}^{2}\delta\Phi~{}=~{}0~{}.\end{array}\right.$$

$$\displaystyle\left[\alpha^{2}+3H\alpha+\frac{k^{2}}{a^{2}}+\Big{(}\frac{2\dot{\Phi}}{\Phi}\Big{)}\alpha\right]\left[\alpha^{2}+3H\alpha+\frac{k^{2}}{a^{2}}+U^{\prime\prime}(\Phi)-\omega^{2}\right]+4\omega^{2}\left[\alpha-\frac{\dot{\Phi}}{\Phi}\right]\alpha=0~{}.$$

$$\frac{k^2}{a^2}\left[\frac{k^2}{a^2}+U^{\prime \prime }(\Phi )-{\omega }^2\right]=0,\mathrm{with}\alpha =0.$$

$$U^{\prime \prime }(\Phi )-{\omega }^2<0,$$

$$\underline{\mathrm{unstable}\mathrm{modes}}:\alpha (k)>0,\mathrm{for}0<k<k_{\max }=\sqrt{{\omega }^2-U^{\prime \prime }(\Phi )}.$$

$$U^{\prime \prime }(\Phi )<0.$$

$$R_Q\sim k_{\mathrm{nl}}^{-1}\sim f_Q{H}^{-1},$$

$$\overline{n}\sim\bigg{(}\dfrac{k_{nl}}{2\pi}\bigg{)}^{3}~{},$$

$$\delta =\frac{\rho -⟨\rho ⟩}{⟨\rho ⟩}>0,$$

$${⟨{\rho }_{\text{PBH}}⟩}_V={⟨B{\rho }_{\mathrm{M}}⟩}_V={⟨BM/V⟩}_V=\frac{1}{V}\int dMP(M|V)B(M,V)M,$$

$$⟨{\rho }_{\text{PBH}}⟩=\int d(\ln V){⟨B{\rho }_{\mathrm{M}}⟩}_V=\int \frac{dV}{V^2}\int dMP(M|V)B(M,V)M.$$

$$P(M|V)=\sum _{N}P(M|N)P(N|V),$$

$$P(M|V)=$$

$$=$$

$$P(M|V)=\left(\prod _i\sum _{Q_i}\right)\left(\prod _j\int d{\xi }_j\right)P(M|\{Q\},\{\xi \})P(\{Q\},\{\xi \}|V).$$

$$P(M|V)=\left(\prod _i\sum _{Q_i}P(Q_i|V)\right)\left(\prod _j\int d{\xi }_{j}P({\xi }_j|V)\right)P(M|\{Q\},\{\xi \}).$$

$$P_{\mathrm{id}}(M|\{Q\},\{\xi \})=\delta (M-M_f(m|\{Q\},\{\xi \}))P(m|\{Q\},\{\xi \}),$$

$$P_{\mathrm{id}}(M|N)=\int dm\delta (M-mN)P(m)=\frac{P(M/N)}{N},$$

$$P_{\mathrm{iid}}(M|N)$$

$$=\int \frac{d\mu }{2\pi }{e}^{-i\mu M}\int dM{e}^{i\mu M}\left(\prod _{i=1}^N\int d{m}_{i}P(m_i)\right)\delta \left(M-\sum _i{m}_i\right)$$

$$=\int \frac{d\mu }{2\pi }{e}^{-i\mu M}\left(\prod _{i=1}^N\int d{m}_{i}P(m_i)e^{i\mu m_i}\right)$$

$$=\int \frac{d\mu }{2\pi }{e}^{-i\mu M}{\tilde{P}}^N(\mu ).$$

$$P_{\mathrm{s}}(m)=\frac{e^{-(m-m_0{)}^2/2{\sigma }^2}}{\sqrt{2\pi {\sigma }^2}},$$

Full text

††preprint: IPMU19-0063, YITP-19-31

1 Introduction

Standard astrophysical black holes form as a result of stellar collapse at relatively recent cosmological times, after the onset of star formation. However, black holes could also form in the early universe. Such primordial black holes (PBHs) can account for all or part of the dark matter (DM) (e.g. Zeldovich:1967 ; Hawking:1971ei ; Carr:1974nx ; GarciaBellido:1996qt ; Khlopov:2008qy ; Frampton:2010sw ; Kawasaki:2016pql ; Cotner:2016cvr ; Carr:2016drx ; Inomata:2016rbd ; Pi:2017gih ; Inomata:2017okj ; Garcia-Bellido:2017aan ; Inoue:2017csr ; Georg:2017mqk ; Inomata:2017bwi ; Kocsis:2017yty ; Ando:2017veq ; Sasaki:2018dmp ; Carr:2018rid ). PBHs can be associated with a variety of astrophysical phenomena, including the recently discovered Abbott:2016blz ; Abbott:2016nmj ; Abbott:2017vtc gravitational waves (e.g. Nakamura:1997sm ; Clesse:2015wea ; Bird:2016dcv ; Raidal:2017mfl ; Eroshenko:2016hmn ; Sasaki:2016jop ; Clesse:2016ajp ; Takhistov:2017bpt ), formation of supermassive black holes Bean:2002kx ; Kawasaki:2012kn ; Clesse:2015wea as well as $r$-process nucleosynthesis Fuller:2017uyd , gamma-ray bursts and microquasars Takhistov:2017nmt from compact star disruptions.

In conventional PBH production mechanisms (see e.g. Ref. Carr:2016drx ; Sasaki:2018dmp for review) BHs form when density fluctuations are of $\mathcal{O}(1)$ at the horizon crossing during re-entry after the end of inflation and reheating. The usual curvature perturbations from inflation are expected to be nearly scale invariant and are constrained by cosmic microwave background (CMB) observations Aghanim:2018eyx at large scales. To allow formation of PBHs in this class of models, one must dramatically enhance the perturbations on some smaller scale, as determined by the PBH mass. This can be accomplished by introducing new features within the inflaton potential, usually at the expense of significant fine-tuning GarciaBellido:1996qt ; Kawasaki:2018daf ; Germani:2018jgr . A hallmark of this class of models is that the PBH mass is close to the horizon mass at the time when large perturbations reenter the horizon. Some production mechanisms suggested invoking metastability of the Standard Model vacuum Espinosa:2017sgp . We will focus on a very different class of scenarios.

Scalar fields that have large vacuum expectation values (VEVs) in de Sitter universe (i.e. during inflation) eventually relax to a minimum of the effective potential after the end of inflation. This process often leads to instabilities and subsequent fragmentation of the initially homogeneous scalar field condensate into solitonic lumps and free particles in plasma, which was studied in detail in connection with scalar fields predicted by supersymmetry (SUSY) Kusenko:1997si ; Dine:2003ax . The basic physics of such instabilities is similar to the origin of Jeans instability, which develops due to the attractive nature of gravity. In the case of a scalar field, the attractive self-interactions inside the condensate may be stronger than gravity, while the same interactions in empty space do not generate long-range forces since the scalar mediator is massive. Solitonic lump formation from the scalar field instabilities requires the scalar potential to be shallower than quadratic, ensuring that attractive field self-interactions are present.

The scalar condensate lumps formed in this process are sizable, and their distribution is stochastically generated. Hence, certain clusters of lumps have a sufficient overdensity to become building blocks for primordial black holes Cotner:2016cvr ; Cotner:2017tir ; Cotner:2018vug . Since the source of perturbations generating PBHs on small scales is decoupled from the inflationary perturbations seeding large scale structure, it is not necessary to fine-tune the inflaton potential for the purpose of producing PBHs. Unlike the standard PBH formation scenarios, PBH formation from fragmentation can occur after inflation either before or after the reheating has happened. While relativistic particles and other forms of radiation have energy density scaling of $\rho\propto a^{-4}$, with $a$ being the cosmological scale factor, the population of solitons scales as matter, i.e. $\rho\propto a^{-3}$. Thus, the matter in the form of the fragments can dominate the energy density at some later time. Eventually, the scalar lumps decay, but the temporary matter-dominated era during which the density fluctuations are large, unrelated to the large-scale structure, and unconstrained by the CMB allows for production of PBH. The timeline common to this class of models is illustrated in Fig. 1, where the matter in the form of scalar lumps (oscillons Cotner:2018vug or Q-balls Cotner:2016cvr ; Cotner:2017tir ) dominates the energy density for a short time, during which PBH are produced.

Since PBHs form during a short matter-dominated phase from matter-like lumps, the associated angular momentum leads to PBHs with a higher spin (described by spin parameter $a_{s}$) than expected from the radiation-era production scenarios. We provide a comparative overview of the two types of PBH production scenarios in Table 1.

The perturbations resulting from scalar field fragmentation are a form of isocurvature perturbations. However, unlike the perturbations generated by the inflaton or in models with a spectator curvaton field (e.g. Ando:2017veq ), these perturbations are produced within the horizon after the end of inflation. Therefore, they do not affect the large-scale structure and any growth of perturbations during the temporary matter-dominated era associated with the scalar fragmentation takes place at the scales too small to have any effect on the CMB.

In this work we describe in generality the mechanism of PBH production based on scalar field fragmentation into solitonic lumps. The mechanism can be realized within a broad class of particle physics models. The scalar field that undergoes fragmentation can be either a spectator or the inflaton, and can be associated with a broad range of inflationary models. Previous studies focused on the special cases of a complex scalar field fragmenting into Q-balls Cotner:2016cvr ; Cotner:2017tir and a real scalar field fragmenting into oscillons Cotner:2018vug . We generalize these results and provide a universal analytic description for PBH formation from solitonic energy lumps with an arbitrary number of conserved discrete or continuous quantities, such as field charge. In addition, we elucidate the relation between the model parameters and the sizes of the resulting black holes. In particular, we show that supersymmetric flat directions produce PBH in the range of masses allowed for dark matter.

The paper is organized as follows. The fragmentation of a scalar field is described in Section 2. We then present a general framework for PBH production in Section 3. We summarize in Section 4.

2 Formation and fragmentation of a scalar condensate

Assuming inflation took place in the early universe, followed by reheating and Friedmann–Robertson–Walker expansion, scalar fields with relatively shallow potentials (or relatively small masses) that are present in the theory undergo a non-trivial time evolution. Such fields tend to develop non-zero vacuum expectation values during inflation, followed by a coherent evolution towards the minimum of the effective potential when inflation is over. The latter is a process in which the field may or may not evolve in a spatially homogeneous manner. Depending on the scalar potential, an initially homogeneous field may either remain homogeneous or fragment into lumps. The instability giving rise to the fragmentation occurs when the scalar self-interaction is attractive and some other auxiliary conditions are met. Overall, the fragmentation is a fairly generic phenomenon and is expected to occur in a large variety of models, including supersymmetric theories.

2.1 Spectator fields in de Sitter universe

We assume inflation in the early universe driven by the inflaton field $\Phi_{I}$ with the potential $U_{I}(\Phi_{I})$, which leads to density perturbations consistent with observations. We will not make any additional assumptions about the details of the inflationary model. We shall focus on a different (spectator) scalar field $\phi$ with a potential $U(\phi$) such that it makes a negligible contribution to the total energy density of the universe during inflation and subsequent reheating, that is

[TABLE]

If the scalar potential of the spectator field $U(\phi)$ is relatively shallow, the field develops a non-zero average expectation value when averaged over super-horizon scales Bunch:1978yq ; Linde:1982uu ; Lee:1987qc ; Starobinsky:1994bd ; Dine:2003ax . In de Sitter space, there is a non-zero probability of tunneling from a lower to a higher value of the effective potential Lee:1987qc . Therefore, a field with an initial zero value at the minimum of the effective potential can develop a non-zero value by a quantum tunneling event. The subsequent relaxation of such a field to the potential minimum can occur by a coherent motion. However, if the potential is not very steep, the time-scale associated with this motion is much longer than the rate of tunneling events, implying that the resulting average VEV of $\phi$ remains non-zero.

In the simplest case of a zero potential, $U(\phi)=0$, the field undergoes a random walk forming a Bunch-Davies vacuum: $\langle|\phi|\rangle\propto t^{1/2}$ Bunch:1978yq . For non-zero potentials the expectation value is such that every scalar degree of freedom carries energy density determined by the Hubble parameter or the “Gibbons-Hawking temperature” Bunch:1978yq ; Linde:1982uu ; Lee:1987qc ; Starobinsky:1994bd ; Dine:2003ax :

[TABLE]

where

[TABLE]

and $M_{pl}$ is the reduced Planck mass. This average value of $\langle\phi\rangle$ serves as the initial condition for post-inflationary physics. As seen from Eq. (2), a shallower potential leads to a larger initial VEV. The field maintains this value until the Hubble parameter $H$ decreases and the time-scale for the coherent motion becomes shorter than the Hubble time.

2.2 Scalar field instabilities: heuristic arguments

We will now outline several heuristic arguments to provide an intuitive understanding of the general features of the scalar evolution in the presence of self-interactions.

A scalar field with a non-zero VEV at the end of inflation must relax to the minimum of its effective potential. It is insightful to consider the connection between the scalar field instabilities and the gravitational instability of self-gravitating gas. Namely, the scalar instability is similar to gravitational Jeans instability that is studied in astrophysics and cosmology in connection with the collapse of gas clouds and subsequent star formation. The attractive self-interaction of gravity makes the homogeneous distribution of particles unstable with respect to formation of dense clumps. It is well known that the associated instability occurs on length scales larger than the Jeans length ($\lambda_{J}\sim 1/k_{J}$), or for the wave numbers

[TABLE]

where $\rho$ is the energy density, $c_{s}$ is the speed of sound and $G$ is the gravitational constant.

Analogously, a homogeneous scalar condensate can be unstable with respect to fragmentation into lumps in the presence of attractive self-interactions. In the case of a scalar field, self-interaction determined by the scalar potential can either stymie or enhance the instability. A repulsive self-interaction can act to counter gravity and keep the distribution homogeneous, while an attractive interaction can provide an attractive force that is stronger and more important than gravity, leading to fragmentation. The relevant forces exist inside the condensate and are mediated by the exchange of a scalar field. If the scalar is massive, there may not be any long-range forces in empty space, while the presence of long-range or sufficiently long-range modes inside the condensate could play a more important role than gravity when it comes to stability and fragmentation.

Let us consider a complex scalar field $\phi(x)$ with a Lagrangian

[TABLE]

A spatially homogeneous solution of the equations of motion can be characterized by density $\rho$ and pressure $p$, which represent the diagonal components $T_{0}^{0}$ and $T_{i}^{i}$ (no sum) of the energy-momentum tensor, respectively:

[TABLE]

For a potential $U(\phi)={\rm const}|\phi|^{n}$, the pressure is related to energy density by a simple relation Turner:1983he

[TABLE]

that has a straightforward physics interpretation. A non-interacting scalar field has potential $U(\phi)\propto\phi^{2}$ and $n=2$, so the scalar condensate behaves as a gas of particles at rest with zero pressure. A potential that is effectively rising faster than the second power of the field (i.e. $n>2$) represents a repulsive self-interaction, which generates a positive pressure. Finally, a potential that rises effectively slower than the second power of the field (i.e. $n<2$) describes an attractive interaction, resulting in negative pressure.

In the presence of gravity, a classical pressureless scalar field is subject to gravitational instability on all scales. The reason is that as long as $p=0$ the speed of sound is zero, resulting in $k_{J}=\infty$. However, due to Heisenberg uncertainty principle, quantum effects result in a finite Jeans mass for a non-interacting scalar field Khlopov:1985jw ; Bianchi:1990mha , hence $k_{J}=(32\pi^{2}GU^{\prime\prime}\rho)^{1/4}$. For a repulsive self-interaction, for example in case of $U(\phi)=m^{2}|\phi|^{2}+\lambda|\phi|^{4}$ with $\lambda>0$, the repulsive forces can counter the gravitational attractive force and either reduce the range of the instability or eliminate it completely Khlopov:1985jw .

An attractive self-interaction corresponds to a potential that effectively grows slower than the second power of the field. Such interactions result in negative pressure, which increases the value of $k_{J}$, opening a wider spectral window for the instability. When the self-interaction is stronger than gravity inside the scalar condensate, it is the scalar interaction that plays the leading role in driving the instability.

The above heuristic arguments highlight why the spatially homogeneous solution can be unstable in presence of attractive self-interactions. However, as we will discuss below, the shape of the potential, its symmetries, and the nature of the homogeneous solution play a central role Kusenko:1997si .

2.3 Which instabilities can lead to black hole formation?

When considering PBH formation, one needs to know the end result of the instability. An initially homogeneous condensate can fragment into stable or relatively stable Q-balls Kusenko:1997si or oscillons  Bogolyubsky:1976nx ; Gleiser:1993pt ; Copeland:1995fq ; Fodor:2006zs ; Kasuya:2002zs ; Honda:2001xg ; Amin:2011hj ; Hong:2017ooe , which scale as cosmological matter and can lead to formation of PBHs. Alternatively, the result of the instability may be a highly inhomogeneous state in which there are no stable high-density lumps and, in fact, the higher density regions may exhibit a higher rate of the condensate decay.

At the level of Feynman diagrams, one can think of an attractive scalar self-interaction as an exchange of a scalar mediator via a tri-linear vertex (see Fig. 2). A cubic interaction of the form $U(\phi)\supset\lambda_{3}\phi^{3}$ is attractive, and it can cause an instability. However, the same interaction is responsible for particle-number changing processes $\phi\phi\rightarrow\phi$. Such reactions can take place in the condensate at the rate that is unsuppressed, as long as the phase space occupied by the quanta in the condensate allows the $2\rightarrow 1$ processes to proceed unimpeded. While a larger value of $\lambda_{3}$ is beneficial for instability and fragmentation, the same coupling also controls the decay of the scalar lumps, which can be detrimental for the PBH formation. This obstacle to PBH formation can be circumvented if the particle non-conserving processes are suppressed either by a symmetry that enforces the particle number conservation (the case of Q-balls), or if the scalar condensate fragments into solitons in which the quanta occupy a peculiar phase space suppressing the $2\rightarrow 1$ processes (the case of oscillons).

Q-balls. If a theory possesses a U(1) symmetry, there is a stable Q-ball solution in a scalar potential that allows an attractive interaction Coleman:1985ki . Since odd powers do not respect the global U(1) symmetry $\phi\rightarrow e^{i\theta}\phi$, it was thought that renormalizable theories in 3+1 dimensions do not allow for Q-balls. However, such solitons are allowed in theories that contain more than one field, if some of the fields have a zero U(1) charge Kusenko:1997zq . The requisite attractive interaction can be given, for example, by an effective vertex $U(\phi)\supset\lambda_{\chi\phi\phi}\chi\phi^{\dagger}\phi$ that respects the U(1) symmetry $\phi\rightarrow e^{i\theta}\phi,\chi\rightarrow\chi$ and, therefore, allows for a particle number conservation. Such self-interaction can, indeed, lead to fragmentation of a scalar condensate into Q-balls, which can be stable or long-lived and can merge into PBH.

Oscillons. The longevity of oscillons is supported by a dynamical invariant that exists for a class of potentials for a real scalar field Kasuya:2002zs . These solitons are related to Q-balls in that an oscillon can be described as a Q-ball in a related theory with a complex scalar field and an approximate U(1) symmetry in the non-relativistic limit Mukaida:2014oza ; Ibe:2019vyo . Therefore, fragmentation of a real scalar field into oscillons bears some similarities to the complex field fragmentation into Q-balls discussed below.

2.4 Instability of scalar field with a U(1) symmetry

We will now discuss some general features of the relevant instabilities, which can be studied by means of Floquet analysis, following Ref. Kusenko:1997si . Let us consider a complex scalar field $\phi(x,t)$ that does not dominate the energy density. Generalizations to a real scalar field and the inflaton are analogous. To explore the scalar potential $U(\phi)$ with a minimum at $\phi=0$ and that respects a U(1) symmetry corresponding to a conserved particle number, it is helpful to reparametrize the field as a real scalar field with a time-dependent complex phase: $\phi=\Phi(x,t)e^{i\Theta(x,t)}$. The spatially homogeneous solution to the equations of motion is described by

[TABLE]

The field starts at some initial value, determined by high-scale physics and inflation, and then relaxes to the potential minimum. Let us examine the stability of such a solution with respect to small perturbations. The classical equations of motion in the spherically symmetric metric $ds^{2}=dt^{2}-a(t)^{2}dr^{2}$ are

[TABLE]

where dots denote the time derivatives, and the space coordinates are labeled by the Latin indices that run from 1 to 3. The Hubble parameter $H=(\dot{a}/a)$, where $a(t)$ is the scale factor, is proportional to $t^{2/3}$ or $t^{1/2}$ for the matter or radiation-dominated universe, respectively.

To study the stability of a homogeneous solution, we add small perturbations $\delta\Phi,\delta\Theta\propto e^{S(t)-i\vec{k}\vec{x}}$ to the system and look for growing modes, $({\rm Re}\,\alpha>0)$, where $\alpha=dS/dt$. The wave number $k$ in the comoving frame is red-shifted with respect to the physical wave number $\tilde{k}=k/a(t)$ in the expanding background.

From the equations of motion, one can derive the equations for small perturbations and linearize them by neglecting any terms quadratic in $\delta\Phi,\delta\Theta$:

[TABLE]

The above implies the following dispersion relation, where we define $\omega\equiv\dot{\Theta}$ :

[TABLE]

We want to identify the range of modes $k$ for which the Lyapunov exponents $\alpha(k)$ are positive, leading to an exponentially growing instability. Examining the zeros of $\alpha$, Eq. (16) becomes

[TABLE]

If the potential has a negative second derivative $U^{\prime\prime}<0$ or if the homogeneous solution has a large enough $\omega>U^{\prime\prime}$, then the following inequality is satisfied

[TABLE]

leading to two non-negative values of $k$ that solve Eq. (16), namely $k=0$ and $k=k_{\rm max}\equiv a\sqrt{\omega^{2}-U^{\prime\prime}(\Phi)}$. Therefore, the function $\alpha(k)$ has only two roots and $\alpha(k)$ must have a constant sign for $0<k<k_{\max}$. By examining the sign of $\alpha$ in the limit $k\rightarrow 0^{+}$, one can see that $\alpha(k)>0$ over the entire range $0<k<k_{\max}$. This implies there is a band of unstable modes given by

[TABLE]

If the condition of Eq. (18) is not satisfied, then all the modes are stable. In the case of a real field or a complex field with no time-dependent phase, when $\omega=\dot{\Omega}=0$, the instability discussed above reduces to the well known tachyonic instability:

[TABLE]

We conclude that a broad class of potentials leads to an instability and fragmentation into Q-balls. These potentials include the tachyonic potentials (e.g. the logarithmically growing potentials in gauge-mediated supersymmetry models) as well as potentials with a positive second derivative, which can also lead to a fragmentation into Q-balls if the condensate has a sufficiently high $\dot{\Theta}$, implying a sufficiently high particle number density $n=\dot{\Theta}\Phi^{2}$.

2.5 Solitonic fragments

The instability described above in the linear regime can develop and become non-linear (i.e. $\delta\phi\gtrsim\overline{\phi})$ if it grows fast enough in comparison with the expansion rate of the universe. The end result of this evolution can be a spatially inhomogeneous state consisting of Q-balls Coleman:1985ki or oscillons Bogolyubsky:1976nx ; Gleiser:1993pt ; Copeland:1995fq ; Fodor:2006zs ; Kasuya:2002zs ; Honda:2001xg .

For a complex scalar field with a potential that respects a global U(1) symmetry, the end result of the evolution is a distribution of non-topological solitons, Q-balls, which minimize the energy for a fixed particle number Coleman:1985ki . These can be thought of as coherent bound states of scalar particles.

All supersymmetric extensions of the Standard Model contain Q-balls in their spectrum, made of scalar baryons and leptons predicted by supersymmetry Kusenko:1997zq . They can be entirely stable if the scalar field carrying the baryon number extends along a flat direction of the scalar potential Dvali:1997qv ; Kusenko:1997si . As already noted, the size of the Q-balls or other lumps resulting from fragmentation is determined by the fastest growing mode that becomes non-linear the earliest, denoted $k_{\rm nl}$. Based on numerical simulations111These numerical simulations have been performed for scalar fields that arise in supersymmetric models, and the fragmentation was considered in the context of Affleck-Dine baryogenesis and dark matter production. While the flat directions associated with supersymmetry are well suited for PBH production, both analytical and numerical results are applicable to a broad range of scalar fields. The simulations have been performed in both 1+1 as well as in full 3+1 dimensions and taking the expansion of the Universe into account.  Kasuya:2000wx ; Multamaki:2002hv ; Kusenko:1997si , the size $R_{Q}$ of a typical lump constitutes a few percent of the horizon size:

[TABLE]

where $f_{Q}\sim 10^{-2}-10^{-1}$. Since the resulting soliton number density at fragmentation time is approximately

[TABLE]

we expect the total number of fragments filling the volume to be in the range $\sim 10-10^{6}$. The matter component composed of scalar lumps generated at the time of fragmentation can later come to dominate the energy density. The nature of the fragments depends on the potential.

Oscillons, pseudo-solitonic bound field configurations Bogolyubsky:1976nx ; Gleiser:1993pt ; Copeland:1995fq ; Fodor:2006zs ; Kasuya:2002zs ; Honda:2001xg , can form as a result of a real scalar field fragmentation. The oscillons arise in many well motivated theories with scalar fields, such as models of inflation Amin:2011hj ; Hong:2017ooe , axions Kolb:1993hw or moduli Antusch:2017flz . The oscillons are localized and metastable. The stability of an oscillon is due to an approximate adiabatic invariant Kawasaki:2015vga ; Kasuya:2002zs , which can be described as an approximately conserved particle number in a non-relativistic limit of a related theory with a complex scalar field Mukaida:2014oza ; Ibe:2019vyo . Early universe oscillons have been recently studied in connection with primordial gravity waves Antusch:2016con as well as baryogenesis Lozanov:2014zfa .

The solitonic lumps we are considering are significantly smaller than the size of the horizon at fragmentation time. This allows us to treat them effectively as point particles and renders the effects of gravitational self-interactions for each individual fragment formation as negligible. Hence, gravitational effects become only relevant for a cluster of lumps and the subsequent collapse to a black hole, while the dynamics of formation of individual lumps (e.g. instabilities and fragmentation) are governed solely by the field’s self-interactions not associated with gravity (potential shallower than quadratic). In the presence of gravitational field interactions and curved space-time non-linear field theories can also admit massive macroscopic gravitational soliton configurations, which can be collectively referred to as “soliton stars” (for review see e.g. Ref. Lee:1991ax ; Schunck:2003kk ). They have been studied extensively in the literature under various names and in different specific realizations, including boson and fermion stars (e.g. Ruffini:1969qy ; Lee:1991ax ; Henriques:1989ez ; Kaup:1968zz , oscillatons UrenaLopez:2001tw , axion stars and the related mini-clusters (e.g. Kolb:1993zz ; Braaten:2018nag ). Collapse to a black hole as well as general stability of individual soliton stars has been extensively explored (e.g. Ref. Helfer:2016ljl ; Braaten:2015eeu ). Due to their size, we do not expect a significant amount of such soliton stars to be formed within a single horizon. Hence, in such case, it will be a non-trivial task to achieve large density perturbations due to a cluster of soliton stars within some volume in order to produce PBHs via the mechanism described in this work. We shall thus focus solely on non-gravitationally coupled soliton configurations.

3 Primordial black hole formation from lumps of scalar field

3.1 General framework

We are interested in finding the expected PBH energy density $\braket{\rho_{\text{PBH}}}$ as a function of PBH mass $M$, i.e. the PBH * mass function*, resulting from gravitational collapse of a collection of soliton fragments. Solitons residing in a group (a cluster) whose local density $\rho$ is larger than the soliton background density $\langle\rho\rangle$, when the fractional overdensity fluctuation is

[TABLE]

tends to infall under the action of gravity and could form a black hole.

If some fraction $B(M,V)$ of soliton clusters of mass $M$ residing within an initial volume $V$ will eventually collapse to black holes, then the expected average energy density of resulting BHs is given by

[TABLE]

where $\rho_{\rm M}=M/V$ denotes the cluster energy density and $P(M|V)$ is the probability of finding a cluster of mass $M$ within a given volume $V$. This is the contribution due to a single volume scale as well as corresponding BHs of a particular size. To get the total average energy density, we sum over all volume scales and effectively over PBHs of all sizes, which would contribute to the process. The total resulting PBH energy density is

[TABLE]

The above treatment allows to capture the behavior of soliton clusters with an arbitrary mass-energy distribution $P(M|V)$ and further input assumptions are required to make progress. In a typical situation the soliton cluster mass distribution is related to the number of individual solitons composing the cluster and whose number distribution is specified. Then, $P(M|V)$ can be decomposed as

[TABLE]

where $P(M|N)$ is a distribution of soliton cluster mass $M$ given a soliton number $N$. Here, $P(N|V)$ is the number distribution of individual solitons within $V$, which we will assume to be Poisson. The decomposition of Eq. (26) is employed to treat PBH formation from oscillons.

In principle, other quantities aside the individual soliton number can have an effect on the mass-energy distribution of the soliton clusters. In the case of Q-balls, the additional dependency originates from the conserved Q-ball charge. This leads to a further conditional decomposition of $P(M|V)$ as

[TABLE]

We note that summation in $N$ accounts for soliton number being a discrete quantity, while the integration in $Q$ accounts for charge being continuous.

In the most general case, if the soliton cluster mass-energy depends on an arbitrary number of discrete quantities $Q_{i}$ forming a set $\{Q\}$ and continuous quantities $\xi_{j}$ forming a set $\{\xi\}$, then the respective mass distribution will be given by

[TABLE]

If the quantities $\{\xi\}$ and $\{Q\}$ are distributed independently, then $P(\{Q\},\{\xi\}|V)$ can be separated further as in Eq. (3.1), yielding

[TABLE]

The expressions above can be simplified if we assume that all solitons within a given volume $V$ are identical (denoted “id”). That is, while the quantities in $\{Q\}$ and $\{\xi\}$ can vary independently, the resulting value of each quantity is the same for every soliton. Then,

[TABLE]

where $P(m|\{Q\},\{\xi\})$ is the probability of an individual soliton to have mass $m$ given $\{Q\}$ and $\{\xi\}$ and the soliton cluster mass is determined by a mass function of $m$, $M_{f}(m|\{Q\},\{\xi\})$, which depends on the values of the quantities in $\{Q\}$ and $\{\xi\}$. In the simplest case of $N$ identical solitons of mass $m$ and no other dependencies, i.e. $\{Q\}=N$ and $\{\xi\}=\emptyset$, one obtains

[TABLE]

where we have assumed that distribution of mass for an individual soliton is independent of the number of solitons present, i.e. $P(m|N)=P(m)$.

In contrast, if the mass of solitons follows an identical distribution $P(m)$ but is determined for each individual soliton separately (i.e. solitons are independently and identically distributed, “iid”), then the mass distribution of a soliton cluster composed of $N$ solitons follows

[TABLE]

The difference between the above two assumptions (id vs. iid) can be readily illustrated for the case of soliton masses being normally distributed, with a PDF for single soliton being

[TABLE]

where $m_{0}$ is the mean mass of an individual soliton and $\sigma^{2}$ is the mass variance. The mass distribution of a soliton cluster with $N$ identical solitons then follows

[TABLE]

The resulting expectation value $E$ and variance Var of a cluster mass, compared to single soliton case, are

[TABLE]

While the id vs. iid mean cluster mass grows linearly in the number of solitons, the cluster mass spread is distinct for the two cases. This is just the standard probability theory result for a collection of $N$ iid variables with a Gaussian distribution.

The field fragmentation process is stochastic, as supported by numerous simulation results (e.g. Kusenko:1997si ; Kasuya:2000wx ; Multamaki:2002hv ). Hence, we can infer the soliton probability distribution within a given volume, $P(N|V)$, from statistical fluctuations alone. Given an average number density $\overline{n}$ of uniformly distributed objects, the probability of finding $N$ objects within a volume $V$ follows the usual Poisson distribution

[TABLE]

The Poisson distribution originates from binomial distribution of $N$ independent random events in the large event limit. This choice of $P(N|V)$ constitutes the most general description of fragmentation assuming that events are uncorrelated and further it allows for an analytic treatment. A deviation from this distribution would be highly model dependent and we thus do not consider this possibility further.

The resulting PBH mass function can be readily obtained from Eq. (25) as

[TABLE]

This gives the average energy density that goes into forming PBHs at the time of scalar field fragmentation. Once the black holes are formed and the universe becomes radiation dominated, the PBH mass distribution is relatively unchanged. To connect with observational constraints, the above PBH energy density must be redshifted into the future.

3.2 Conditions for black hole formation from extended objects

Since solitons are extended objects, additional conditions are required to ensure that they can form BHs. If fragmentation takes place at a time $t_{f}$, corresponding to a scale factor of $a_{f}$, solitons can come to matter-dominate the energy density at a later time $t_{\rm Q}$, corresponding to $a_{\rm Q}$, ending with decays resulting in a radiation-dominated era at $t_{R}$, when the scale factor is $a_{\rm R}$. In terms of redshift $z$ and time, this is

[TABLE]

where the factors $r_{f}=t_{Q}/t_{f}$ and $r=t_{R}/t_{Q}$ have been defined for future convenience.

A collection of $N$, mass $M$ solitons, can form a BH with mass $M_{\rm BH}=NM_{s}$ at $a=a_{Q}$ if they are located within a Schwarzschild radius $R_{\rm BH}(N)=2NM_{s}/M_{\rm pl}^{2}$, with $M_{\rm pl}$ being the reduced Planck mass. The solitons must be sufficiently compact to fully fit inside the BH horizon and if one is to treat them particle-like. We require that each individual soliton exceeds the size of a BH with the same mass, but a collection of $N$ solitons can be fully enclosed in a volume of BH with mass $NM_{s}$ without solitons overlapping in space. Assuming all solitons are spherical with identical radius $R_{s}$, these restrictions translate to

[TABLE]

In principle, the r.h.s. condition of Eq. (41) should also reflect the theoretical limit of how dense non-overlapping spheres can be packed into a given volume, which for 3-D Euclidean space constitutes $\sim 75\%$ of the volume (see Kepler’s theorem Hales:2015 ). We do not account for this restriction, which will not drastically alter our conclusions and that for our problem at hand might be further modified by GR-related corrections. Further, while the resulting solitons are identically randomly distributed and won’t have the exactly same radius, from simulations (e.g. Kusenko:1997si ; Kasuya:2000wx ; Multamaki:2002hv ) we expect that their size is strongly concentrated around the mean value, which we take to be the $R_{s}$ above. We also note that, strictly, Eq. (41) is for resulting BHs with negligible rotation (i.e. spin parameter $a_{s}\sim 0$). Since PBHs from matter-dominated formation will typically have a non-negligible spin, the resulting BH radius $R_{\rm BH}=(1+\sqrt{1-a_{s}})GM_{\rm BH}$ can be up to a factor of 2 smaller. This effect will not drastically alter our analysis results and is not included.

Furthermore, particles separated by distances $\sim\,d$ can form a black hole if they are gravitationally bound and do not participate in the Hubble flow. That is, the kinetic energy associated with the Hubble flow velocities $v\sim Hd$ is smaller than the gravitational potential energy:

[TABLE]

Since the solitons are non-overlapping, the separation distances at the time of formation $d_{f}$ cannot be smaller than the soliton radius, $d_{f}>R_{s}$. The expansion of the universe stretches these distances at least by factor $r_{f}^{1/2}$ until they can form a gravitationally bound system during the matter dominated era $d>R_{s}r_{f}^{1/2}$. Hence, condition of Eq. (42) implies

[TABLE]

which for some parameter values can supersede the r.h.s. of condition in Eq. (41).

The above conditions constrain the allowed mass-radius relations for solitons that are to compose BHs and for a given theoretical model they can be translated into restrictions on the allowed range of model parameters. For Q-balls, conditions of Eq. (41) and Eq. (43) can be directly translated into conditions for parameters of the potential through Eq. (3.3), resulting in

[TABLE]

and

[TABLE]

Since for the scenario we will consider below $3\beta-\alpha=0$, the above constraint will be independent of $Q$.

3.3 Primordial black holes from supersymmetric Q-balls

We illustrate the application of the general framework presented in Section 3 by outlining in detail PBH formation from supersymmetric Q-balls, following analysis of Ref. Cotner:2016cvr ; Cotner:2017tir (see Ref. Cotner:2018vug for PBH formation from oscillons).

As mentioned, Q-balls appear in all SUSY extensions of the SM Kusenko:1997zq . For instance, in gauge-mediated supersymmetry breaking (parametrized by the scale $M_{\rm SUSY}$), supersymmetric flat directions are lifted by a potential

[TABLE]

We note that this potential has a negative second derivative at large VEVs, so that the condition of Eq. (18) is automatically satisfied for any value of $\dot{\Omega}$. If the flat direction carries a conserved $U(1)$ number, such as a baryon or lepton number, then the mass of a Q-ball and its profile are determined by the amount of baryonic or leptonic charge $Q=B,L$ as

[TABLE]

More generally, the mass and radius of a Q-ball are given by some mass scale $\Lambda$ and parameters $\alpha$ and $\beta$, with $0<\alpha$ and $\beta<1$, which depend on the shape of the scalar potential

[TABLE]

While Q-balls are stable with respect to the decay into the original scalar field particles, they could decay, for example, into lighter fermions Cohen:1986ct ; Kawasaki:2012gk ; Enqvist:1998xd or if the associated $\mathrm{U}(1)$ symmetry is broken by the higher-dimension operators Kusenko:2005du ; Kasuya:2014ofa ; Cotner:2016dhw ; Kawasaki:2005xc . We parametrize the total model-dependent Q-ball decay width as

[TABLE]

where $\tau_{Q}$ denotes the Q-ball lifetime and contributions of all channels are included. Q-balls have found many applications in particle physics, e.g. they play an important role in Affleck-Dine baryogenesis Dine:2003ax .

As discussed above, supersymmetry provides a strong motivation for considering PBH formed from the scalar field fragmentation. In particular, the electroweak-scale supersymmetry naturally results in the PBH masses for which there are no observational constraints and which can account for all dark matter. While different models of supersymmetry-breaking are possible, for definiteness, let us consider gauge-mediated supersymmetry breaking at the scale $M_{\rm SUSY}\sim 100\ {\rm TeV}$. At the time of fragmentation, when the energy density of the condensate is a fraction $f\sim r_{f}^{-1/2}\lesssim 1$ of the total energy density ($\rho_{\phi}\sim M_{\rm SUSY}^{4}\equiv f\times\rho_{\rm tot}$), the mass inside the horizon is

[TABLE]

If the number of Q-balls per the horizon is $N_{H}\sim 100$ and if a typical PBH has a mass consisting of 10 to 100 Q-ball masses, then the mass of a typical PBH is

[TABLE]

which is consistent with the open window for dark matter in the form of PBH at masses

[TABLE]

The estimates above are crude, and the precise PBH mass function will be calculated below.

It is intriguing that the supersymmetry just above the electroweak scale predicts the masses of PBH consistent with the current bounds for dark matter.

3.3.1   Background energy density of Q-balls

In order to determine the collapse probability of a Q-ball soliton cluster to a BH, we first need to obtain the average background Q-ball density $\braket{\rho_{Q}(t_{f})}$ at fragmentation time $t_{f}$, which enters into the expression for overdensity fluctuations $\delta$ of Eq. (23). This can be done by finding the total Q-ball mass-energy within the horizon $M_{H}=M_{H}(t_{f})$ and then averaging over the horizon volume $V_{H}=V_{H}(t_{f})$, where all the relevant quantities are evaluated at the fragmentation time. Using previous notation, this corresponds to

[TABLE]

which approximates $M_{H}/V_{H}$ in the large $V_{H}$ limit. Here, $M(Q,N)=M_{Q}N=\Lambda|Q|^{\alpha}N$ is the total horizon mass containing $N$ identical Q-balls of mass $M_{Q}$. We assume that the total field charge $Q_{H}$ within the horizon is distributed evenly among all the Q-balls and that this charge is conserved, which we enforce through $P(Q|V_{H})$. As before, we take $P(N|V_{H})$ to be Poisson distribution, centered around $\overline{N}_{H}=\overline{n}V_{H}$ amount of Q-balls present within horizon. Then,

[TABLE]

In the very large soliton number limit, the expression in square brackets of Eq. (3.3.1) approaches $\braket{N^{1-\alpha}}\approx\overline{N}_{H}^{1-\alpha}$. The energy density then simply becomes

[TABLE]

which is the form we will use in further computations below.

3.3.2   Cosmological evolution of energy density

In order to later make connection with observations at present time, we must comment on the cosmological evolution of the energy density and compute the relevant scale factor $a$ after field fragmentation. As a starting point, we assume an initial period of inflation and reheating, resulting in a radiation-dominated epoch with a uniformly distributed charged scalar field condensate as a sub-dominant component of the energy density. At $t_{f}$, the scalar field fragments into Q-balls, which scale as matter and come to dominate the energy density at a later time $t_{Q}$. Primordial black holes are produced during the Q-ball matter-dominated era. As the Q-balls decay at $t_{R}$, radiation comes to dominate again and the resulting PBH density is frozen in and evolves as a non-relativistic matter to present time. The quantities corresponding to these epochs are denoted with associated subscripts.

Initial radiation-dominated era:   At the end of inflation, the inflaton coherently oscillates at the bottom of potential, resulting in a matter-dominated era. At $t_{\rm RH}=1/\Gamma_{I}$ the inflaton decays and reheats the universe, onsetting a radiation-dominated era with temperature of $T_{\rm RH}=0.55g_{\ast}^{-1/4}(\Gamma_{I}M_{\rm pl})^{1/2}$ and radiation energy density

[TABLE]

where $\Gamma_{I}\sim 1/t_{\rm RH}$ is the decay width of the inflaton and $g_{\ast}(T)$ denotes the relevant number of temperature-dependent relativistic degrees of freedom. During this period, using Eq. (56), the radiation redshifts as

[TABLE]

At $t_{f}$, the scalar field condensate fragments into Q-balls, with the energy density given by Eq. (55). The time evolution of Q-balls follows that of decaying non-relativistic matter

[TABLE]

Q-ball matter-dominated era:   As Q-balls evolve as matter, at $t_{Q}$ they will come to dominate the energy density, with $\rho_{R}(t_{Q})=\braket{\rho_{Q}(t_{Q})}$. Since Q-balls decay, the temperature associated with radiation will decrease less slowly compared to expansion alone. Following Scherrer:1985 , the radiation density in the presence of reheating decays is given by

[TABLE]

where $x=\Gamma_{Q}t$, $x_{0}=\Gamma_{Q}t_{Q}$ and $z=(x/x_{0})^{2/3}$. The Q-balls continue to redshift as well as decay, with the energy density of

[TABLE]

With Q-balls decaying, radiation eventually comes to dominate the universe again, defined by $\rho_{R}(t_{R})=\braket{\rho_{Q}(t_{R})}$. At $t_{R}$, the radiation density is given by

[TABLE]

where $\Gamma(n,x_{0},x)=\int_{x_{0}}^{x}dy\,y^{n-1}e^{-y}$ is the generalized incomplete gamma function. Equating this to the Q-ball density $\braket{\rho_{Q}(t_{R})}$, one has

[TABLE]

which can be restated in terms of $r=t_{R}/t_{Q}$ and $r_{Q}=\tau_{Q}/t_{Q}$ as

[TABLE]

We can now solve numerically222We find that Eq. (63) produces more numerically stable results than Eq. (62), while also requiring specification of one less variable. for $r_{Q}$ as a function of $r$. Specifying ($r_{f}=t_{Q}/t_{f}$, $t_{f}$, $r$), allows us to calculate the rest of the parameters $t_{Q}$, $t_{R}$, $\tau_{Q}$ and $\Lambda Q_{H}^{\alpha}$ via

[TABLE]

from which we can calculate all other quantities.

Standard cosmological era:   After Q-ball decay, return to the radiation era onsets the standard cosmology. The additional extended early-time matter-dominated era shifts the time at which the standard cosmological events (BBN, matter-radiation equality or recombination) take place and evolution of PBHs and radiation should be determined based on thermal history, where cosmological events occur at a specific temperature. Then, the ratio of two scale factors is given by

[TABLE]

We consider the evolution from the era at the end of the matter-dominated epoch, with temperature $T_{R}$ (defined at $\rho_{R}(t_{R})$, following Eq. (56)), to the present epoch with temperature $T_{0}=2.7\text{ K}=2.3$ meV. To avoid entropy injection during BBN, we require $T_{R}>T_{\rm BBN}\sim$ MeV. The resulting present-day scale factor relevant for PBH is then given by

[TABLE]

In above, $a(t_{f})=1$ has been implicitly assumed.

3.3.3   PBH mass function

Using the formalism of Section 3, the redshifted PBH mass function is obtained from evaluation of

[TABLE]

The probability that a soliton cluster of mass $M$ will result in a BH is given by Kodama:1986ud 333In our earlier studies Cotner:2016cvr ; Cotner:2017tir , for concreteness, we have assumed that BH collapse probability in matter-dominated era follows the analyses of Ref. Polnarev:1986bi ; Kokubu:2018fxy . This does not fully capture the behavior of sub-horizon isocurvature perturbations present in our scenario. Further, we have previously assumed that overdensities will grow in the matter-dominated era. However, at the linear level of analysis, no overdensity growth is expected. It is a non-trivial task to properly capture analytically the non-linear behavior of the perturbations and we thus omit here this assumption.

[TABLE]

where the initial overdensity fluctuation is given by

[TABLE]

and the horizon mass at the beginning of matter/Q-ball domination within horizon volume $V_{Q}=V_{H}(t_{Q})$ is

[TABLE]

The step function $\theta$ selects regions where the overdensity at Q-ball matter-dominated era is larger than the critical overdensity threshold $\delta_{c}\sim\mathcal{O}(0.1)$ necessary for a black hole collapse. Additional features of a soliton collection such as asphericity, inhomogeneity and angular momentum Kokubu:2018fxy ; Harada:2017fjm ; Harada:2017fjm will affect the collapse and the associated BH formation probability. However, a detailed study of how these considerations will affect our scenario is beyond the scope of this work. We parametrize these effects in Eq. (68) by a phenomenological prefactor $K$. The probability of collapse as given by Eq. (68) is further subject to additional constraints, such as those of Section 3.2, which can be translated into effective limits on the allowed range of values for soliton number $N$ and charge $Q$ that will lead to the formation of black holes.

The probability density of soliton clusters in Eq. (67) can be further decomposed as

[TABLE]

where in the last step we have integrated out the delta-function associated with Q-ball charge as $P(M|N,Q)$ is the only part of the integral that depends on $Q$. The indirect dependence of $B(M,V)$ on $Q$ will be also automatically accounted for. In the above we have assumed that all Q-balls within a cluster are identical and that all clusters have the same overall charge density.

Combining the above, the resulting PBH mass function contains a delta-function restricting volumes $V$

[TABLE]

To ensure that the volume set by the delta-function is within the integration range, we require that

[TABLE]

with

[TABLE]

The lower bound $V_{\rm min}$ sets the smallest volume that could collapse and form a PBH, which we take to be on the order of the volume occupied by a single Q-ball. The upper bound $V_{\rm max}$ is the smaller of either a) the horizon volume at fragmentation (since this sets the max length scale of density perturbations by causality), or b) the largest volume of a cluster (given fixed mass) that will collapse to form PBH. The latter originates from trading the step-function for a limit on $V_{\rm max}$, found from using the definition of $\delta_{0}$ to solve for $V$. This can be intuitively understood as a maximum volume where the total overdensity is larger than the average density within that volume, as increasing it will dilute the density below the threshold.

Using the obtained scale factor from Eq. (66), the present-day differential fraction of dark matter in PBH is given by

[TABLE]

where $r_{f}=t_{Q}/t_{f}$ denotes ratio of time of the beginning of the Q-ball matter-dominated era $t_{Q}$ and the time at fragmentation $t_{f}$.

3.3.4   Conditions for collapse to a black hole

In deriving the final expression for PBH mass-spectrum Eq. (75), the model-specific quantities, energy scale $\Lambda$ of the potential and charge of a single Q-ball $Q$ have been taken as free parameters and only the product $M_{\text{Q-ball}}=\Lambda Q^{\alpha}$ is fixed (see Eq. (64)). The BH collapse restrictions of Section 3.2 will constrain the allowed $\Lambda$ and $Q$ parameter space.

As an illustration, we show how these restrictions affect the parameter space of several concerete models using computations of Ref. Cotner:2016cvr ; Cotner:2017tir and input parameters given in Table 2, with their associated derived quantities shown in Table 3. In particular, by using the definitions of $M$, $M_{H}$ and Eq. (55), Eq. (44) can be recast as a constraint on $N$

[TABLE]

The above is further simplified once we consider the special case of Q-balls from gauge-mediated SUSY, with $\alpha=3/4$ and $\beta=1/4$, which sets $(\alpha-\beta)/\alpha=1$. The other BH formation condition of Eq. (45) can be viewed as simply a restriction on $\Lambda$. The result of imposing these constraints is shown in Figure 3, where each line corresponds to model A, B, C depicted on the constraints plot on the left.

4 Discussion and conclusion

In this work we have developed a general analytic description of PBH formation from scalar field fragmentation in the early universe.

Formation of PBHs from scalar field dynamics appears to be very common. In addition to the Higgs field, whose existence has been decisively confirmed, theories beyond the Standard Model generally predict appearance of multiple scalar fields. Any one of these fields can develop a large expectation value and undergo fragmentation under some very mild constraints on the shape of the potential. More specifically, the potential should be shallow enough to allow for a large VEV at the end of inflation, in accordance with Eq. (2), and the second derivative of the potential should be small or negative, in accordance with Eq. (18). A good example is supersymmetry. All supersymmetric generalizations of the Standard Model contain multiple “flat directions” along which the potential is zero in the limit of exact supersymmetry and which are lifted by supersymmetry breaking terms Gherghetta:1995dv . The soft supersymmetry breaking terms usually result in a flat direction that grows slowly at large VEV, so that it develops a large VEV at the end of inflation. In the case of gauge-mediated supersymmetry breaking, the condition of Eq. (18) is also satisfied at any scale above the supersymmetry breaking scale. Therefore, such flat directions are well suited for producing PBHs.

Given the preponderance of scalar fields with the requisite properties in theories beyond the Standard Model, PBH formation appears to be a fairly generic phenomenon. This makes PBHs an appealing candidate for dark matter in theories with supersymmetry or other theories involving scalar fields.

Acknowledgements.

The work of E.C., A.K. and V.T. was supported by the U.S. Department of Energy (DOE) Grant No. DE-SC0009937. The work of M.S. was supported by the Ministry of Education, Culture, Sports, Science and Technology of Japan (MEXT) and by Japan Society for the Promotion of Science (JSPS) KAKENHI Grants No. 15H05888 and No. 15K21733. A.K. and M.S. were also supported by the World Premier International Research Center Initiative (WPI), MEXT, Japan. Part of this work was performed at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.