Early formation of supermassive black holes via dark matter self-interactions

TL;DR

This paper investigates how a small fraction of self-interacting dark matter can lead to the early formation of supermassive black holes at high redshifts using N-body simulations, challenging previous fluid-based models.

Contribution

It introduces N-body gravitational simulations to study dark matter self-interactions and their role in early supermassive black hole formation, providing a more detailed analysis than prior fluid models.

Findings

Self-interacting dark matter can seed supermassive black holes at high redshifts.

A specific parameter space allows for early black hole formation consistent with observations.

Both elastic and dissipative scattering processes influence black hole growth.

Abstract

The existence of supermassive black holes at high redshifts ($z\sim7$) is difficult to accommodate in standard astrophysical scenarios. It has been shown that dark matter models with a subdominant self-interacting component are able to produce early seeds for supermassive black holes through the gravothermal catastrophe. Previous studies used a fluid equation approach, requiring some limiting assumptions. Here we reconsider the problem using $N$-body gravitational simulations starting from the formation of the initial dark matter halo. We consider both elastic and dissipative scattering, and elucidate the interplay between the dark matter microphysics and subsequent accretion of the black hole needed to match the properties of observed high redshift supermassive black holes. We find a region of parameter space in which a small component of self-interacting dark matter can produce the observed high redshift supermassive black holes.

Tables

Table 1. Table 1 : A summary of the three SIDM models considered in this work. The last column shows the results of the simulations for the approximate value of the fraction λ 𝜆 \lambda of total SIDM mass that forms the SMBH within each model (see figures  1 and  3 ).

Model	Description	$\lambda =M_{\mathrm{SMBH}}/M_{\mathrm{SIDM}}$
$A$	Elastic scattering	$6\times {10}^{-3}$
$B$	Inelastic above cut-off $v_c=\sqrt{2\Delta E/m}$, elastic below $v_c$	$1\times {10}^{-3}$
$C$	Totally inelastic scattering, elastic scattering once bound state is formed	$6\times {10}^{-4}$

Table 2. Table 2 : The redshifts and masses of the three highest- z 𝑧 z SMBHs, which we use to compare our results to observations.

Galaxy	Redshift	$M_{\mathrm{SMBH}}$
J1342+0928[1]	7.54	$7.8\times {10}^8$
J1120+0641[2]	7.09	$2.0\times {10}^9$
J2348–3054[3]	6.89	$2.1\times {10}^9$

Equations

$$M(t)=M_0\exp \left(\frac{1-{ϵ}_r}{{ϵ}_r}\frac{t}{0.45\mathrm{Gyr}}\right).$$

$$t_d=r_c/v_{\mathrm{rms}}\ll t_r,$$

$$\frac{\sigma }{m}\lesssim 0.2-1{\mathrm{cm}}^2/\mathrm{g}.$$

$$\rho (r)=\frac{{\rho }_s}{{\left(\frac{r}{R_s}\right)}^{\gamma }{\left(1+\frac{r}{R_s}\right)}^{3-\gamma }},$$

$$t_r=\frac{m}{af\sigma {\rho }_s{v}_s},$$

$$v_s=\sqrt{4\pi G{\rho }_s}{R}_s,$$

$$t_r=0.28\mathrm{Myr}\left(\frac{1\mathrm{c}{\mathrm{m}}^2/g}{f\sigma /m}\right){\left(\frac{10^{10}{M}_{\odot }/{\mathrm{kpc}}^3}{{\rho }_s}\right)}^{3/2}\left(\frac{1\mathrm{kpc}}{R_s}\right).$$

$$\frac{\sigma f}{m}\lesssim \frac{1}{{\rho }_s{R}_s}.$$

$$R_s=11.1\mathrm{kpc},{\rho }_s=1.49\times 10^6{\mathrm{M}}_{\odot }{\mathrm{kpc}}^{-3}$$

$$\rho(r)=\left\{\begin{array}[]{cl}\rho_{i}>\rho_{\rm crit},&r<r_{i}\\ \rho_{o}<\rho_{\rm crit},&r_{o}>r>r_{i}\\ \rho_{\rm crit},&r>r_{o}\end{array}\right.$$

$$v_f=\{\begin{array}{ll}{v}_i& v_i<v_c\\ \sqrt{v_i^2-2\Delta E/m}& v_i>v_c.\end{array}$$

$$v_{\mathrm{circ},\max }\cong 0.244\sqrt{G{M}_0/R_s}$$

$${\rho }_s$$

$$R_s$$

$$v_{\mathrm{circ},\max }=506km/s$$

$$M_{\mathrm{SMBH}}=M_{\mathrm{seed}}\exp \left(\frac{t_{\mathrm{obs}}-t_{\mathrm{col}}}{50\mathrm{Myr}}\right).$$

$$M_{\mathrm{seed}}\approx \lambda f{M}_0$$

$$\frac{\sigma }{m}\sim 10^8\mathrm{c}{\mathrm{m}}^2/g!$$

$$t_{dd}=\frac{\mathcal{E}}{dP/dV},$$

$$T_{\mathrm{vir}}=\frac{1}{5}\frac{G_N{M}_{\mathrm{vir}}m}{R_{\mathrm{vir}}}=5.9\times 10^{-8}m.$$

$$t\approx 6\times 10^3\mathrm{Gyr}\left(\frac{0.1}{f}\right)\left(\frac{1\mathrm{c}{\mathrm{m}}^2/g}{\sigma /m}\right)\left(\frac{10^{-7}}{\Delta E/m}\right)$$

Full text

1 Introduction

Supermassive black holes (SMBHs) are now known to be ubiquitous in the centers of Milky way-like and larger galaxies. Although our own galaxy’s SMBH is quiescent, those in active galactic nuclei (quasars) are highly luminous due to radiation from accretion, outshining their entire host galaxy. In recent years, quasars containing SMBHs with masses of order $10^{9}\,M_{\odot}$ have been discovered at redshifts of up to 7.5 [1, 2, 3]. In standard scenarios for structure formation, it is difficult to account for these large masses at such early times, since the progenitors must start out significantly lighter and only acquire their observed masses through accretion. The rate of accretion is bounded by the Eddington limit, which is the maximum allowed by the balance of gravitational force versus radiative pressure. This restricts the rate of growth to an e-folding time of order $50\,\rm{Myr}$ [4].

Alternative astrophysical mechanisms have been proposed for producing early SMBHs, that typically rely upon boosting the mass of the progenitor to order $\gtrsim 100\,\rm{M}_{\odot}$, so that less accretion time is needed. These include early Population III stars, collisions of stellar-mass black holes and stars in stellar clusters to form black holes with mass $\sim 10^{3}$-$10^{4}\,\rm{M}_{\odot}$, or the direct collapse of low metallicity gas clouds into black holes. For a review of these mechanisms, see ref. [5].

The mass of an accreting black hole as a function of time is given by [4, 5]

[TABLE]

The radiative efficiency is typically taken to be $\epsilon_{r}\approx 0.1$[5]. Then a black hole seed of $M_{0}=10^{2}\,\rm{M}_{\odot}$ would take at least 0.81 Gyr to develop into a $10^{9}\,\rm{M}_{\odot}$ SMBH, whereas a seed with $M_{0}=10^{5}\,\rm{M}_{\odot}$ would take 0.46 Gyr. The age of the universe at $z=7$ is approximately $0.76\,\rm{Gyr}$. This means the seeds must either form very early, e.g., $z=13.5$ in the case of $M_{0}=10^{5}\,{\rm M}_{\odot}$, or be very large, presenting a challenge even for the above mechanisms.

An alternative mechanism is the gravothermal collapse of a self-interacting dark matter (SIDM) halo, as shown in ref. [6], hereafter called PSS14. Gravothermal collapse is the process believed to be the origin of globular clusters, through gravitational interactions that eject more energetic stars, allowing the gravitationally bound system to contract [7]. Such systems have negative specific heat, and the process can run away unless halted by some interaction that prevents further outflow of energy. In the case of globular clusters, formation of binary systems may halt runaway collapse.

Self-interactions of dark matter (DM) can cause the analogous process in DM halos. In this case, there need not be anything that halts the collapse, which results in a black hole. Several early studies of halo formation with SIDM considered this process [8, 9, 10], in the context of using SIDM to solve the core-cusp problem of halo density profiles, rather than trying to explain SMBH formation. Refs. [9, 10] showed that, with proper cosmological boundary conditions applied to the halo, gravothermal collapse would not occur within a Hubble time unless the cross section per DM mass is much larger than that required to match observations of halo profiles, or allowed by constraints from the Bullet Cluster [11, 12], $\sigma/m\sim 1\,\rm{cm}^{2}/\rm{g}$, where $\sigma$ is the elastic scattering cross section and $m$ is the DM mass.

Nevertheless, a subdominant component of strongly interacting DM could still initiate collapse of SMBH seeds while remaining consistent with such bounds, as was first claimed by PSS14, in a study limited to the effects of elastic scattering. More recently refs. [13, 14] investigated this general idea within the framework of mirror dark matter, assuming a large fraction $f\sim 0.2$ of dissipative SIDM. However the mechanism of collapse explored in these works is not the gravothermal catastrophe, but rather a modified version of ordinary SMBH formation, accelerated by lowering the temperature of the dark sector.

In PSS14, the gravothermal collapse was modeled using a set of fluid equations for spherically symmetric distributions of mass, temperature, velocity dispersion and radiated heat. To implement the fluid approach with two DM components, it was necessary for ref. PSS14 to make some simplifying assumptions: first that the initial density for the dominant component followed the usual Navarro-Frenk-White (NFW) profile [15], despite the possible influence of the SIDM component, and second that during the subsequent evolution the two densities should maintain the same profile shape, apart from the different normalizations. One might question whether these assumptions are really innocuous as regards the main features of gravothermal collapse, and to what extent they are borne out in a more exact treatment.

To overcome the limitations of the fluid approach, in this work we reconsider the problem by simulating the gravothermal collapse of a partially SIDM halo using an $N$-body code, initially developed in ref. [10]. We aim for a generic, model-independent treatment, exploring the effects of both elastic and dissipative scattering for the production of SMBHs. Our simplified models of dissipative interactions are designed to mimic energy loss through excitation followed by emission of dark radiation, or the formation of DM bound states.

In section 2, we review the process of gravothermal collapse, introduce the framework of two-component dark matter and summarize the previous results of ref. [6] (hereafter referred to as PSS14) on SMBH formation from elastically scattering DM. In section 3 we describe our $N$-body simulation methodology and present the results of simulations for an elastically scattering subdominant DM component. We show that it is not consistent to assume an initial NFW profile, and that one must instead simulate the full halo formation process. Moreover we show that elastic scattering cannot produce early SMBHs unless the cross section is large, $\sigma/m\gtrsim 10^{3}\,\rm{cm}^{2}/\rm{g}$. In section 4 we turn our attention to two simplified models of dissipative DM, which greatly speeds up the process of collapse, allowing smaller $\sigma/m$ to explain high-redshift SMBHs. In section 5 we combine these results with a model of subsequent accretion to illustrate a range of possible working parameters in the three classes of interactions considered, comparing to the properties of three observed high-redshift SMBHs. We briefly consider the possible formation of black holes in smaller systems, namely dwarf galaxies. Discussion of these results is given in sect. 6 and conclusions in sect. 7.

2 Gravothermal collapse and the gravothermal catastrophe

Gravothermal collapse can occur when heat and matter are transferred out of a virialized, gravitationally bound system of point masses. The virial theorem states that $U=-2T$, where $U$ is potential and $T$ kinetic energy, so that the total energy is $E=U+T=-T$. Such systems therefore have a negative specific heat: when energy is added they become less strongly gravitationally bound (and therefore the kinetic energy, or temperature, decreases), and when energy is removed they become more strongly bound, increasing the temperature.

In a halo with a negative radial temperature gradient, heat and mass will flow radially outward as it evolves towards equilibrium. This causes the inner part of the halo to shrink and further increase in temperature. If the specific heat of the outer halo is smaller than that of the inner, eventually the two regions reach equilibrium and the inner halo stops contracting. If it is larger, the process instead continues in a runaway fashion known as the gravothermal catastrophe [7]. Collapse occurs on a timescale related to the relaxation time $t_{r}$, the average time between collisions for a particle in the halo.

During the contraction, particles may eventually reach relativistic speeds and form a black hole through the radial instability. This occurs on a dynamical timescale,

[TABLE]

where $r_{c}$ is the core radius and $v_{\rm rms}$ the core r.m.s. speed. Once the core reaches relativistic speeds it very quickly collapses into a black hole [16].

This process requires the conduction of heat, which can happen through elastic scattering. A classic example is globular clusters, where heat is transferred by the gravitational interactions of stars, in particular when a higher-energy star is scatterered outward to a larger radius while the lower-energy star falls inward toward the center of the halo (increasing its kinetic energy in the process). In contrast, the DM particles in a cold dark matter (CDM) halo are typically not massive enough for gravitational self-scattering to lead to gravothermal collapse. But the nongravitational self-interactions of SIDM can be much stronger, as we discuss next.

2.1 Self-interacting dark matter

While standard CDM is defined to be collisionless, self-interacting models have garnered much interest in recent years. DM scattering with cross sections per DM mass of order $\sigma/m\sim 1.0\,\rm{cm}^{2}/\rm{g}$ have been shown to ameliorate several problems in CDM small scale structure predictions, including the cusp/core and missing satellite problems [17, 18, 19, 20]. The former refers to the tendency of CDM simulations to produce ‘cuspy’ halos whose densities diverge at small radii [21, 15, 22, 23, 24, 25], in contrast to observations of dwarf and low surface brightness galaxies that indicate a flattening density profile at small radii (cored) [26, 27, 28, 29, 30, 31, 32]. The latter refers to the observation that CDM, while correctly predicting large scale structure and the number and distribution of large halos, predicts far more small satellite halos than are observed [33].

The required cross section for SIDM to solve the small-scale structure problems is of the same order as the upper bound coming from observations of the Bullet Cluster [11, 12] and other colliding galaxy clusters [34], the inner density profile of the Draco dwarf spheroidal galaxy [35], and brightest cluster galaxy offsets [36]. These studies give limits in the range

[TABLE]

However since larger values are needed for gravothermal collapse at early times [37, 10], we are motivated to consider models with two components of DM, that make it possible to evade (2.2), by making the strongly self-interacting component sufficiently subdominant.

2.2 Two-component dark matter

The Bullet Cluster bound (2.2) assumes that all the DM has the same self-interaction cross section, but if DM consists of two (or more) species, the smaller component could have a much larger value of $\sigma/m$. Stemming from observational uncertainities, it is estimated that the colliding DM subcluster could have lost as much as $23\%$ of its mass in the collision [12]. One could then imagine that a fraction of very strongly self-interacting DM as large as $f\sim 0.23$ is allowed. On the other hand it is possible that the allowed fraction is a function of $\sigma/m$; no explicit study of this question, which is outside of the scope of the present work, has so far been done. We will assume that $f$ as large as $0.1$ is allowed, regardless of how large $\sigma/m$ is.

A stronger, complementary bound of $f<0.05$ arises if the DM is significantly coupled to dark radiation, which could lead to dark acoustic oscillations in the matter power spectrum for large scale structure [38]. This however is more model-dependent and can be evaded if dark radiation is absent or suppressed.

In PSS14, a two-component scenario is investigated using a fluid approach, starting from an initial NFW profile and evolving it according to the gravothermal fluid equations. A generalized NFW profile can be defined as:

[TABLE]

with $R_{s}$ the scale radius and $\rho_{s}$ the scale density. The parameter $\gamma$ controls the extent to which the profile is cuspy or cored, with $\gamma=1$ corresponding to the original NFW profile.

The results from the fluid formalism are given in terms of the relaxation time,

[TABLE]

where $a=4/\sqrt{\pi}$ for hard-sphere interactions and $v_{s}$ is the velocity dispersion at the characteristic radius,

[TABLE]

For reference, we will ultimately be interested in halos with mass $\sim 10^{12}\,M_{\odot}$ and NFW parameters $\rho_{s}\sim 10^{10}M_{\odot}/{\rm kpc}^{3}$, $R_{s}\sim 1\,$kpc, leading to $v_{s}=2300\,$km/s and a relaxation time of

[TABLE]

The initial choice of an NFW profile is justified so long as the halo is optically thin at its scale radius,

[TABLE]

This follows from demanding that the relaxation time (2.4) is greater than the dynamical timescale for the halo, $R_{s}/v_{s}$ (analogous to that for the core, eq. (2.1)), ensuring that the initial halo structure is not strongly perturbed by the SIDM component.111 Although SIDM has been shown to result in the formation of a core over time, it is argued in PSS14 that if the relaxation time (the average time between scatterings of a typical particle) is much greater than the dynamical time (the timescale of halo formation), then the average SIDM particle has not scattered at all during formation, and therefore the resulting initial profile should not be far from that of collisionless dark matter. Eq. 2.7 implies that the optical depth of the halo to DM self-interactions is larger than the halo size. PSS14 finds that

• •

The gravothermal catastrophe occurs (and therefore the SMBH forms) after approximately $450\,t_{r}$ regardless of cross section or SIDM fraction $f$. Therefore, the time taken depends only on the combination $\sigma f$.

• •

The SMBH contains $2.5\%$ of the SIDM component. Therefore for a halo with mass $M_{0}$, $M_{\rm SMBH}=0.025fM_{0}$.

• •

There is a region of parameter space in which SMBHs of the correct size may form early enough to accommodate observations ($z=7$).

In the following we will obtain different results: the gravothermal catastrophe occurs after approximately $480\,f^{-2}\,t_{r}$, greatly increasing the time until collapse for halos with a small SIDM fraction, and the SMBH contains a smaller fraction of the total SIDM component, $M_{\rm SMBH}/M_{\rm SIDM}\approx 0.6\%$. Due to the additional dependence on $f$ of the time of collapse, we will find that although there is still a region of parameter space in which SMBHs of the correct size form by $z=7$, the scattering cross sections required are much larger, unless dissipative interactions are introduced. For these larger elastic cross sections, the consistency requirement (2.7) is no longer satisfied, invalidating the assumption of an initial NFW halo.

3 $N$-body simulations of elastically scattering two-component

dark matter

Both $N$-body and hydrodynamical simulations are frequently used to study the collapse of a DM halo. The former have the disadvantage of being quite computationally expensive, as the gravitational potential must be calculated for a large number ($N\sim 500000$ in our case) of particles, which must then be individually evolved forward in time. Scattering probabilities between neighbouring particles must be calculated, along with the resulting velocities if a scattering does occur [10].

Hydrodynamical simulations instead discretize space into a series of radial shells, keeping track of the amount of DM in each shell. This formulates the problem as a set of coupled partial differential equations. When key constants have been correctly calibrated, it can reproduce the results of $N$-body simulations [37, 10]. In this formalism it is difficult to accommodate two-component DM, which is crucial to the formation of high redshift SMBHs. Each DM component requires its own set of shells since the self-interactions differ between the two, but when computing the gravitational potential one would have to interpolate between the shells. Errors in interpolation grow quickly between successive timesteps, making this approach impractical. To circumvent these difficulties, ref. PSS14 applied the hydrodynamical simulation to the SIDM component only, while assuming a gravitational potential consistent with an NFW profile, i.e., the SIDM component does not significantly affect the overall gravitational potential or distribution of CDM. The validity of this assumption is not obvious, motivating our use of $N$-body simulations that are not limited in this way.

3.1 Simulation of gravothermal collapse from an initial NFW halo

As a first step we employed the GADGET $N$-body simulation code [39, 40] to simulate the gravothermal collapse of an initial NFW halo. The main motivation for doing so is to be able to cleanly compare our results with those of PSS14, which used this as an initial condition. In sections 3.2 and 5 we will drop this simplifying assumption and consider formation of the halo starting from a primordial overdensity. GADGET is capable of simulating both noninteracting DM and baryonic gas. Baryonic simulations are much more computationally intensive; because of limited computer time we consider only DM. For ref. [10], GADGET was modified to include DM self-interactions between nearest neighbor particles, and the modified code is available online [41]. We further developed it to allow for dissipative (in addition to elastic) scattering of a subdominant DM component.

To test the code we first considered a single DM component with hard-sphere scattering, having a velocity independent cross section of $\sigma/m=38\,\rm{cm}^{2}/\rm{g}$, to facilitate comparison with previous work [10, 6] that used this value. The initial conditions are that of an isolated NFW halo, as used in ref. [10], which has a total mass $M_{0}=10^{11}\,\rm{M}_{\odot}$, and NFW parameters

[TABLE]

and a maximum halo radius $R_{\rm max}=100\,R_{s}$, at which we place a reflective boundary, reversing the radial velocity of particles which exceed this value. This is chosen to be a sufficiently large cut-off that it has no effect on the dynamics and evolution of the inner halo. From eq. (2.6), the relaxation time is $t_{r}=0.37/f\,$Gyr, which is too long to allow for SMBH formation by $z\sim 7$, for realistic values of $f$. We will consider more promising examples later, in section 5.

As the halo evolves, mass flows inward as expected for gravothermal collapse, until the central density begins to very rapidly increase and causes the timestep $\Delta t$ to approach zero. This occurs because $\Delta t$ goes inversely to the density in the modified code, $\Delta t\sim 1/\rho(r)$, and $\rho(0)$ diverges as the core collapses. For practical purposes, we identify the time at which $\Delta t$ falls to $10^{-5}$ of its initial value as marking the onset of the gravothermal catastrophe, and formation of the black hole seed. At this moment, the inner part of the density profile increases quite suddenly, following a long period of slow evolution. The mass in the central region quickly contracts, leading to a flattening of $M(r)$, the mass enclosed within radius $r$, shown in figure 1). These qualitative observations are consistent with the results of hydrodynamical treatments, where the halo shows very little change over most of its history, followed by a sudden contraction [6, 42].

Our results roughly agree with those of PSS14 for the limiting case of single-component SIDM, $f=1$, with the gravothermal catastrophe occurring after approximately 480 relaxation times (close to their result of $\sim 450\,t_{r}$), as can be seen in figure 1. But for smaller values $f=0.5$ and $f=0.1$, with the combination $f\sigma$ held constant, we find that the gravothermal catastrophe occurs after 1980 or 49000 relaxation times respectively. The results of these two simulations are shown in figure 1. The dependence upon $f$ has a simple form, expressed by the empirical observation that if $f^{3}\sigma$ is held fixed, the time of SMBH formation remains nearly constant.222To achieve greater numerical accuracy for small values of $f$, which would have large statistical fluctuations if the number of SIDM particles was simply reduced, we simulate the normal and SIDM components using equal numbers of particles, but with the SIDM mass adjusted so that the total SIDM mass if only a fraction $f$ of the total DM mass, and $\sigma/m$ is also rescaled accordingly. The $N$-body code is designed to treat these configurations as being physically equivalent.

We thus find that the time of collapse does not simply scale with the relaxation time (2.4), but rather as $1/(f^{3}\sigma)$. This is at first surprising, since one would naively expect that the scattering rate of the SIDM component, proportional to $f^{2}\sigma$ (also at variance with the findings of PSS14), should control heat conduction through the halo. But this heat takes the form of kinetic energy of the SIDM particles, which also scales with their total mass, bringing an additional factor of $f$.333Another way of understanding the additional factor of $f$ could be that gravothermal collapse proceeds through the formation of a shrinking core. If the SIDM is only a small fraction of the overall halo, the core cannot become as massive; its mass scales as $f$. This weakens the gravitational potential of the core proportionally to $f$, slowing its growth and resulting in a total proportionality of $f^{3}$.

The final fraction of the SIDM mass that becomes part of the supermassive black hole is $M_{SMBH}/M_{SIDM}\sim 0.6\%$. This can be seen in figure 1, where the interior mass eventually levels off at small radii, showing that a fixed amount of the SIDM has collapsed to a central region smaller than our minimum resolvable radius.444 The gravitational smoothing, which roughly corresponds to the minimum resolvable radius, is taken to be $0.01{R_{s}}$ for $f=1$ or $f=0.5$, and $0.06{R_{s}}$ for $f=0.1$ to compensate for the much greater computational time required at small $f$.   The SMBH mass, $M_{\rm SMBH}$, is defined as the mass inside this radius at the time of its formation. The fraction of the SIDM that forms the SMBH is independent of $f$. This $f$-independence agrees with the results of PSS14, except that the final value of $M_{\rm SMBH}$ is smaller than their estimate of $2.5\%\times M_{0}$.

Combining these results, we can compare to the limit of $f\sigma/m\geq 0.336\,\rm{cm}^{2}/\rm{g}$ advocated by PSS14 to explain observations of high redshift SMBHs. This has some overlap with the constraint from eq. 2.7, that implies $f\sigma/m<0.425\,\rm{cm}^{2}/\rm{g}$. Our numerical values scale with the relative number of relaxation times before collapse, $480/450=1.07$, but more importantly, our required value for SMBH formation scales as $f^{3}\sigma$, in contrast to the optical depth bound which goes as $f\sigma$. Since $f\lesssim 0.2$ from the Bullet Cluster constraint, there is no longer any overlap between the two inequalities. Hence the assumption of an initial NFW halo with common shape for both the CDM and SIDM components cannot be justified, since the SIDM scatterings could alter both distributions. This motivates our subsequent investigation, where we model the collapse of the halo to determine the impact of violating (2.7) on the initial halo profile.

3.2 Halo formation in a two-component universe

Since the assumption of an initial NFW profile may not be justified, we use GADGET to simulate the formation of a two-component halo using a simple spherical collapse model [43]. An initial spherically symmetric overdensity in the early universe is given by

[TABLE]

Well outside the overdense region, the universe behaves as a flat expanding universe, whereas inside it acts like a closed universe that undergoes expansion to a maximum local scale factor. The density contrast at the time of maximum expansion is ${\rho}/{\rho_{\rm crit}}=5.55$, after which the overdensity begins to collapse.

We simulate these conditions by implementing periodic boundary conditions within a cube of length $L=(2000\,{\rm kpc})/(1+z)$ on each side. Within the cube is a spherical region of uniform density with $r_{i}=(372\,{\rm kpc})/(1+z)$ and $\rho=5.55\,\rho_{\rm crit}$. Outside the sphere, the density is chosen such that the total average density within the cube is $\rho_{\rm crit}$. Due to the periodic boundary conditions, far from the overdense region the universe is effectively flat. The size of the cube and overdensity are chosen such that the latter contains $10^{11}\,\rm{M}_{\odot}$ of DM, facilitating comparison with our prior simulations, that used the same halo mass. We begin the simulation at $z=63$,555 This value is sufficiently early that the halo virializes by $z\sim 15$. Other simulations were done beginning at redshifts of $z=40$, $30$, and $20$. The results are shown in figure 2 to be largely insensitive to the choice of starting redshift. and the initial condition file is constructed using the GADGET initial condition generator [39, 40]. We fit the results to a generalized NFW profile using nonlinear least squares, minimizing over the three halo parameters $\rho_{s}$, $\gamma$, and $R_{s}$. This ansatz is flexible enough to give good fits to our numerical profiles.

The simulation is allowed to continue until $z=15$, by which point the halo will have virialized into an NFW profile. This expectation is borne out by the Milli-Millennium database [44], derived from Millennium Simulation [45] structure formation results for $\Lambda$CDM universes. The largest halo in the dataset at $z\sim 7$ has total mass $M_{0}\gtrsim 10^{12}\,\rm{M}_{\odot}$, and formed at $z\sim 15$. We therefore expect that smaller halos will also have virialized by $z=15$. The results for both CDM and the two-component model are shown in figure 2. In the CDM-only simulation, the DM halo collapses into a NFW profile with $\gamma=1.4$ (see equation 2.3) by $z=15$.

We then performed a two-component simulation with $\sigma/m=380$ cm${}^{2}/$g and $f=0.1$. The scaling law for the time of SMBH formation found above, $t\sim m/(f^{3}\sigma)$, shows that this is nearly the minimum value expected to produce a SMBH by $z=7$, given our choice of halo parameters. Fig. 2(b) shows that the CDM component again collapses into a NFW profile by $z=15$, but the influence of the SIDM leads to a less cuspy profile for the CDM with $\gamma=1.0$. The SIDM component itself is far more cored, and is poorly fit by an NFW profile. Hence for the interesting region of parameter space where $f^{3}\sigma/m\gtrsim 1\,\rm{cm}^{2}/g$, the full collapse of the halo must be simulated, rather than assuming an NFW profile. Given that the two components evolve very differently from each other, the hydrodynamical approach may not be well suited to modelling the gravothermal collapse of a two-component DM halo. A proper treatment would require separate sets of mass shells for the two components, not implemented in PSS14.

4 Dissipative dark matter

We have found that large elastic cross sections $\sigma/m\gg 1\,\rm{cm}^{2}/g$ are required for early SMBH formation, but one expects that gravothermal collapse could be accelerated by instead using dissipative (inelastic) scattering. Such processes can greatly increase the heat flow from the inner halo to the outer, hastening the collapse of the DM halo, for example through the emission of dark radiation. Ref. [13] showed that a subdominant mirror sector could effectively seed SMBHs during structure formation.666Upper limits on dissipative scattering were obtained by ref. [42], in the context of a single component of DM. There is one important caveat: if the dark radiation exerts a significant pressure on the collapsing halo, it can slow or even halt the collapse. In the present work we circumvent this potential issue, by assuming that any radiation or light particles produced during inelastic collisions are free to exit the halo: the optical thickness is larger than the halo size. In this section we continue to use an initial NFW profile for purposes of comparison with previous work using hydrodynamical equations [6]; our final results in section 5 do not rely upon this simplification.

4.1 Dissipative dark matter models

In the interests of making a model-independent analysis, we consider two simplified models of inelastic scattering, that could plausibly capture the essential features of more realistic models. We will refer to them as models $B$ and $C$, with $A$ denoting simple elastic scattering.

In Model $B$, the SIDM loses a fixed quantity of kinetic energy in each scattering event, if sufficient energy is available. This can approximate the effect of creating an excited DM state, that subsequently decays by of radiation or a light particle. Such a transition occurs in multi-state DM models [46], and dark atom models [47], in which collisions between the dark atoms could result in hyperfine excited states with fast radiative decays.777If the temperature of the dark sector is high enough, dark atom-dark electron collisions could also lead to excitation of the dark atoms. While our toy model does not capture this effect, it has been explored in Refs. [48, 49, 50], which also discuss how in combination with other dissipation mechanisms it can lead to collapsed halo substructures. Ref. [47] notes that selection rules require both atoms to become excited. Accordingly, we assume that the SIDM scatters elastically if its center of mass (c.m.) kinetic energy per particle is $<\Delta E$, and inelastically otherwise, in which case each particle loses energy equal to $\Delta E$ in the c.m. frame. The final c.m. speed of the SIDM particles after scattering inelastically is given by:

[TABLE]

The cutoff velocity $v_{c}=\sqrt{2\Delta E/m}$ plays an important role: to have any inelastic collisions, it must be less than the velocity dispersion $v_{s}$ of the halos of interest. At the other extreme, if $v_{c}$ is too low, very little energy is lost in the collisions, making the inelasticity less effective. This could lead to gravothermal collapse in dwarf galaxies or low surface brightness galaxies (LSBs) resulting in cuspy DM profiles [42] contrary to perceptions that these systems have cored profiles [26, 27, 28, 29, 30, 31, 32]. On the other hand there is evidence suggesting that not all dwarf spheroidals are cored [51, 52, 35]. In the present work we are primarily concerned with much more massive galaxies where SMBHs have been observed, so we confine our investigation to the range $200\,{\rm km}/{\rm s}\lesssim v_{c}\lesssim 500\,\rm{km}/\rm{s}$. More details are given below.

To make the simulations scale-independent, it is useful to express the $v_{c}$ in units of the maximum circular velocity of the halo, which for an NFW halo extending to $\sim 100\,R_{s}$ (as in our initial conditions in section 3) is

[TABLE]

In the second simplified model, denoted $C$, the DM interacts completely inelastically, as through forming a bound state, whose subsequent scatterings are assumed to be purely elastic, taking the same cross section for simplicity. This could mimic mirror dark matter models in which the formation of dark H2 molecules is the primary mechanism for dissipating energy [13]. A summary of the models is given in table 1.

Having established the scaling of gravothermal collapse time with the SIDM fraction $f$ in the previous section, we can reduce the noise associated with large relative fluctuations in the scattering rate by taking $f=1$, since this choice maximizes the probability for scattering. The results of three such dissipative simulations, starting from the same initial halo as in section 3, are shown in figure 3. We find that the SMBH forms within $\sim 3$-$7\,t_{r}$ for Model $B$ and $\sim 0.4\,t_{r}$ for Model $C$, in contrast to the elastic scattering result $\sim 450\,t_{r}$. The inelastic scenarios however result in smaller SMBHs, with mass approximately $0.1\%$ of the SIDM total mass for Model $B$ and $0.06\%$ for Model $C$.

The time required for collapse is thus greatly reduced relative to that found for elastically scattering DM, consistent with the results found by refs.  [13, 42]. However direct comparison with previous studies is hampered by key differences between the approaches. In ref. [13] the SIDM component was taken to be a perfect mirror sector of the Standard Model (SM) with fraction $f\cong 0.2$. Only because the mirror sector is taken to have a lower temperature than the SM, the mirror baryons can behave differently than their SM particle counterparts. The main dissipative process is formation of mirror H2 molecules by H$+e^{-}\to H^{-}+\gamma^{\prime}$ and H${}^{-}+H\to$ H${}_{2}+e^{-}$, which is sensitive to the dark photon temperature and cannot be adequately modeled by our simplified treatment.

Ref. [42] also considered the gravothermal collapse of a halo of dissipative DM, but for a single-component model with $f=1$. Constraints on the cross section are derived by demanding that gravothermal collapse does not occur in dwarf galaxies and low surface brightness galaxies (LSBs), which would create cuspy density profiles unlike those that are observed in some systems. There are two means by which SIDM could avoid having a strong impact on smaller galaxies, while still accelerating the formation of SMBHs. The first is by taking the SIDM fraction to be sufficiently small, so that even if the SIDM component undergoes gravothermal collapse it will have little impact on the combined profile. Exactly how small it should be remains a problem for further investigation. The second is by adjusting the cutoff velocity $v_{c}$ appropriately in a model with a threshold for inelasticity, like our model $B$, as mentioned in section 3.2. For observed SMBHs, we are interested in halos of mass $M_{0}\sim 10^{12}\,{\rm M}_{\odot}$ that form by $z=15$, giving a scale radius of $R_{s}\sim 1\,{\rm kpc}$, and a maximum circular velocity of $v_{\rm circ,max}=506\,{\rm km}/{\rm s}$ (see eq. (4.2)). For sufficiently large values of $v_{c}$, we can evade the bounds placed by ref. [42], as the constraints disappear for $v_{c}>200\,{\rm km}/{\rm s}$ (their parameter $v_{\rm loss}$ coincides with $v_{c}$). We therefore will confine our investigation to values such that $v_{c}\gtrsim 0.40\,v_{\rm circ,max}$.

We can also compare our predicted timescale for collapse $t_{\rm col}$ with that of ref. [42], which like us finds accelerated collapse from dissipative relative to elastic interactions, modelling dissipation similarly to our Model $B$. They determine the time reduction relative to elastic scattering (Model A), and for $v_{c}=0.12\,v_{s}$ ($v_{c}=0.25\,v_{\rm circ,max}$ for our halo) they find that $t_{\rm col}$ is reduced by a factor of 90, whereas we obtain the somewhat smaller factor of 68. For $v_{c}=0.39\,v_{s}$ ($v_{c}=0.8\,v_{\rm circ,max}$), however, the discrepancy is larger, $t_{\rm col}$ being reduced by a factor of 600 in [42] versus our value of 166. The difference may be due to the fact that we consider only dissipative scattering for $v>v_{c}$ and elastic scattering for $v<v_{c}$ rather than allowing high velocity particles to scatter both elastically and inelastically. Moreover in ref. [42] the DM scatters only if its velocity in the halo rest frame is $v>v_{c}$, whereas we impose the weaker requirement $v_{\rm rel}>2v_{c}$.

5 Comparison to observations

We now discuss simulations similar to those described in section 3 to constrain the parameters $f$ and $\sigma/m$ with respect to seeding SMBHs like those observed at high redshifts [2, 3, 1]. Because of limited computational resources, we restrict this preliminary study to a unique initial halo mass, subject to the varying scattering scenarios of our models $A$, $B$, $C$. The most favorable initial condition for explaining the observed SMBHs is a very massive halo that virializes sufficiently early. Since our focus is on comparing to the three most massive SMBHs observed, it may be reasonable to assume that these are outliers of a larger distribution, that correspond to the most massive initial halos.

The Milli-Millennium database includes a publicly available subset ($1/512$ fraction of the total volume) of the data [44] from the Millennium Simulation [45], a large-scale structure formation simulation using $\Lambda$CDM cosmology. The largest halo in the dataset at $z\sim 7$ has total mass $M_{0}\gtrsim 10^{12}\,\rm{M}_{\odot}$. Its history suggests that it virializes by $z\sim 15$. We take this to be the most favorable candidate for early SMBH formation. The distribution of halo masses from the Milli-Millenium database is shown in figure 4. The halo is atypical, having a much higher $\rho_{s}$ and smaller $R_{s}$,

[TABLE]

than halos of similar mass that form later. Eq. (4.2) gives a maximum circular velocity of

[TABLE]

We therefore simulate halo formation starting at $z=63$ (see footnote 5) from an overdensity with mass $10^{12}\,\rm{M}_{\odot}$, that will produce a halo of this mass before $z\sim 7$. Model $B$ requires a choice of $v_{c}$, that we take to be $v_{c}=0.25\,v_{\rm circ,max}$ and $v_{c}=0.6\,v_{\rm circ,max}$, using eq. (5.2). The simulations are carried out on a grid in the plane of $\sigma/m$ versus $f$, at $f=0.01$, $0.02$, $0.05$, $0.1$, $0.5$ and integer values of $\log_{10}\sigma$, for models $A$, $B$ and $C$. For each simulation, the redshift of SMBH formation is calculated, leading to contours labeled by $z$ as shown in fig. 5.

Table 2 lists the properties of the three high-$z$ SMBHs that we would like to explain by the simulations. To do so requires taking account of a degeneracy: the observed SMBH mass can be partly due to accretion after its initial formation. Taking the commonly assumed value $\epsilon_{r}=0.1$ for the radiative efficiency in eq. 1.1, this growth is described by

[TABLE]

where $t_{\rm col}$ is the time of collapse.888This accretion rate could be affected by the dissipative interactions, an effect which we do not consider here but which has been explored in [53]. In sect. 3 we saw that the black hole seed mass is a fixed fraction $\lambda$ of the total SIDM mass, depending on the model; see table 1):

[TABLE]

where $M_{0}$ is the total mass of the host halo.

Because of possible accretion, an observed SMBH can be explained by values of $f$ and $\sigma/m$ lying on curves, parametrized by the number of $e$-foldings of growth following the collapse. These are shown in fig. 5, with heavy dots marking successive $e$-foldings for the three observed SMBHs. Since the timescale for growth is faster than the Hubble rate, these curves cross the constant-$z$ contours at a shallow angle. Points where the trajectories are terminated by stars indicate the limiting cases where no accretion has occurred and the observed mass is entirely due to the initial collapse. These curves should be interpreted as lower limits on the cross section needed to explain a given SMBH, since they assume that the rate of accretion saturates the Eddington limit, and we ignore disturbances such as mergers or tidal stripping by dwarf galaxies or sub-halos that could slow SMBH formation by revirializing the halo.

It is encouraging that the trajectories for the three different SMBHs are nearly coincident, which need not have been the case. It suggests the possibility that all three systems could be explained by a single DM model, albeit with different amounts of accretion. In particular, J1342+0928 requires significantly less growth for given values of $f$ and $\sigma$ than the others because of its smaller mass. This is to be expected, since it was observed at a significantly higher redshift and thus had less time to accrete.

6 Discussion

We have demonstrated the existence of regions of SIDM parameter space that can consistently explain early SMBH formation. Figure 5 shows that, depending upon the accretion history, it is possible to match the masses and formation times of the three observed earliest-forming SMBHs for SIDM cross sections and abundances spanning several orders of magnitude. In realistic settings, one could expect larger values of $f\sigma/m$ than our idealized simulations will be required, since accretion may be less efficient than assumed in eq. (5.3). For example gas may become depleted within the vicinity of the SMBH, interrupting accretion. This may explain why not all halos with the minimal properties develop early SMBHs, making them rare events. From a particle physics perspective, very large values of $\sigma/m$ (compared for example to the Bullet Cluster limit) need not strain credulity. Atomic dark matter generically has $\sigma$ of order $\gtrsim 10\,a_{0}^{2}$, where $a_{0}$ is the mirror Bohr radius [54]. For an exact mirror of the standard model, this gives

[TABLE]

Although SIDM-induced gravothermal collapse is capable of forming very massive high-$z$ SMBHs, it will not necessarily do so in all galaxies. Our simulations assumed an isolated halo corresponding to a galaxy in the field, but most galaxies form in more chaotic environments. Mergers and the stripping of the SIDM could slow or even halt the gravothermal collapse of the halo by injecting energy and revirializing the halo, leading to the delayed formation of a SMBH. Moreover we took a special case in which the halo forms unusually early.

The SIDM mechanism of SMBH formation is not mutually exclusive with others. For example, Population III stars are able to form large black holes ($\sim 100\,\rm{M}_{\odot}$) at high redshifts, but unless they form extraordinarily early, super-Eddington accretion is required to grow them to $\sim 10^{9}\,\rm{M}_{\odot}$ by redshift $\sim 7$ [55]. The gravothermal collapse of a SIDM cloud provides a natural mechanism for super-Eddington accretion, as the radiation pressure can be far smaller or absent in the dark sector. Simulation of the accretion of an SIDM halo onto a pre-existing SMBH could be interesting for a future study,

6.1 Connection to CDM small-scale structure

Another interesting question is whether two-component SIDM is capable of addressing the small-scale structure problems of CDM that provided one of the original motivations for (single-component) SIDM [17, 18, 19, 20]. Although one may suspect that with small enough fraction $f$ there should be little effect on the central part of the DM density profile, this could depend upon $\sigma/m$ for the SIDM component, and thus far no $N$-body studies have been carried out to address this question for typical halos. It is therefore possible that the scenario we present could also have an impact on the cusp-core problem.

In fact, our inelastic models $B$ and $C$ can produce SMBHs even for $\sigma<1.0\,\rm{cm}^{2}/\rm{g}$ with $f=1$, which obeys the Bullet Cluster constraint. (Although Model $B$ ostensibly requires $f\lesssim 0.8$ to form J1342+0928 in figure 5, considering a slightly smaller initial halo would likely resolve this discrepancy.) Of course another simple way to combine the two mechanisms is to allow the principal DM component to have elastic $\sigma/m\sim 1\,$cm${}^{2}/$g, which would match the usual requirements of one-component SIDM without invalidating our results, since the dominant component would experience gravothermal collapse only on a timescale of $500\,f^{-2}\,t_{r}$, much greater than the Hubble time.

6.2 Dark disk formation

An aspect of dissipative matter that has been vigorously studied is its propensity to form a DM disk, that would overlap with the baryonic disk in Milky-Way-like galaxies [56, 57]. The existence of a dark disk in the Milky Way (MW) is strongly constrained by an analysis of recent Gaia data [58, 59], though this constraint can be evaded if the local MW halo is out of equilibrium, for example through a recent tidal disruption. So far no $N$-body simulations of dissipative DM have been done to investigate formation of a dark disk.

Ref. [56] studied dark disk formation assuming the SIDM component consisted of ionized dark atoms, leading to dissipation via bremsstrahlung interactions amongst the massive dark particles and inverse Compton scattering off a dark photon background. Here we make an order of magnitude estimate for the timescale $t_{dd}$ for dark disk formation, within our models $B$ and $C$. Defining ${\cal E}$ to be the kinetic energy density of the SIDM component and $dP/dV$ to be the kinetic energy lost per unit time and volume.

[TABLE]

We take $dP/dV=2n_{\chi^{\prime}}^{2}\sigma v\Delta E$, where $n_{\chi^{\prime}}$ is the average SIDM number density inside the virial radius, $v$ is its average velocity which we estimate as $v=\sqrt{3T_{\rm vir}/m}$, $E=(3/2)T_{\rm vir}n_{\chi^{\prime}}$, and $2\Delta E$ is the average energy lost in each collision. The MW has a virial mass of approximately $1.5\times 10^{12}M_{\odot}$ (taking the overdensity constant $\Delta_{c}$ = 200) [60] corresponding to a virial radius of 240 kpc and hence

[TABLE]

Combining these relations we determine that for the Milky Way

[TABLE]

which can be longer than the age of the universe, $13.8$ Gyr, for values of $f\Delta E/m$ that are compatible with early SMBH formation as discussed in sec. 4.1. For example with $v_{c}=200\,$km/s to evade constraints of ref. [42], $\Delta E/m\cong 2\times 10^{-7}$.

More realistic SIDM scenarios than our simplified models could have interactions between dark atoms and a dark radiation bath that might change this conclusion, but such effects are model-dependent and beyond the scope of this work. Such models must have a dark sector temperature substantially lower than that of the visible sector, to avoid dark acoustic oscillations and modifications of the matter power spectrum [38].

7 Conclusion

We have conducted the first $N$-body study of gravothermal collapse of a subdominant fraction $f$ of self-interacting dark matter, coexisiting with a dominant component of cold dark matter, as a means of seeding the early formation of supermassive black holes. This was motivated by technical limitations of an earlier hydrodynamical study, PSS14, that artificially required the normal CDM and SIDM components to maintain proportional density profiles, and which also confined its investigation to elastic scattering. Although we validate their results for the limiting case $f=1$, we find an important difference in the timescale for collapse, going as $f^{-3}$ instead of $f^{-1}$. Moreover we extended our study to include simplified models of dissipative interactions, showing that they are more effective than elastic scattering, at a fixed cross section $\sigma$.

We find that three observed SMBH’s with masses $\sim 10^{9}\,\rm{M}_{\odot}$ and redshifts $z\sim 7$ can be simultaneously explained with reasonable values of $f\sigma/m$, allowing for different amounts of accretion subsequent to collapse. Moreover, if the scattering is dissipative, a possible choice is $f=1$, $\sigma/m\cong 1$ cm2/g, which can be marginally consistent with constraints from the Bullet Cluster, while addressing puzzles about small scale structure in CDM, like the core-cusp problem.

There are several simplifying assumptions that could be improved upon in a future study. We incorporated simplified models of dissipation that are meant to capture the main features of more realistic models, where DM might form bound states, or excited states that decay by radiative emission. Our results are based upon a rare initial halo that is exceptionally large and early-forming, although still realistic in that it was taken from a large-scale cosmological simulation. We took an idealized model of subsequent accretion to describe the SMBH after initial collapse, assuming saturation of the Eddington limit, and ignoring complications such as mergers or collisions that could interrupt the SMBH growth.

More generally, the effects of dissipative interactions on structure formation is a subject that has not yet been explored in a very quantitative way, in the context of $N$-body simulations. Issues like the formation of a dark disk or distinctive effects of inelastic collisions on the small-scale structure problems represent interesting targets for future study.

Acknowledgments. We thank Guido D’Amico, JiJi Fan, Jun Koda, Roya Mohayaee, Paolo Panci, Jason Pollack, Matt Reece, Takashi Toma, Ran Huo, and Yiming Zhong for very helpful discussions or correspondence. We acknowledge Calcul Québec (www.calculquebec.ca) and Compute Canada (www.computecanada.ca) for supercomputing resources. JC thanks the Niels Bohr International Academy for its hospitality during the inception of this work, which was also supported by the McGill Space Institute and the Natural Sciences and Engineering Research Council of Canada.