Titans of the early Universe: The Prato statement on the origin of the first supermassive black holes

Abstract

In recent years, the discovery of massive quasars at $z\sim7$ has provided a striking challenge to our understanding of the origin and growth of supermassive black holes in the early Universe. Mounting observational and theoretical evidence indicates the viability of massive seeds, formed by the collapse of supermassive stars, as a progenitor model for such early, massive accreting black holes. Although considerable progress has been made in our theoretical understanding, many questions remain regarding how (and how often) such objects may form, how they live and die, and how next generation observatories may yield new insight into the origin of these primordial titans. This review focusses on our present understanding of this remarkable formation scenario, based on the discussions held at the Monash Prato Centre from November 20 to 24, 2017, during the workshop ‘Titans of the Early Universe: The Origin of the First Supermassive Black Holes’.

1. Introduction

Most massive galaxies host a central supermassive black hole (SMBH); however, the origin of these objects remains uncertain. The canonical summary of possible formation pathways was first laid out in the Halley Lecture of 1978 at Oxford University (Rees 1978). Within the so-called Rees diagram (see also Figure 1), supermassive stars (SMSs) (Hoyle & Fowler 1963; Iben 1963), dense stellar clusters (e.g., Begelman & Rees 1978), and a host of other objects were laid out as possible intermediaries. Many of these hypothesised progenitor objects were initially suggested to actually be the sources powering the emission seen in active galactic nuclei, before mounting observational evidence made it clear that these were in fact accreting SMBHs (Rees 1984). Common to all of the SMBH progenitor channels in Figure 1 is the concentration of a large quantity of gas in a sufficiently small volume, leading to runaway black hole growth. How often each channel may be realised in nature, however, remains an outstanding problem. The majority of these scenarios remain plausible but unproven today.

Figure 1. Original diagram from Rees (1978, 1984), outlining the possible formation pathways for supermassive black holes. In this review, as in the conference, our focus is on the left side of the diagram.

By far the greatest challenge to any theory of SMBH formation has been the discovery of luminous ($$ \mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} {10^{13}}{{\rm{L}}_ \odot }$$) quasars at $z\sim7$, when the Universe was only ${\sim}800\,\mathrm{Myr}$ old (e.g., Mortlock et al. 2011; Wu et al. 2015; Bañados et al. 2018). The masses of these objects are all $\mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 10^{9}\,\mathrm{M}_{\odot}$, inferred from the breadth of the observed Mg ii 2 798 Å line (e.g., McLure & Dunlop 2004), and consistent with their luminosities being near the Eddington limit. Among the most troubling examples, SDSS J010013.02+280225.8 is a known redshift ~6.3 quasar that is already $1.2 \times\,10^{10}\mathrm{M}_{\odot}$ (Wu et al. 2015), while even earlier in the Universe, ULAS J134208.10+092838.61 is a $0.8\times10^{9}\,\mathrm{M}_{\odot}$ quasar at $z=7.54$. How did these black holes reach of order 1–10 billion solar masses in the first billion years of the Universe?

The problem is best illustrated if we consider the optimistic case of persistently Eddington-limited accretion for the entire prior lifetimes of these objects. A black hole may only grow in this way from an initial ‘seed’ mass $M_{0}$ to a given mass $M_{\rm{BH}}$ in a time:

(1) \begin{equation} t_{\rm{growth}} \approx 0.45 \frac{\epsilon}{1 - \epsilon}\ln\left(\frac{M_{\rm{BH}}}{M_{0}}\right) \rm{Gyr,} \label{eqn1} \end{equation}

where $\epsilon\sim 0.1$ is the typical radiative efficiency for thin-disk accretion (see, e.g., Shakura & Sunyaev 1973, for discussion). Even in this most favourable scenario, producing a ${{\gt}}10^{9}\,\mathrm{M}_{\odot}$ quasar from a typical $$\sim10 - 100\,{{\rm{M}}_ \odot }$$ Pop III remnant would require an accretion time greater than the age of the Universe at $z\sim7$, unless significantly lower radiative efficiencies may be invoked (i.e., strongly ‘super-Eddington,’ accretion, see, e.g., Natarajan 2014; Inayoshi et al. 2016, for further discussion). Even so, numerical simulations suggest that most of such stellar-mass Pop III black holes were likely to have been ‘born starving’, and unable to grow substantially via accretion early in the Universe, particularly due to their strong ionising feedback and possible ejection from their halos via dynamical 3-body interactions (e.g., Johnson & Bromm 2007; Whalen Fryer 2012; Smith et al. 2018).

These simple considerations provide a strong hint that very rapid accretion rates and a massive ‘seed’ are necessary ingredients in the origin of the most massive high-z quasars, although relatively lower-mass progenitors may yet be plausible for these objects (see, e.g., Natarajan & Volonteri 2012; Pezzulli et al. 2017, and references therein), and can reproduce the observed M-σ relation in the local Universe (see, e.g., Taylor & Kobayashi 2014, 2015). The relative abundance of light and massive seeds, and their role in the origin of all SMBHs, depends sensitively on the prevalence of halos where the formation of massive black hole seeds is possible (e.g., Lodato & Natarajan 2007; Agarwal et al. 2012; Dijkstra, Ferrara, & Mesinger 2014; Habouzit et al. 2016a).

The first question then becomes what massive seed formation channels are viable. The collapse of dense stellar clusters to form a massive protostar has been considered for decades (e.g., Begelman & Rees 1978). The necessary inefficiency of this process, however, in which much of the original energy and angular momentum of the system are shed with stellar mass ejected in 3-body interactions, strongly limits its ability to produce extremely massive seeds (see, e.g., Latif, Schleicher, & Hartwig 2016a, and references therein). More exotic channels, such as the growth of primordial black holes (Zel’dovich & Novikov 1967), or intermediate-mass black holes (IMBH) formed from dissipative dark matter (D’Amico et al. 2018), remain somewhat speculative and require further study. A promising model for producing both large seed masses and rapid accretion rates has emerged in the atomically cooled halo scenario (e.g., Dijkstra et al. 2008), in which exposure to a strong Lyman–Werner flux from an adjacent Pop III halo destroys the molecular hydrogen in a primordial halo which has not yet undergone fragmentation (Regan et al. 2017). This prevents ‘normal’ Pop III star formation and allows infall rates of up to $$\sim0.1 - 1.0\,{{\rm{M}}_ \odot }/{\rm{yr}}$$. Numerical simulations have consistently shown such high accretion rates leading to the formation of a nuclear-burning, SMS, which will later undergo collapse through a relativistic instability (referred to in the literature as a ‘direct collapse black hole’, DCBH), leaving a massive black hole remnant (Hosokawa et al. 2013; Umeda et al. 2016; Woods et al. 2017; Haemmerlé et al. 2018a).

High baryonic streaming velocities relative to dark matter within some regions in the early Universe may provide another means to suppress star formation in low mass halos (Tseliakhovich & Hirata 2010a) and promote the growth of a DCBH seed (Tanaka & Li 2014a; Hirano et al. 2017), although such a mechanism may primarily act as a catalyst within the atomically cooled halo model (e.g., Schauer et al. 2017). Other alternative channels for producing massive seeds, including massive high-z galaxy mergers (Mayer et al. 2010; Mayer & Bonoli 2018) and the disruption of dense stellar clusters (e.g., Begelman & Rees 1978; Katz, Sijacki, & Haehnelt 2015; Sakurai et al. 2017), have also been the subject of significant investigation, although questions remain regarding the nature of the remnants they produce (see, e.g., Latif & Ferrara 2016, for a thorough discussion).

The upcoming launch of next-generation observatories such as the James Webb Space Telescope (JWST) and Euclid, as well as ground-based extremely large telescopes (ELTs), promise to yield great insight into the nature and number of the most luminous objects in the early Universe. Already, observational case studies such as the Lyman-α emitter CR7 have provided invaluable tests for our ability to distinguish between models for luminous high-z objects (Sobral et al. 2015; Pallottini et al. 2015; Hartwig et al. 2016; Agarwal et al. 2017b; Bowler et al. 2017b; Pacucci et al. 2017a). The ongoing search for IMBHs in the nearby Universe may soon allow us to constrain the distribution of black hole seed masses (Mezcua 2017). With this progress in mind, the workshop Titans of the Early Universe: The Origin of the First Supermassive Black Holes was held in late November of 2017 in Prato, Italy, with the goal of determining what theoretical and observational progress can now be made in understanding the plausibility, origin, and nature of SMSs and their direct collapse as massive black hole seeds, as well as the role of alternative formation channels. With respect to Figure 1, then, our focus is on the left side of the diagram, and in particular all paths which pass through a ‘SMS’ phase.

In the following, we begin in Section 2 by providing a summary of some of the latest theoretical developments in understanding how the extremely rapid accretion rates ($$\sim0.1 - 1\,{{\rm{M}}_ \odot }/{\rm{yr}}$$) needed to fuel the growth of massive black hole seeds and very high-z quasars may be achieved, what effects may moderate or enhance this growth, and how progress can be made in future simulations. Then, in Section 3, we outline the recent numerical results which have reached a consensus on the initial product of these extreme accretion rates: high-entropy, hydrostatic nuclear-burning SMSs, which collapse through a relativistic instability to produce massive (${\sim}10^{5}\, \mathrm{M}_{\odot}$) black hole seeds. In Section 4, we discuss the relative formation rates of light and massive seeds, the influence of still-uncertain parameterised physics, and future prospects for encapsulating the underlying physics of massive seed/DCBH formation in simplified prescriptions for use in population modelling. Finally, in Section 5, we discuss observational prospects for detecting evidence of black hole seed formation channels (and, eventually, determining the initial black hole ‘seed mass function’) using next-generation observatories, before concluding in Section 6.

2. How could SMSs form?

Independent of whether the seeds of the most massive high-z quasars are light or heavy, there is a growing consensus that sustained extremely high accretion rates are a vital ingredient in their origin. The central problem is that the growth of black holes via mergers alone is unlikely to produce ${\sim} 10^{9}\,\mathrm{M}_{\odot}$ black holes by $z\sim7$ (e.g., Sesana et al. 2007; Tanaka & Haiman 2009; Natarajan 2011, see Section 4 for further discussion). Rapid accretion rates ($$ \mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 0.1{\mkern 1mu} {{\rm{M}}_ \odot }{\mkern 1mu} {\rm{y}}{{\rm{r}}^{ - 1}}$$) are also essential if any DCBH seeds are to pass through a supermassive (${\sim} 10^{5}\,\mathrm{M}_{\odot}$) stellar evolutionary phase, as the Eddington limit imposes a maximum nuclear burning lifetime of ${\sim}2\,\mathrm{Myr}$ for these objects.

Over the last two decades, numerical simulations of the collapse of primordial halos have grown increasingly sophisticated. Early numerical work focussed on investigating the collapse of low angular momentum halos (e.g., Loeb & Rasio 1994; Begelman, Volonteri, & Rees 2006; Lodato & Natarajan 2006). It was soon understood that the dissociation of molecular hydrogen by Lyman–Werner radiation from the first stars could produce primordial halos in which collisionally excited atomic line transitions dominate cooling, allowing the isothermal collapse of such halos at ${\sim} 10^{4}\,\mathrm{K}$ (Haiman, Rees, & Loeb 1997). Since both the Jeans mass and the infall rate scale with the cube of the sound speed, and hence as $T^{3/2}$, where T is the temperature, such ‘atomically cooled halos’ provide a natural means for producing the conditions necessary for the formation of early quasars from massive seeds (e.g., Volonteri 2010, and references therein). Initial efforts to simulate the collapse of atomically cooled halos presumed the prior destruction of $\rm{H}_{2}$ (e.g., Bromm & Loeb 2003; Wise, Turk, & Abel 2008; Regan & Haehnelt 2009a, b). Subsequent work, however, has focussed on determining the initial conditions under which such a halo may arise, including the critical intensity for Lyman–Werner radiation, required to completely suppress ${\rm{H}}_{2}$ cooling within a halo (e.g., Dijkstra et al. 2008; Shang, Bryan, & Haiman 2010). This quantity (‘$J_{\rm{crit}}$’) is conventionally parameterised in the literature in terms of the specific intensity at the Lyman limit, $J_{\rm LW}$, expressed in units of $10^{-21}$ erg ${\rm{s}}^{-1}$ cm$^{-2}$ ${\rm{Hz}}^{-1}$ ${\rm{sr}}^{-1}$ (e.g., Omukai 2001). The overall picture which has emerged is that a strong Lyman–Werner flux, provided by a nearby halo which has recently undergone Pop III fragmentation, can indeed produce an atomically cooled halo (e.g., Dijkstra et al. 2008; Regan et al. 2017). The number density of DCBHs which may be formed in this way is understood to vary as $J_{\rm LW}^{-4}$ (e.g., Dijkstra et al. 2008; Inayoshi & Tanaka 2015; Chon et al. 2016). Varying estimates for $J_{\rm{ crit}}$ correspond to the ~8 orders of magnitude variation in existing estimates for the number density of DCBHs (see discussion in Section 4.2); whether a strong Lyman–Werner flux alone is primarily responsible for triggering DCBH formation is increasingly in doubt however, and remains an area of active study.

This ‘synchronised pairs’ model (Visbal, Haiman, & Bryan 2014b) implies very strict constraints on the timing, evolution, and separation of any pair of primordial halos which may provide a site for the formation of a massive black hole seed (Visbal, Haiman, & Bryan 2014a; Dijkstra et al. 2014; Regan et al. 2017). Many questions remain, however, regarding the precise conditions under which atomically cooled halos could form and plausibly reproduce the masses and space density of high-z quasars (Latif & Ferrara 2016). To what extent does fragmentation pose a problem for the formation of SMSs (e.g., Chon et al. 2016; Regan & Downes 2018a)? Can SMSs form without an atomically cooled halo, e.g., with lower-mass star formation suppressed via high baryonic streaming velocities relative to dark matter in the early Universe, or is this only complementary to the standard scenario (Schauer et al. 2017)? These questions, as well as the next steps needed to make further progress in numerical simulations of atomically cooled halos, are discussed in the following section.

2.1. ${\rm{H}}_{2}$ dissociation in primordial halos

2.1.1. Important reactions

The birth of the first stars (Pop III) in the Universe also marked the onset of stellar radiative feedback. This must then be accounted for in any self-consistent treatment of subsequent star formation. In addition to hydrogen-ionising photons (with $E{\gt}13.6\,\mathrm{eV}$), stellar populations emit photons in two other energy bands which are of vital importance in determining the thermal and chemical state of primordial star-forming halos. Following Agarwal (2018), we describe each of these in turn, before returning to address the treatment of irradiating stellar spectra in assessing the conditions necessary for the formation of atomically cooled halos (and, perhaps, massive black hole seeds).

Photons in the Lyman–Werner band, with energies $$11.2\,{\rm{eV}} \lt E \lt 13.6\,{\rm{eV}}$$, can dissociate molecular hydrogen ($\mathrm{H}_{2}$). They have insufficient energy, however, to ionise atomic hydrogen, allowing Lyman–Werner photons which escape from nascent halos to travel substantial path lengths in the early Universe, destroying ${\rm{H}}_{2}$ in other primordial minihalos and potentially delaying further Pop III star formation (Haiman, Rees, & Loeb 1996; Ciardi, Ferrara, & Abel 2000). Radiative de-excitations of excited ${\rm{H}}_{2}$ provide a powerful source of cooling in primordial halos. In its absence, the next available strong cooling term becomes radiative transitions of collisionally excited atomic hydrogen, requiring gas temperatures of ${\sim} 8\,000\,\mathrm{K}$. For the subsequent isothermal collapse, such high temperatures imply a Jeans mass of ${\sim} 10^6\,\mathrm{M}_\odot$ at $n\sim10^3\,{\rm cm}^{-3}$. Therefore, atomically cooled halos undergoing isothermal collapse permit the birth of truly supermassive objects, potentially leading to the formation of, e.g., DCBHs and/or quasi-stars (see subsequent sections).

Previous efforts to model the effect of a photo-dissociating background on star formation in the early Universe have focussed on the strength of this LW radiation ($J_{\rm LW}$, see above). Numerical studies in this vein (e.g., Shang et al. 2010; Latif et al. 2014c) suggest that it is indeed possible to almost completely dissociate ${\rm{H}}_{2}$ (and therefore suppress ‘normal’ Pop III star formation) in primordial halos, requiring however very high values of the Lyman-Werner (LW)-specific intensity (with the extent of the problem particularly dependent on the assumed spectrum, see below).

In addition to near-UV photons, stellar radiation in the infrared (IR) can also provide an important contribution to feedback in the early Universe (Wolcott-Green & Haiman 2012). In particular, photons in the energy range $h\nu \geq 0.76\,\mathrm{eV}$ can unbind the extra electron in $\mathrm{H}^{-}$, which otherwise catalyses ${\rm{H}}_{2}$ formation at low densities via the reactions:

\begin{eqnarray} \nonumber \rm H + e \rightarrow H^- + \gamma, \\ \rm H^- + H \rightarrow H_2 + e, \nonumber \end{eqnarray}

thus playing a key role in determining the equilibrium ${\rm{H}}_{2}$ fraction. Therefore, cooler low-mass stars (in particular, second generation ‘Pop II’ stars) can also play an important role in regulating early star formation. ${\rm{H}}_{2}$ formation can also be catalysed via the reaction chain

(2) $$\eqalign{ & {{\rm{H}}^{\rm{ + }}}{\rm{ + }}\,{\rm{H }} \to {\rm{ H}}_{\rm{2}}^{\rm{ + }}{\rm{ + }}\gamma , \cr & {\rm{H}}_{\rm{2}}^{\rm{ + }}{\rm{ + H}} \to {{\rm{H}}_2}{\rm{ + }}{{\rm{H}}^{\rm{ + }}}. \cr} $$

However, this is considerably less effective than forming ${\rm{H}}_{2}$ via the $\mathrm{H}^{-}$ ion, and so these reactions play an important role in SMS formation only in unusual circumstances (Sugimura et al. 2016).

Taking all this together, the photodetachment of $\mathrm{H}^{-}$ and the photo-dissociation of ${\rm{H}}_{2}$ may be parametrised as

(3) \begin{eqnarray}k_{\rm di} = C_{\rm di}\,\alpha\, J_{\rm LW}, \label{eqn3}\end{eqnarray}

(4) \begin{eqnarray} k_{\rm de} = C_{\rm de}\,\beta\, J_{\rm LW}, \label{eqn4}\end{eqnarray}

where $C_{\rm di} = 1.38\times10^{-12}\ \rm s^{-1}$ and $C_{\rm de} = 1.1\times 10^{-10}\ \rm s^{-1}$ are rate constants, α and $\beta$ are rate parameters that depend on the spectral shape of the irradiating source, and $J_{\rm{LW}}$ is the LW-specific intensity, as defined earlier. Therefore, it is essential that both the spectral shape of the irradiating source(s) and their specific intensity be accounted for in modelling the equilibrium ${\rm{H}}_{2}$ fraction in primordial halos.

2.1.2. Is there a single critical Lyman–Werner flux?

Much of the work that has been done to investigate the impact of radiation from the earliest star forming systems on the later collapse of primordial gas has adopted a highly simplified description of the spectral energy distribution (SED) of the radiation, often describing it in terms of a single temperature black body or a power-law spectrum. In particular, the use of a $T = 10^{5}\,\mathrm{K}$ black body spectrum (hereafter a T5 spectrum) to model emission from Pop III dominated systems and a $T = 10^{4}\,\mathrm{K}$ black body spectrum (hereafter a T4 spectrum) to model emission from Pop II dominated systems is particularly common, although some studies have also considered black-body spectra with intermediate temperatures (Sugimura, Omukai, & Inoue 2014; Latif et al. 2015).

The earliest study to examine the impact of extremely large levels of incident Lyman–Werner radiation on the gravitational collapse of primordial gas was carried out by Omukai (2001).^a He showed that for a sufficiently large incident flux, ${\rm{H}}_{2}$ cooling could be completely suppressed in the collapsing gas. In its absence, cooling from Lyman-α and Balmer series emission and from processes such as $\mathrm{H}^{-}$ formation results in a quasi-isothermal collapse of the gas at a temperature $T \sim 7\,000\,\mathrm{K}$, dropping only gradually with increasing density, thereby leading to a very high Jeans mass of ${\sim} 10^{6}\,{\rm M}_{\odot}$ at $n \sim 10^{3}\,{\rm cm}^{-3}$. In a low-mass dark matter minihalo, where the total gas reservoir often does not exceed this value, this would be enough to completely suppress the collapse of the gas. However, in an atomic cooling halo with a virial temperature $T_{\rm vir} \sim 10^{4}\,\mathrm{K}$ and with a much larger gas reservoir, collapse would not be suppressed. However, this high temperature would inhibit fragmentation, which leads to the idea that in such systems the gas may collapse directly without fragmenting to form a single SMS (or quasi-star) with a mass of $10^{4}$–$10^{5}\,{\rm M}_{\odot}$.

A key parameter in this scenario is the strength of the Lyman–Werner radiation field required in order to suppress ${\rm{H}}_{2}$ cooling, commonly referred to as $J_{\rm{ crit}}$. The earliest studies (Omukai 2001; Bromm & Loeb 2003) found that extremely high values of $J_{\rm {crit}}$ were required, but Shang et al. (2010) later showed that this was a consequence of the simplified treatment of ${\rm{H}}_{2}$ collisional dissociation adopted in these studies, and that simulations with a more accurate treatment of this process yielded significantly lower values for $J_{\rm {crit}}$. The issue of the correct value of $J_{\rm{crit}}$ has subsequently been addressed by a number of other authors, but despite this considerable uncertainties remain. Single-zone models of the thermal and chemical evolution of the gas, similar to those studied by Omukai (2001), typically find $J_{\rm{crit}} \sim 10$ for a T4 spectrum and $J_{\rm{crit}} \sim 1\,000$ for a T5 spectrum (see, e.g., Sugimura et al. 2014; Glover 2015a).^b However, there is a factor of 2–3 uncertainty in these values that arises from uncertainties in the input rate coefficient data (Glover 2015b). In addition, there may be as much as an order of magnitude uncertainty arising from how ${\rm{H}}_{2}$ self-shielding is treated (Sugimura et al. 2014; Glover 2015b), with different self-shielding prescriptions yielding substantially different results.

In addition, three-dimensional simulations of the collapse of highly irradiated primordial gas typically find values of $J_{\rm{crit}}$ that are much larger than those found in one-zone models, with $J_{\rm{crit}} \sim 400$ for a T4 spectrum and $J_{\rm{crit}} \sim 10^{4}$ for a T5 spectrum (Shang et al. 2010; Regan, Johansson, & Haehnelt 2014; Latif et al. 2015). The significant difference between one-zone and 3D model predictions for the value of $J_{\rm{crit}}$ is not completely understood, but may be a consequence of the large number of weak shocks that one finds in the 3D models, which produce a pervasive low level of ionisation in the collapsing gas that is missing in the one-zone models.

Additional uncertainties arise once one considers other sources of radiation. For example, if a high redshift X-ray background is present, then this can significantly increase $J_{\rm{crit}}$ if the strength of the X-ray background is high enough (Inayoshi & Omukai 2011; Inayoshi & Tanaka 2015; Glover 2016), although the effect is much greater in one-zone models than in 3D models (Latif et al. 2015). Accounting for H⁻ detachment by the Lyman-α photons produced during the collapse of the gas can also lead to a factor of a few reduction in $J_{\rm{crit}}$, potentially corresponding to up to a ~100-fold increase in the anticipated number density of DCBHs (Johnson & Dijkstra 2017).

Using a single black body spectrum to represent the integrated light from all of the stars in a galaxy is, however, a very crude approximation. A more realistic treatment necessitates incorporating more detailed models for SEDs from either single (Starburst99, Leitherer et al. 1999), or more recently available binary (BPASS, Eldridge & Stanway 2016) stellar population and spectral synthesis codes. This allows one to compute the rate parameters α and β for models of early stellar populations as a function of stellar mass, star formation rate (SFR), metallicity (Z), and age of the stellar population. Carrying out this analysis, Agarwal and Khochfar (2015) and Agarwal et al. (2017a) found that α and β, and therefore $k_{\rm de}$ and $k_{\rm di}$, vary over several orders of magnitude during the lifetime of an instantaneous burst of star formation, though less than an order of magnitude for a continuous SFR. This implies that for galaxies with realistic spectra, one cannot summarise the effects of their radiation with a single $J_{\rm{crit}}$ parameter, since the value of this will depend on α and β, and hence on the properties of the individual stellar populations. Instead, the radiation constraint for DCBH formation can be better understood in terms of the interplay between $k_{\rm de}$ and $k_{\rm di}$.

A particularly useful way to express the radiation constraint (Sugimura et al. 2014; Agarwal et al. 2016a; Wolcott-Green et al. 2017) is in the form of a critical curve in the $k_{\rm de}$–$k_{\rm di}$ parameter space (see Figure 2) which separates those sets of reaction rates that lead to efficient ${\rm{H}}_{2}$ cooling and fragmentation from those that lead to the suppression of ${\rm{H}}_{2}$ cooling and SMS formation. This curve is a consequence of the chemo-thermal evolution of the irradiated gas, and although it is affected by many of the same uncertainties that hamper our efforts to determine $J_{\rm{crit}}$ for idealised spectra, it has the significant advantage that it does not depend on the spectral properties of the source of the radiation. We can therefore study what types of sources are likely to lead to the suppression of ${\rm{H}}_{2}$ cooling by examining where the photochemical rates produced by different sources lie with respect to this critical curve.

Figure 2. The critical curve(s) in $k_{\rm de}$–$k_{\rm di}$ phase space, characterising the photo-detachment of H⁻ and the photo-dissociation of ${\rm{H}}_{2}$ in an irradiated halo (Sugimura et al. 2014; Agarwal et al. 2016a; Wolcott-Green, Haiman, & Bryan 2017). Each curve divides the parameter space into two regions, depending on the equilibrium ${\rm{H}}_{2}$ fraction resulting from the given values of $k_{\rm de}$ and $k_{\rm di}$. Rates above the curve lead to collapse into a DCBH and below result in fragmentation into Pop III stars. The difference between the curves is due to the difference in the Jeans length and self-shielding of ${\rm{H}}_{2}$. The blue curve (Agarwal et al. 2016a) is obtained with a Jeans length that is twice that of the one used in the green curve (Wolcott-Green et al. 2017).

In summary, no single-valued critical Lyman–Werner intensity $J_{\rm{crit}}$ appears to provide a suitable criterion for the formation of an atomically cooled halo under arbitrary circumstances. Alternative prescriptions based on the availability and photo-destruction of relevant species, such as that formulated by Sugimura et al. (2014), Agarwal et al. (2016a), and Wolcott-Green et al. (2017), are more useful for understanding the equilibrium ${\rm{H}}_{2}$ fraction in primordial halos irradiated by an external Lyman–Werner flux. Whether this irradiation provides the primary catalyst for DCBH formation remains an open question, with recent results adding further evidence that, e.g., structure formation dynamics may in fact play a more significant role (Wise et al. 2019). This and other complicating factors are addressed in the following sections.

2.2. Fragmentation

One of the key requirements for the formation of a supermassive primordial star (SMS) with $M_\star \mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 10^5~{\rm M}_\odot$ is to suppress vigorous fragmentation of its parent gas cloud. It is anticipated that for a nearly isothermal cloud in an atomic-cooling halo, where ${\rm{H}}_{2}$ is dissociated by intense LW radiation and/or collisions with atomic hydrogen (see the previous section), the gas undergoes monolithic collapse without a major episode of fragmentation. Significant progress has been made in the past decade to understand the details of isothermal collapse of such an ${\rm{H}}_{2}$-free gas in atomic cooling halos, both via analytical frameworks and numerical simulations. Three-dimensional hydrodynamical simulations show that the gas cloud is unlikely to fragment into gravitationally bound clumps during the runaway collapse phase and forms a central protostar surrounding a massive accretion disk (Latif et al. 2013b; Inayoshi, Omukai, & Tasker 2014; Becerra et al. 2015).

However, as material falls in, conservation of angular momentum leads to the build-up of a centrifugally supported accretion disk. As the disk grows in mass it may become gravitationally unstable, fragment, and form multiple objects (Clark et al. 2011a, 2011b; Smith et al. 2011; Greif et al. 2012; Regan et al. 2014; Sakurai et al. 2016). This is illustrated in Figure 3 for the case of an ${\rm{H}}_{2}$ dominated disk which will form normal-mass Pop III stars. High resolution simulations suggest that fragmentation mostly occurs within the central parsec, but may also happen on larger scales in some cases, depending on the larger environment, merger history, spin, and so forth. Even if fragmentation occurs within the disk on small scales, many clumps may eventually migrate inward and merge with the central protostar. This is because the migration timescale in a self-gravitating disk is as short as the orbital timescale (Clark et al. 2011a; Smith et al. 2012; Inayoshi & Haiman 2014; Latif & Schleicher 2015a). This leads to intermittent mass accretion onto the central object, with mass inflow rates up to $0.1~{\rm M}_{\odot}~{\rm yr}^{-1}$ or more (Sakurai et al. 2016; Latif, Schleicher, & Hartwig 2016b; Regan & Downes 2018a; Becerra et al. 2018). In some cases, tidal force caused by nearby galaxies strongly influences the morphology of a massive ${\rm{H}}_{2}$ cloud that potentially forms an SMS, leading to the potential formation of a smaller cluster of SMSs under a strong tidal field (Chon et al. 2016; Chon, Hosokawa, & Yoshida 2018). Clearly, this area requires further attention.

Figure 3. Fragmentation of the accretion disk in the centre of a primordial halo for conditions leading to the formation of normal-mass Pop III stars. The image is adopted from Clark et al. (2011a).

The efficiency of fragmentation depends strongly on the thermal evolution of the gas (Klessen & Glover 2016), and so an adequate description requires one to take into account all relevant processes regarding radiative cooling and chemical reactions. For this purpose, recent simulations have employed detailed chemical models in order to explore the impact of H⁻ cooling (${\rm H}+{\rm e}^- \rightarrow {\rm H}^- + \gamma$; ${\rm H}+{\rm e}^- \rightarrow {\rm H} + {\rm e}^- +\gamma$), and opacities due to H⁻ bound-free/free-free transition and Rayleigh scattering of hydrogen atoms. In fact, H⁻ cooling decreases the temperature with increasing density, i.e., $\gamma_{\rm eff} {\lt}1$, where $\gamma_{\rm eff}$ is the specific heat index ($P\propto \rho^{\gamma_{\rm eff}}$). Based on linear stability analysis, a runaway-collapsing cloud, the so-called Larson–Penston self-similar solution (Larson 1969), is unstable against non-spherical perturbations when $\gamma_{\rm eff}\mathbin{\lower.3ex\hbox{$\buildrel\lt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 1.1$. As a result, the collapsing region is elongated with the density and fragments into clumps (Lai 2000; Hanawa & Matsumoto 2000; Sugimura et al. 2017). However, since the growth timescale is sufficiently long, turbulence in the collapsing gas would suppress the unstable modes. On the other hand, when the density increases to ${\simeq} 10^6$–$10^{13}~{\rm cm}^{-3}$, cold clumps with $T\sim 10^3\,\mathrm{K}$ are produced by the thermal instability induced by the combination of the adiabatic cooling due to turbulent expansion and the subsequent H⁻ and ${\rm{H}}_{2}$ cooling (Inayoshi et al. 2014). Since the cold clumps are not massive enough to be gravitationally bound, the evolution of the central collapsing region is not affected. We note that the thermal instability would play an important role after a rotationally supported disk forms at the centre because the gravitational collapse is significantly delayed. A similar outcome is observed in less conventional direct collapse scenarios, such as in the simulations of multi-scale inflows driven by mergers between massive galaxies at $z \sim $ 8–10 (Mayer & Bonoli 2018). In the latter case, gas is already metal-enriched but the extremely high optical depth in the inner regions relegates rapid cooling and sporadic fragmentation into transient clumps just to the outer rims of the disk.

At densities above $\rm 10^{16} ~cm^{-3}$, the central region becomes opaque to all radiative processes and the temperature rises above $\rm 10^4$ K (Inayoshi et al. 2014; Van Borm et al. 2014; Latif et al. 2016b; Becerra et al. 2018). The opaque core (i.e., protostar) further increases the density and temperature adiabatically via mass accretion from the surrounding material. Due to the rapid input of high-entropy material, the central protostar evolves with an expanding envelope during the earliest stages of gas accretion (Inayoshi et al. 2014) (see also Section 3). In the late phase, gas supply onto the protostar is dominated by episodic bursts because the surrounding accretion disk is unstable against its self-gravity, forming multiple fragments. In the quiescent phases of the episodic accretion, the protostar contracts via radiative diffusion and starts to emit ionising photons. Chon et al. (2018) found that even though the ionising photons produced from the proto-SMS heats the surrounding gas, the feedback effect cannot reverse the gas inflow to the centre. This is in line with results from radiation hydrodynamic simulations of present-day star formation (see, e.g., Krumholz et al. 2009; Peters et al. 2010, 2011; Kuiper et al. 2011). However, radiative heating can reduce the level of disk fragmentation, and so further numerical simulations with a high spatial resolution and multi-frequency radiative transfer are required to fully quantify this effect. In addition, studies with a larger sample of halos in a cosmological volume will allow us to address the fate of gas collapse, fragmentation, and formation of SMSs in atomic-cooling halos.

So far, the main focus of these studies has been on understanding the isothermal collapse and its implications for the formation of an SMS. This scenario demands rather special conditions and mandates that the halo should be metal free and the formation of $\rm H_2$ remains suppressed (see previous section for details). These requirements make the sites of SMSs less abundant if not rare but this is still an open question. However, less idealised conditions such as halos polluted with trace amount of metals $Z/{\rm Z}_{\odot} \leq 10^{-5}$, halos irradiated by a moderate strength UV flux below the critical limit, or halos with large baryonic streaming velocities (see the following, dedicated subsection on this topic) may still form an SMS (Latif & Schleicher 2015b; Latif et al. 2016b; Hirano et al. 2017; Schauer et al. 2017; Inayoshi, Li, & Haiman 2018).

In fact, halos irradiated by a moderate UV flux of strength $\sim 1\,000 J_{21}$ can still provide large inflow rates of ${\sim} 0.1\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$, sufficient to produce an SMS (Hosokawa, Omukai, & Yorke 2012; Latif et al. 2014d; Latif & Schleicher 2015b; Regan & Downes 2018a). Under these conditions, gas remains hot with a temperature $\sim8\,000\,\mathrm{K}$ at scales above ${\sim}1$ parsec, similar to the isothermal case, while in the core of the halo sufficient $\rm H_2$ can form to allow the temperature to fall to a few hundred K. Similar conditions are observed for halos irradiated by an intense UV flux, but polluted with trace amounts of metals ($Z/\mathrm{Z}_{\odot} \leq 10^{-5}$). For further discussion regarding the role of streaming motions in the birth of SMS, see the following subsection.

Less idealised conditions are more prone to fragmentation on scales below 1 pc, but clump migration in combination with dynamical/viscous heating at these scales stabilises the disk. Therefore, conditions still seem conducive for the formation of a massive central object. Unfortunately, a better understanding of these conditions is still limited due to the computational constraints; however, they may provide a potential pathway for the formation of SMSs at $z {\gt} 10$.

2.3. Streaming velocities

Before recombination, baryons and photons are tightly coupled. This leads to acoustic waves (Sunyaev & Zel’dovich 1970) that also result in oscillations between baryons and dark matter. At $z \approx 1\,000$ the relative velocities are ${\sim} 30\,$km$\,\textit{s}^{-1}$, with coherence lengths of several 10–100 Mpc (Silk 1986; Tseliakhovich & Hirata 2010b). This streaming velocity decays linearly with z and reaches values of about 1 km$\,\textit{s}^{-1}$ at a redshift of $z \approx 30$. It is thus comparable to the virial velocity of the first halos to cool and collapse, and has the potential to strongly influence gas dynamics and star formation within these halos (Tseliakhovich, Barkana, & Hirata 2011; Fialkov et al. 2012).

Early studies of this process, based on numerical simulations, indeed suggest that the kinetic energy provided by the streaming velocity provides additional stability. It reduces the baryon overdensity in low-mass halos and delays the onset of cooling, altogether resulting in a larger critical mass for collapse (Greif et al. 2011; Stacy, Bromm & Loeb 2011; Maio et al. 2011; Naoz, Yoshida, & Gnedin 2012, 2013; O’Leary & McQuinn 2012; Latif, Niemeyer, & Schleicher 2014b; Schauer et al. 2017) and different angular momentum evolution of dark matter and dense gas (Chiou et al. 2018; Druschke et al. 2018). They may also have substantial impact on the resulting $21\,$cm emission (Fialkov et al. 2012; McQuinn 2012; Visbal et al. 2012). Specific to the subject of this review, it has been suggested that the presence of large streaming velocities may create the ideal conditions for the formation of SMBHs by preventing fragmentation and the formation of normal Population III stars (Tanaka & Li 2014b). Indeed, Visbal et al. (2014b) and Schauer et al. (2017) indicate that streaming velocities may suppress the formation of ${\rm{H}}_{2}$, reduce the amount of cooling, and allow halos to reach the mass limit for atomic cooling to set in. However, their results also demonstrate that once collapse has set in, ${\rm{H}}_{2}$ will rapidly form in the centre of the halo and lead to fragmentation and the onset of a more normal mode of star formation. Therefore, we may conclude that external irradiation is still needed in the presence of streaming velocities to enable the formation of SMBHs by direct collapse. Indeed, Schauer et al. (2017) argue that streaming velocities can facilitate the formation of synchronised halo pairs. In particular, if two halos form in close proximity to each other, then they can reach the atomic cooling phase at roughly the same time and start to go into collapse separated only by a few million years or less. The one forming stars earlier can then provide enough ultraviolet radiation to suppress ${\rm{H}}_{2}$ cooling in the other. If this second halo goes into collapse before the stars in the first halo explode as supernovae or before the metals from these explosions have reached it, (i.e., if strong cooling and fragmentation can be prevented), then this collapse may lead directly to the build-up of an SMBH (Agarwal et al. 2017b, 2018).

2.4. Magnetic fields

The presence of dynamically important magnetic fields has the potential to significantly alter the picture presented so far. The current Universe is highly magnetised (Beck et al. 1996), but our knowledge of magnetic fields at high redshifts is very sparse. Possible sources of magnetic fields in the early Universe are battery processes (Biermann 1950), the Weibel instability (Lazar et al. 2009; Medvedev et al. 2004), or thermal plasma fluctuations (Schlickeiser & Shukla 2003). Other theories resort to phase transitions occurring during cosmic inflation (Sigl, Olinto, & Jedamzik 1997; Grasso & Rubinstein 2001; Banerjee & Jedamzik 2003; Widrow et al. 2012). The resulting fields are typically too weak to have any dynamical impact, and so, magnetohydrodynamic effects are usually neglected when studying SMBH formation.

This simplifying assumption could be incorrect, however, due to the small-scale turbulent dynamo that is active during all phases of structure formation. This is able to efficiently amplify even extremely small primordial seed fields to the saturation level. The amplification timescale is comparable to the eddy-turnover time on viscous or resistive length scales (Kazantsev 1968; Subramanian 1998; Sur et al. 2010, 2012; Schober et al. 2012a, 2012b) depending on the Prandtl numbers, and in either case is much shorter than the dynamical or collapse timescale. The resulting magnetic field configuration is extremely complex, as illustrated in Figure 4, with the energy density being dominated by fluctuations on the smallest scales accessible to the system. Once backreactions become important, the growth rate slows down, and saturation is reached within a few large-scale eddy-turnover times (Schekochihin et al. 2002, 2004; Schober et al. 2015), which again is shorter than the typical lifetimes of protostellar accretion disks. The field in this phase is dominated by the more coherent larger-scale rotational structure of the accretion disk. The magnetic energy density at saturation is thought to reach levels of 0.1% up to a few × 10% of the kinetic energy density (Federrath et al. 2011a, 2014).

Figure 4. Magnetic field lines (yellow: 0.5 ${\rm \mu}$G, red: 1 ${\rm \mu}$G) illustrating the complex magnetic field structure in the centre of a contracting primordial halo resulting from the small-scale dynamo. The image is adopted from Federrath et al. (2011b).

Dynamically significant magnetic fields can strongly affect the evolution of protostellar accretion disks. They remove angular momentum from the star-forming gas (Machida, Matsumoto, & Inutsuka 2008; Machida & Doi 2013; Bovino et al. 2013; Latif et al. 2013a, 2014a), drive protostellar jets and outflows (Machida et al. 2006), and most importantly for our discussion, they reduce the level of fragmentation in the disk (Turk et al. 2011; Peters et al. 2014) and thus may facilitate the formation of SMBHs by direct collapse. To better understand this process, and in particular, to quantitatively determine how magnetic fields influence the dynamical evolution of the gas, remains an area of very active research.

2.5. Black hole seed growth: simulation requirements

Chon et al. (2018) and Regan and Downes (2018b) recently modelled the formation of SMSs in cosmological environments with semi-analytical recipes for the evolution of the SMS using the Gadget and Enzo codes, respectively (Springel 2005; Bryan et al. 2014). Both simulations included radiative transfer to approximate stellar feedback from the SMS and thus make progress towards modelling the correct thermal state of the gas surrounding the SMS. The SMS is thought in most cases to collapse to a black hole at the end of its life because it either encounters the general relativistic (GR) instability above a few hundred thousand $\rm{M_{\odot}}~$ (Woods et al. 2017) or it exhausts its nuclear fuel.

The DCBH is likely born in the high gas densities in which the SMS had grown, with infall rates of $10^{-2}$–$10^{-1}$ $\rm{M_{\odot}}~$ yr${-1}$. Accretion onto the black hole is regulated by a number of processes on scales down to nearly its event horizon, which cannot be resolved in current cosmological simulations. Therefore, the growth of the black hole is usually approximated by Bondi–Hoyle (B–H) accretion. B–H accretion assumes that flows onto the black hole are spherical, and it is characterised by the B–H radius, $R_{\mathrm{B-H}}$:

(5)\begin{equation} R_{\mathrm{B-H}} = 2\,G\,M_{\mathrm{BH}}/c^2_\mathrm{s}, \label{eqn5} \end{equation}

where $M_{\mathrm{BH}}$ is the black hole mass and $c_\mathrm{s}$ is the sound speed in the surrounding medium. Although flows onto the black hole are not spherical, the B–H model provides a reasonable benchmark for accretion rates. In numerical simulations, the accretion rate is usually just taken to be the inflow rate through $R_{\mathrm{B-H}}$ and the luminosity of the black hole is computed as $\epsilon \dot{m}c^2$, where $\epsilon$ is the radiative efficiency, the fraction of the mass that is converted to energy upon being swallowed up by the black hole. Past simulations have simply deposited this luminosity locally as heat (e.g., Di Matteo et al. 2012; Feng et al. 2014; Costa et al. 2014; Hirschmann et al. 2014) but the newest models now perform multi-frequency radiation transport from the black hole (Aykutalp et al. 2014; Smidt, Wiggins, & Johnson 2016; Smidt et al. 2018) as well as mechanical feedback in the form of jets (Regan et al. 2018c).

If the medium is not heated by the SMS, a 50 000 $\rm{M_{\odot}}~$ DCBH will have a B–H radius of ${\sim}10^{-2}$ pc at birth. However, this radius can drop by up to three orders of magnitude when the black hole begins to accrete because X-rays and outflows heat the gas and raise its sound speed by a similar factor. The B–H radius of an SMS can be resolved in cosmological simulations, but doing so for a DCBH would restrict the simulation to such small time steps that it would be difficult to evolve the black hole or study its effect on its own subsequent growth. Consequently, resolving the B–H sphere of an SMBH seed will remain a challenge for the foreseeable future.

2.5.1. Modelling accretion and feedback from a seed black hole below the Eddington limit

As already noted above, a rapidly accreting black hole seed will generate significant feedback, which must be accurately modelled by the simulation in order to correctly estimate the resulting growth. Feedback from a black hole can be broken down into two different components—radiative and mechanical. Furthermore, the regimes of black hole accretion are thought to lead to different feedback characteristics. Below an accretion rate of $\dot{M} \sim 10^{-3} \dot{M}_\mathrm{Edd}$ radiative feedback is expected to be minimal as the material surrounding the disk is optically thin and radiative feedback ineffective. In this case, mechanical feedback in the form of radio jets is thought to be the dominant feedback mechanism. Above this limit but below the Eddington limit (i.e., $10^{-3} \dot{M}_\mathrm{Edd} \mathbin{\lower.3ex\hbox{$\buildrel\lt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} \dot{M} \mathbin{\lower.3ex\hbox{$\buildrel\lt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} \dot{M}_\mathrm{Edd}$), radiative feedback is expected to be the dominant feedback mechanism. Then for super-Eddington flows, the radiative feedback component may become trapped by the rapid inflow; radiative feedback may then once again become inefficient, though jets may again dominate (see Figure 5).

Figure 5. Figure taken from Sądowski and Narayan (2016) showing the regimes where both radiative and mechanical feedback are expected to dominate. $a_*$ is the normalised spin of the black hole. Radiative efficiency is plotted on the y-axis against accretion rate on the x-axis.

However, it should be noted that these regimes are contentious and jets may in fact be active at all phases of the accretion scale (e.g., Sądowski & Gaspari 2017).

Radiative feedback from the accretion disk surrounding a black hole is the dominant form of feedback in the regime between radio-mode accretion and super-Eddington accretion. The accretion onto a black hole and the radiative feedback resulting from this accretion is usually captured by the radiative efficiency parameter, $\epsilon$. For a geometrically thin but optically thick (Shakura–Sunyaev, Shakura Sunyaev 1973) disk accreting close to, but below, the Eddington limit, the radiative efficiency is $\epsilon = 0.057$ for a non-spinning black hole. Radiative feedback from the accretion disk can be modelled in a straightforward manner, by determining the appropriate SED for a given black hole mass and accretion rate (e.g., Done et al. 2012).

Thermal feedback is often used as an approximation to model the impact of both radiative and mechanical feedback. This is simpler and less computationally expensive to invoke, but may miss the essential physics of radiative and/or mechanical feedback. However, in the regime where a black hole is accreting above the Eddington limit, the SED is unknown and may not necessarily follow that of the relatively well-understood thin disk. In this case modelling the radiative feedback is problematic as we do not know the SED and the only recourse, at the present time, may be thermal feedback modelling.

2.5.2. Super-Eddington accretion: slim disk accretion

The slim disk model of accretion onto a black hole surface is derived from numerical integration of the Navier–Stokes equations, where the point at which the gas velocity exceeds the local sound speed is set as a critical point (Abramowicz et al. 1988). The slim disk solutions can be viewed as a generalisation of the thin disk solutions derived by Shakura and Sunyaev (1973). Although different models for super-Eddington accretion exist (e.g., Jiang, Stone, & Davis 2014) we focus here on a discussion of the slim disk model as it is widely applied and captures the essential physics of super critical accretion flows.

The slim disk model is valid for accretion flows which exceed the Eddington limit, a regime in which the thin disk solutions break down (in fact the thin disk model breaks down at approximately $\dot{m} \sim 0.5 \dot{m}_\mathrm{Edd}$). Black holes accreting at rates above the Eddington limit have been shown to generate powerful jets (McKinney, Tchekhovskoy, & Blandford 2013; Sądowski et al. 2014, 2016; Jiang et al. 2014) which may significantly disrupt the accretion flow. While the jets have been shown to be highly collimated at the scales of the high resolution general relativistic radiative magneto-hydrodynamic (GRRMHD) simulations used to model them, the scales at which these simulations operate ($\mathbin{\lower.3ex\hbox{$\buildrel\lt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}$1 000 Schwarzschild radii) are typically orders of magnitude separated from the scales at which astrophysical codes can operate. Furthermore, as the bow shock from the jet propagates outwards, it will broaden and disrupt gas perpendicular to the disk. Simulations of jets from accreting seed black holes which use the results of much higher resolution, much smaller scale GRRMHD simulations are only now being attempted.

Furthermore, as the accretion flow onto a black hole exceeds the Eddington rate, the radiative efficiency of the feedback has been shown to reduce—in some cases, significantly. Thus in the super-Eddington regime, the dominant form of feedback can be mechanical. For black holes accreting above the Eddington limit, Madau, Haardt, & Dotti 2014 provide an empirical fit to the GRRMHD simulations of Sądowski (2009), where the radiative efficiency of a black hole accreting at rates above the critical Eddington limit has indeed been shown to decrease, thus potentially allowing more accreted material to reach the black hole surface compared to the Eddington limited case. Of course, the net effect of the jet needs to be accounted for to determine whether the jets positively or negatively impact the accretion onto the black hole.

In accretion disk simulations, the degree to which the radiative efficiency can decrease is, however, not yet established, as it appears to depend on the specific way by which radiation transport is treated, as well as other aspects of the numerical implementation (see review chapter by Mayer and Bonoli (2018), for a detailed discussion). In particular, the recent radiative magneto-hydrodynamic simulations of Jiang, Stone, and Davis (2017), which are among the few to study accretion on massive black holes rather than stellar mass black holes, exhibit radiative efficiencies that are an order of magnitude higher than those in the slim disk model, albeit being still almost an order of magnitude lower than standard viscous disk efficiencies. Interestingly, radiative efficiencies are also found to decrease with increasing accretion rates, being as small as $1\,\%$ of $\dot M c^2$ for accretion rates above $100\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$.

Finally, assuming that rapid accretion onto a seed black hole can be achieved, how easily can the black hole accrete the inflowing material? Sugimura et al. (2018) investigated the angular momentum barrier for black holes attempting to accrete at high-z. Using 2D axially symmetric radiation hydrodynamic simulations, they found that the interplay of radiation feedback and angular momentum transport significantly suppresses accretion. They concluded that in order for the black hole to efficiently accrete the incoming material, the angular momentum of the gas must be significantly suppressed within a centrifugal radius $\mathbin{\lower.3ex\hbox{$\buildrel\lt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 0.01 \times R_{\rm{Bondi}}$. At face value this appears to introduce yet another significant barrier to achieving rapid black hole growth onto seeds at high-z. However, the much smaller scale GRRMHD simulations which model the magneto-rotational instability in a self-consistent manner can efficiently transport angular momentum outwards (e.g., Sądowski & Gaspari 2017) showing that rapid accretion is possible, at least at small scales. Lupi et al. (2016) carry out 3D hydrodynamical simulations of massive circum-nuclear disks, a few pc in size, which are likely present in gas-rich galaxies at high redshift. Using both adaptive mesh refinement (AMR) grid-based codes and Lagrangian mesh-less particle codes, they show that such systems become quickly clumpy and gravitoturbulent. The large scale turbulence arising from gravitational instability drives efficient angular momentum transport through the disk down to sub-pc scales. By initialising stellar mass black hole seeds and assuming a slim disk model as a sub-grid recipe to treat black hole accretion, these authors suggest that there is indeed a bottleneck in the accretion flow from large to small scales, even when radiative heating of the ambient medium originating from accretion is taken into account. These simulations, however, still fail to resolve the gas flow at scales significantly below a few parsecs.

Further high resolution simulations, incorporating both the physics of SMS formation, evolution, and transition into a black hole, and the feedback effects of the resulting accretion onto the black hole, will be required to determine how attractive this pathway to eventual SMBH formation may be. As with the case of growing Pop III remnants (the so-called light seeds), accretion can be shut-off and it is as yet unclear what conditions are required to achieve optimal growth possibly including super-Eddington accretion.

3. How do SMSs live and die?

SMSs were originally conceived as a model for quasars themselves, rather than their progenitors (Hoyle & Fowler 1963 Iben 1963; Fowler 1964). Early analytic developments for the most part ignored how an SMS may have originated, assuming a single monolithically contracting cloud, and neglecting the impact of its formation history on its subsequent evolution. In this case, many of the physical characteristics of SMSs can be computed analytically, being well approximated as polytropes of index 3. Such stars are, however, extremely prone to instability (Chandrasekhar 1939). It was soon understood that even in the relatively modest gravitational potential of SMSs, post-Newtonian corrections to the equation of state would place a firm upper limit on its mass of ${\sim} 10^{5}\,\mathrm{M}_{\odot}$, above which the star would collapse on its dynamical timescale (Chandrasekhar 1964). This led to a sustained interest in their evolution as potential progenitors of SMBHs (Rees 1984). Later, numerical experiments using the detailed stellar evolution code Kepler (Weaver, Zimmerman, & Woosley 1978) suggested that for primordial composition SMSs, any collapse inevitably led to the formation of a black hole, preceded at lower masses by stable hydrogen burning (Fuller, Woosley, & Weaver 1986, see also Figure 6). More recently, it has been suggested that there could be a narrow range of masses around $50\,000\,\mathrm{M}_\odot$ for which primordial, ‘monolithic’ SMSs may instead undergo a thermonuclear explosion, leaving no gravitationally bound remnant (Chen et al. 2014). It remains uncertain, however, whether this result is robust against, e.g., details in assumptions and modelling of the underlying physics. The narrowness of the mass window may make such explosions very rare and hence hard to exclude by direct observations, but see Section 5 for further discussion.

Figure 6. Evolutionary outcomes for monolithically formed supermassive stars as a function of mass andmetallicity Z. Adapted from Fuller et al. (1986), although see Chen et al. (2014). Note that no rapidly accreting models have been found to explode (e.g., Umeda et al. 2016; Woods et al. 2017).

Producing a single, supermassive, contracting cloud as a precursor to an SMS, e.g., through the collapse and coalescence of a dense stellar cluster, appears unlikely, however, to produce stellar masses greater than approximately a few $10^{3}\,\mathrm{M}_{\odot}$ (see Latif & Ferrara 2016, and references therein), though other channels may be viable (see, e.g., Mayer et al. 2010; Mayer & Bonoli 2018). In the atomically cooled halo scenario, a supermassive proto-star instead grows (relatively) gradually at accretion rates of order $0.1$–$10\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$. This has a profound effect on the structure of the star: it never reaches thermal equilibrium (Begelman 2010), and therefore detailed stellar evolutionary calculations are needed to predict the star’s observational characteristics and fate (Hosokawa et al. 2012). In particular, although such stars still encounter the post-Newtonian instability and collapse, simple polytropic models underestimate the final mass (see, e.g., discussion in Woods et al. 2017). For non-rotating SMSs, the general properties of accreting SMSs are now well understood (Hosokawa et al. 2013; Haemmerlé et al. 2018a). In particular, the effective temperatures of rapidly accreting SMSs are found to be significantly cooler than in the monolithic case, and no thermonuclear explosions are expected to reverse their collapse to black holes. There is also a broad agreement on the final mass at collapse as a function of infall rate. The largest discrepancies in final masses are found at the highest accretion rates, up to about a factor of 2 (Umeda et al. 2016; Woods et al. 2017; Haemmerlé et al. 2018a).

Beyond this qualitative agreement, however, significant differences remain in the predicted effective temperatures and final evolutionary states of SMSs found by the different numerical studies (Section 3.2). Furthermore, efforts to include rotation in detailed evolutionary calculations have clearly demonstrated a bottleneck in the growth of any SMS; an extremely efficient mechanism for angular momentum removal must be in place if an accreting proto-star is to avoid spinning-up to the mass-shedding limit before reaching the SMS regime (Haemmerlé et al. 2018b). How this may take place remains uncertain, although this problem is known to exist for massive star formation at all metallicities in general. In the following, we discuss the latest progress on these issues, as well as benchmarking stellar evolutionary calculations from the participating groups.

3.1. The structure and evolution of monolithic SMSs

What do we mean when we say a star is supermassive? Attempts at naming schemes for ‘big’ stars based strictly on total mass have been made in the past (e.g., Heger & Woosley 2002; Heger et al. 2003; Heger & Woosley 2010; Heger 2012); here, however, we use an approach based on the relevant physics. In order to arrive at a clear definition, we must first consider the structure of massive stars in general. For a gravitationally bound sphere of ideal gas, hydrostatic equilibrium requires greater central temperatures at greater total masses. This means that for higher and higher stellar masses, more and more of the total pressure in the star comes from radiation. For very high masses, gas pressure is only a small perturbation to the total pressure. In this case, it can be shown that, to a good approximation, the total pressure scales as $\propto \rho ^{4/3}$; more formally, the star is well approximated as a polytrope of index 3 (Kippenhahn, Weigert, & Weiss 2012). The temperature gradient in a radiation-dominated $n\approx3$ polytrope is almost equal to the adiabatic gradient, and therefore nuclear burning easily drives convection throughout the star. Consequently, convective mixing enforces a flat entropy profile; in the high-entropy limit relevant to massive, radiation-dominated stars, the adiabatic pressure gradient may then be well approximated as (Chandrasekhar 1939):

(6)\begin{equation} \Gamma = \left(\frac{\partial\ln P}{\partial\ln \rho}\right)_\mathrm{\!\!ad} = 1 + \frac{1}{n} \approx \frac{4}{3} + \frac{\beta}{6}, \label{eqn6}\end{equation}

where the ratio of gas pressure to total pressure, β, can be approximated in the limit $\beta \ll 1$ as

(7)\begin{equation} \beta = \frac{P_\mathrm{g}}{P_{\rm{tot}}} \approx \frac{4.3}{\mu}\left(\frac{\mathrm{M}_{\odot}}{M}\right)^{1/2}, \label{eqn7} \end{equation}

where ${\rm \mu}$ is the mean molecular weight of the gas, assuming the star is chemically homogeneous. For very large M, the pressure gradient $\Gamma \rightarrow 4/3$. This ‘softness’ of the equation of state leaves the star very unstable against small perturbations (Fowler 1964; Goldreich & Weber 1980; Kippenhahn et al. 2012). For this reason, a number of typically negligible effects must be taken into account in any detailed numerical modelling of extremely massive stars (Fuller et al. 1986), including, perhaps surprisingly, general relativistic corrections to the equations of stellar structure.

The critical role post-Newtonian effects play in the fate of extremely massive, quasi-statically contracting stars and proto-stars was first noted in the works of Hoyle and Fowler (1963) and Iben (1963). To summarise (following Fuller et al. 1986), hydrostatic equilibrium requires that the total energy of a star, a function of its total mass, central density, and entropy, be at an extremum. An initially low-density/high-entropy supermassive proto-star, produced by, e.g., the disruption and collapse of a dense stellar cluster (Begelman & Rees 1978) will lose energy to radiation and contract, increasing its central energy and decreasing its equilibrium energy. The effect of including general relativistic corrections to the equations of stellar structure, however, is that the equilibrium energy will have a minimum. As the star contracts, its central density may exceed a critical threshold $\rho _{\rm{crit}}$, above which energy must be added to achieve hydrostatic equilibrium. Therefore, if a quasi-statically contracting polytrope reaches this point before nuclear burning can stabilise its evolution, it will undergo a dynamical collapse. This instability was first formally derived by Chandrasekhar (1964) and Feynman (1963), and in the literature is variously called the Chandrasekhar, Chandrasekhar-Feynman, or post-Newtonian instability. For a polytrope of index 3, the critical central density for the onset of a dynamical instability is

(8)\begin{equation} \rho _{\rm{crit}} \approx 1.994 \times 10^{18} \left(\frac{0.5}{\mu}\right)^{3}\left(\frac{\mathrm{M}_{\odot}}{M}\right)^{7/2} \mathrm{g}\,\mathrm{cm}^{-3}. \label{eqn8} \end{equation}

Another, perhaps more intuitive, way to express this criterion is as a condition on the critical pressure gradient needed to support the star against collapse; for a fully Newtonian star, this is $\frac{4}{3}$, whereas in the post-Newtonian limit, general relativity requires a slightly steeper gradient:

(9)\begin{equation} \Gamma_{\rm{crit}} \approx 4/3 + 1.12\,\frac{2\,G\,M}{R\,c^{2}}, \label{eqn9} \end{equation}

where the coefficient $1.12$ is valid for the specific case of a polytrope with index $n=3$ (cf., Chandrasekhar 1964). Comparing Eq. (9) with Eq. (6), we see that the adiabatic pressure gradient in a fully convective, radiation-dominated polytrope of $n=3$ falls below the threshold needed to support it against collapse above a characteristic mass:

(10)\begin{equation} M_{\rm{SMS}} \approx \left( 0.32\; \frac{R\,c^{2}}{\mu\,G\,\mathrm{M}_{\odot}}\right)^{2/3}\,\mathrm{M}_{\odot} \sim 10^{5}\,\mathrm{M}_{\odot} \label{eqn10} \end{equation}

This defines the realm of truly SMSs; very massive stars below this regime (roughly $10^{2}$–$10^{4}\,\mathrm{M}_{\odot}$) invariably survive to nuclear-burning and collapse on the electron-positron pair instability (e.g., Bond, Arnett, & Carr 1984). Understanding the boundary between these two regimes, however, and the role nuclear burning can play in supporting the star against collapse, or even reversing this and producing an explosion, requires detailed stellar evolution calculations.

Computational modelling efforts to follow the evolution and fate of SMSs in detail were carried out as early as the 1970s (e.g., Appenzeller & Fricke 1972a, 1972b). The first systematic effort in modelling these objects to include a comprehensive treatment of nuclear burning (including the rp-process, Wallace & Woosley 1981), post-Newtonian corrections, and hydrodynamics was that of Fuller et al. (1986). Beginning with a grid of highly inflated (high-entropy) polytropes of index $n=3$, Fuller et al. (1986) allowed their models to quasi-statically contract and followed the onset of nuclear burning. They found that SMSs of primordial composition either collapsed directly to black holes within approximately the Kelvin–Helmholtz time, or survived to undergo the hydrogen-burning main sequence. Modern simulations incorporating a complete nuclear reaction network place the boundary between these two regimes at ${\approx}155\,000\,\mathrm{M}_{\odot}$ (Woods et al. 2019). It was found that more massive SMSs, however, with an initial metallicity, $Z\mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}1/4\,\rm{Z}_{\odot}$ could be completely disrupted in thermonuclear explosions, not seen in the primordial case (though see Chen et al. 2014).

This behaviour is easily understood. For primordial SMSs, the absence of Carbon-Nitrogen-Oxygen (CNO) elements necessitates that the star contract until the onset of He-burning through the triple-alpha reaction, as hydrogen-burning through the pp-chain alone cannot stabilise massive stars against further contraction (Heger & Woosley 2010; Kippenhahn et al. 2012). This means that vigorous nuclear burning begins in primordial SMSs only after most of their mass has fallen much deeper into their potential wells. For SMSs above some threshold, nuclear burning is unable to reverse the collapse. For SMSs of higher metallicity, with abundant CNO to catalyse nuclear burning, this is not the case, and thermonuclear explosions are possible in some cases (see Figure 6).

It remains uncertain, however, whether conditions under which non-zero metallicity SMSs could form ever arise in the Universe. In particular, the disruption of dense stellar clusters alone does not appear to be able to produce sufficiently massive objects (see Latif & Ferrara 2016, for a review). The merger of massive, gas-rich proto-galaxies may be capable of producing SMSs of arbitrary metallicity; however, it is not clear at present whether the structure of such stars would in any way resemble the monolithic polytropes explored above. It is clear, however, that the rapidly accreting supermassive protostars produced in the atomically cooled halo scenario, feasible only for primordial or extremely low metallicities, emerge with a radically different internal structure (e.g., Hosokawa et al. 2012, 2013; Umeda et al. 2016; Woods et al. 2017; Haemmerlé et al. 2018a). This has profound consequences for their evolution and fate, which we outline in detail in the following subsection.

3.2. Rapidly accreting SMSs

3.2.1. Non-rotating case

Rather than forming ‘all-at-once’ from the contraction of a single monolithic cloud, the prevailing theory for the formation of massive black hole seeds is now the ‘synchronised pairs’ scenario, giving rise to an atomically cooled halo (see Section 2). In this picture, SMSs form under extreme accretion rates of ${\sim} 0.01$–$10.0\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$ onto a small central seed. This means that under typical conditions, the time needed for the star to reach the supermassive regime (${\mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}} 10^{5}\,\mathrm{M}_{\odot}$) is comparable to its nuclear-burning lifetime.

This has profound consequences for the evolution of the rapidly accreting protostar (Begelman 2010). In particular, early in its lifetime the accretion timescale

(11)\begin{equation} t_{\rm{acc}} = \frac{M_{*}}{\dot M} \label{eqn11} \end{equation}

is considerably shorter than the Kelvin–Helmholtz (thermal) timescale

(12)\begin{equation} t_{\rm{KH}} = \frac{GM^{2}_{*}}{R_{*}L_{\rm{int}}} \label{eqn12} \end{equation}

with $L_{\rm{int}}$ the internally generated luminosity of the star and $R_{*}$ its radius. Unable to cool to thermal equilibrium, these stars dramatically expand before reaching the main sequence (Omukai & Palla 2001, 2003). The internal structure of such stars was first investigated by Begelman (2010), who modified the original polytropic formulation in order to incorporate a mass-dependence in the equation of state (the so-called ‘hylotropes’). Since the continuously growing outer layers of a rapidly accreting SMS cannot generally relax within its nuclear-burning lifetime, Begelman showed that the structure of the star on the main-sequence would deviate dramatically from that of an $n=3$ polytrope, with an inner, nuclear-burning convective core surrounded by a high-entropy envelope (see also Figure 7). In the time it takes to become supermassive, the core will have already evolved significantly. Accretion continues to proceed at a rate faster than the envelope can cool, so that the envelope’s entropy profile steeply increases outwards.

Figure 7. Kippenhahn diagram for a supermassive star accreting $1\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$, showing the structure of the star at any given mass coordinate as a function of time. Green single-hashed regions are convective, and blue double-hashed regions are radiative. The blue shading denotes energy generation, with the scale given on the right axis. From Woods et al. (2017).

The first efforts to produce numerical simulations of rapidly accreting SMSs followed shortly thereafter. Hosokawa et al. (2012) followed the early evolution of such an object until reaching ${\sim}1\,000\,\mathrm{M}_{\odot}$. They found that the effective temperature of the inflated protostar fell to an almost constant ${\approx}5\,000\,\mathrm{K}$, approximately the Hayashi limit for primordial stars (Hayashi 1961; Volonteri 2010). The extreme temperature sensitivity of H^– absorption, the primary source of opacity in this regime, prevents the star from cooling and its outermost layers from expanding further. Given that the luminosity of a massive, radiation-dominated SMS throughout its evolution will be very near the Eddington limit ($\propto M_{*}$), the Stefan–Boltzmann law provides a relation between the radius of an SMS and its mass:

(13)\begin{equation} R_{*} \approx 2.3\times 10^{3}\,\mathrm{R}_{\odot} \left(\frac{M}{100\,\mathrm{M}_{\odot}}\right)^{1/2} \label{eqn13} \end{equation}

closely approximating numerical results (Hosokawa et al. 2012).

During the star’s subsequent evolution, (Hosokawa et al. (2013) found that the radius continues to grow monotonically with the mass until the star reaches ${\approx}10^{4}\,\mathrm{M}_{\odot}$. At this point, the ever-lower gas density in the outer envelope falls too low for H⁻ absorption to be an effective source of opacity, and the SMS contracts somewhat (see also Haemmerlé et al. 2018a). Critically, however, the SMS never grows hot enough to become a significant source of ionising photons, in sharp contrast to the monolithic case (Johnson et al. 2012), reaching only ${\sim} 1-2 \times 10^{4}\,\mathrm{K}$ in detailed numerical simulations using different implementations (e.g., Hosokawa et al. 2013; Umeda et al. 2016; Woods et al. 2017; Haemmerlé et al. 2018a). Therefore, in contrast to the case for lower accretion rates (${\sim} 10^{-3}\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$) typical of massive star formation in the present-day Universe (Omukai & Palla 2001), the growing luminosity of accreting SMSs never becomes self-limiting, since the star always remains too cool to ionise a surrounding H II region (Johnson et al. 2012; Hosokawa et al. 2013). This is essential to the viability of SMSs as the progenitors of the first black holes in the atomically cooled halo scenario.

A potential caveat to this, however, is that all of the studies discussed above assumed constant accretion rates throughout the lifetimes of their model SMSs. Under more realistic conditions, however, accretion may be highly variable, e.g., if mediated by a disk which becomes gravitationally unstable and fragments (e.g., Stacy, Greif, & Bromm 2010). Indeed, such disk fragmentation is expected in atomically cooled halos (see, e.g., Regan et al. 2014). The resulting ‘clumps’ of matter would then migrate inward through the disk and merge with the central SMSs in bursts of rapid accretion. Sakurai et al. (2015) were the first to examine in detail the effect this process would have on the evolution of SMSs. They found that in the quiescent phases in between these bursts, the outer envelope of the SMS would begin to thermally relax and contract. Consequently, the effective temperature of the SMS may grow high enough to produce a substantial flux of ionising photons. If left unchecked, this phenomenon could provide a significant source of feedback, potentially halting the growth of the star. Subsequent bursts of accretion, however, will invariably re-inflate the SMS. Sakurai et al. (2015) found that for quiescence durations greater than

(14)\begin{equation} \Delta t_{\rm{quiescent}} \mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} t_{\rm{KH,env}} \sim \text{few}\times 10^{3}\,\mathrm{yr} \label{eqn14} \end{equation}

the ionising photon luminosity of the contracting SMS may grow sufficient to produce an H II region (actually, this is the point in their calculation where one ionising photon is produced by the SMS for every atom accreted, a simplifying assumption). Numerical simulations of the collapse of atomically cooled halos suggest quiescence times in excess of a few 1 000 yr may be expected (e.g., Inayoshi & Haiman 2014). UV feedback may therefore play a role in limiting the growth of SMSs, although more detailed radiative transfer calculations are needed.

Perhaps the most vital question regarding the evolution of SMSs, at least in the context of the formation of high-redshift quasars, is the final mass and fate of such objects. Here again the radically different structure of rapidly accreting SMSs sharply contrasts with the monolithic case. In particular, the lower densities of inflated, rapidly accreting SMSs should allow them to avoid collapse at higher masses than otherwise expected (Begelman 2010; Umeda et al. 2016; Woods et al. 2017). Classically, a monolithic SMS undergoes the post-Newtonian instability once its core density satisfies Eq. (8); the pressure support in an $n=3$ polytrope can no longer hold the star up against collapse. Another way of expressing this constraint is evident from comparison of Eqs. (6) and (9); both the critical pressure support needed to prevent collapse (including the post-Newtonian correction) and the adiabatic pressure gradient are very nearly 4/3, each including only a small perturbation. Equating the two gives us yet another formulation of the condition for instability:

(15)\begin{equation} \frac{\beta}{6} {\lt} 1.12\,\frac{2\,GM}{R\,c^{2}}, \label{eqn15} \end{equation}

which, if satisfied within a fully convective, radiation-dominated star, should signal the onset of collapse.

An SMS, however, is emphatically not well-approximated by an $n=3$ polytrope, with the specific entropy rising steeply throughout its outer envelope. An illustrative example is given in Figure 8, for the extreme case of an SMS accreting $10\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$. Here we have plotted either side of the inequality above; for a fully convective star, the onset of collapse should arise when the post-Newtonian correction to the critical pressure support exceeds that provided by gas pressure. For rapidly accreting stars, however, the structure of the star at this point is not well modelled as an $n=3$ polytrope, even when only the inner convective core is considered. Incorporating implicit hydrodynamics (and post-Newtonian corrections to the equations of stellar structure) allows one to follow the onset of collapse in detail using a detailed stellar evolutionary code (e.g., Fuller et al. 1986; Umeda et al. 2016; Woods et al. 2017). In the case presented in Figure 8, the SMS clearly avoids collapse long after the polytropic criterion has been satisfied within the core; the accreting SMS only encounters the post-Newtonian instability and undergoes hydrodynamic collapse after reaching ${\approx}330\,000\,\mathrm{M}_{\odot}$.

Figure 8. Evaluation of the polytropic criterion for the onset of the post-Newtonian instability for a numerical simulation of the structure of an SMS accreting $10\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$. This is shown for two moments in its evolution: once it has reached $10^{5}\,\mathrm{M}_{\odot}$ (blue lines), and at ${\approx}3.2\times10^5\,\mathrm{M}_{\odot}$ (red lines), the latter being shortly before the collapse of the star (Woods et al. 2017). Solid lines plot $\Gamma_1-4/3\approx\beta/6$, and dashed lines indicate the $n=3$ polytropic criterion for instability, both as a function of mass coordinate. For comparison, the solid line for a true $n=3$ polytrope would be flat with mass coordinate, and the instability would arise once the dashed line rose above it at any point in the star.

Evidently, following the accretion history and evolution of an SMS, with the ability to follow its hydrodynamic response in detail at the onset of instability, is essential in order to accurately predict the final mass at collapse of massive black hole seeds. Several studies have now investigated this question for the case of non-rotating SMSs accreting zero-metallicity gas. For ${{\gt}}0.1\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$, the final mass at collapse is well approximated as a logarithmic function of the accretion rate:

(16)\begin{equation} {M}_\mathrm{SMS,final}\approx\left[0.83\,\log_{10}\left(\frac{\dot{\it M}}{\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}}\right)+2.48\right]\times10^5\,\mathrm{M}_{\odot} \label{eqn16} \end{equation}

although there is some divergence at higher accretion rates between different codes, rising to a factor of ${\sim} 2$ for ${\sim} 10\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$ (see Section 3.3 for further discussion). For the first time then, the seed masses of the first quasars may be described as a function of their initial conditions, albeit only for the non-rotating case. In principle, a rapidly accreting protostar should gain not only mass but angular momentum, quickly reaching very high spin velocities at the surface. In the following section, we turn to the evolution of such rotating SMSs.

3.2.2. Rotation and mass loss

Star formation always requires mechanisms to extract angular momentum from the collapsing pre-stellar gas (Spitzer 1978). The existence of stars proves that nature finds ways to circumvent this angular momentum problem, through magnetic fields, gravitational torques, or viscosity in accretion disks. The efficiency of these mechanisms, however, depends on the properties of the forming star, mainly on its mass and metallicity. SMSs are not immune to this angular momentum problem. Thus their rotation is not only of interest in itself, but is also another mechanism which could prevent SMS formation by rapid accretion. Furthermore, rotation could impact the life and death of SMSs by stabilising the star against radial pulsations (Fowler 1966) and by allowing a small fraction of the stellar mass to stay outside the horizon after the collapse (Baumgarte & Shapiro 1999b).

As a consequence of this angular momentum problem, rotating monolithic models (Section 3.1), which did not address the formation process, assumed that SMSs rotate at the Keplerian velocity, losing mass as they contract. This limit is thus called the ‘mass-shedding’ limit (Baumgarte & Shapiro 1999a). In particular, it implies a strong deformation of the star by rotation. Moreover, since monolithic models are fully convective, the star is assumed to rotate as a solid-body.

This picture changes completely when we include the accretion process. SMSs forming by accretion must be slow rotators, with surface velocities smaller than 10–20% of their Keplerian velocity (Haemmerlé et al. 2018b). At such low velocities, the deformation of the star by rotation is negligible. This constraint is a consequence of the $\Omega\Gamma$-limit (Maeder & Meynet 2000): for stars near the Eddington limit, like SMSs, the critical velocity at which the effective gravity vanishes at the equator is reduced compared to the Keplerian limit, due to the contribution of radiation pressure to hydrostatic equilibrium. The constraint from the $\Omega\Gamma$-limit on the surface velocity of accreting SMSs is illustrated in Figure 9. If at a given stage the surface velocity violates this constraint, accretion stops and the star cannot reach higher masses. Thus SMSs formed by accretion never reach the mass-shedding limit.

Figure 9. Constraint from the $\Omega\Gamma$-limit on the surface rotation velocity of accreting SMSs, as a function of their mass. The blue curve indicates the upper velocity limit imposed by the $\Omega\Gamma$-limit (Haemmerlé et al. 2018b).

Moreover, in contrast to the monolithic models, accreting SMSs rotate highly differentially, with a frequency in the core 4–5 orders of magnitude higher than that at the surface (Haemmerlé et al. 2018b). Indeed, accreting models are mostly radiative and contract in a highly non-homologous way (Section 3.2). Angular momentum transport in radiative regions turns out to be negligible on such short timescales, and non-homology enhances differential rotation in the case of local angular momentum conservation.

For the constraint on surface velocity to be satisfied, the angular momentum of the accreted material must be smaller than about 1% of the Keplerian angular momentum (Haemmerlé et al. 2018b). If we assume accretion to proceed via a Keplerian disk, it implies that 99% of the angular momentum in the inner disk must be extracted. The mechanism by which the accretor gets rid of this angular momentum excess is still unknown. Three of the most common mechanisms invoked in understanding this problem for modern star-formation have also been studied in the context of Pop III SMSs and are thought to be sufficiently efficient:

1. – magnetic fields (Latif & Schleicher 2016);

2. – gravitational torques (Begelman et al. 2006; Wise, Turk, & Abel 2008);

3. – viscosity in the disk (Takahashi & Omukai 2017).

In particular, the efficiency of gravitational torques and viscosity grows as the mass density in the disk increases. Thus if at any stage the constraint from the $\Omega\Gamma$-limit is violated, accretion stops and mass accumulates in the disk, until angular momentum transport becomes efficient enough for accretion to start again. In this picture, the upper velocity limit shown in Figure 9 is not only an upper limit, but the actual velocity expected for SMSs of various masses. Nevertheless, understanding the specific role of each of these mechanisms in the context of SMS formation requires progress in numerical simulations of the accretion process, including a self-consistent treatment of the inner accretion disk, which currently limits our knowledge on this question.

Due to their low surface velocities, accreting SMSs do not lose mass through rotational mass-loss. As well, radiative mass-loss is inefficient because of the low opacity on the Hayashi limit, which results in winds that never reach the escape velocity (Nakauchi et al. 2017). The only mechanism allowing for mass-loss is pulsational instability (Inayoshi, Hosokawa, & Omukai 2013; Hosokawa et al. 2013). However, the mass-loss rate due to pulsations is not expected to exceed about a few 0.001M$\odot$ yr$^{-1}$, i.e., several orders of magnitude smaller than the accretion rate. Thus mass-loss does not prevent SMS formation by accretion.

3.3. Benchmarking results

Numerical studies of the evolution of rapidly accreting SMSs (Section 3.2) agree on several conclusions:

1. – SMSs accreting at the rates predicted for atomically cooled halos ($0.1$–$10\,\mathrm{M}_\odot\,\mathrm{yr}^{-1}$) evolve as red supergiant protostars, along the Hayashi limit. As a consequence, their ionising feedback has negligible impact on accretion.

2. – They collapse at masses of several $10^5\,\mathrm{M}_\odot$, during central H-burning, through the GR instability. The higher the rate, the larger the final mass.

3. – The structure of accreting SMSs is mostly radiative, with a convective core due to H-burning.

Aside from this general agreement, however, questions remain regarding

1. – the exact value of the final mass (and especially reducing the uncertainty at high accretion rates);

2. – the size of the convective core at given mass (in particular at collapse);

3. – the central hydrogen abundance at a given mass.

Quantitative comparisons are given in Table 2 and Figure 10.

Table 2. Properties of SMSs at collapse for various accretion rates, according to Umeda et al. (2016) (U+16), Woods et al. (2017) (W+17), and Haemmerlé et al. (2018a) (H+18).

Figure 10. Final masses of accreting SMSs as a function of their accretion rate, from Umeda et al. (2016), Woods et al. (2017), and Haemmerlé et al. (2018a).

The origin of these divergences is not clear. The differences in the final masses may arise from the differences in the size of the convective core, since the larger the convective core, the closer the structure is to an $n=3$ polytrope, and thus the earlier the instability arises. The differences in the size of the convective core could be due to the treatment of convection and atmosphere, the accretion of entropy, or the implemented opacities. In particular, the various studies used different prescriptions for convection (e.g., use of Schwarzschild vs. Ledoux criterion, inclusion of overshooting, see Kippenhahn et al. 2012; Ledoux 1947) relying on the free parameters intrinsic to the mixing-length theory. Thus a precise knowledge of the structure of rapidly accreting SMSs would require one to go beyond mixing-length theory.

The accretion of entropy (i.e., the thermal boundary condition) could also impact the structure of the star. The simplest prescription is ‘cold-disc accretion’, where the entropy of the accreted material matches continuously that of the outermost layer of the star. This is the lowest limit for accretion of entropy, corresponding to a disk-like accretion geometry, where any entropy excess can be radiated in the polar directions before being advected. On the other hand, some codes include advection of entropy through various prescriptions (Hosokawa et al. 2013; Woods et al. 2017). The advection of entropy, however, is thought to be small compared to the energy budget of the star, due to the high luminosity at the Eddington limit (Hosokawa et al. 2013).

3.4. Discussion and outstanding problems

From the preceding discussion, it is clear that the last decade has seen accelerating progress in understanding the structure and fate of massive black hole seeds; however, many questions remain. Aside from the outstanding discrepancies between numerical solutions outlined in the previous section, a problem of particular interest is the maximum mass attainable for DCBH seeds prior to collapse. Does such a limit exist? From Section 3.2.1, it is clear that the answer depends directly on the maximum accretion rate attainable within a DCBH-forming halo. Given the logarithmic dependence of the final mass at collapse on the accretion rate, however, it appears to be difficult to construct DCBH seeds with initial masses greater than ${\sim}10^{6}\,\mathrm{M}_{\odot}$. A related, and perhaps more worrying question, is the extent to which a massive seed black hole can continue growing without merger events in an isolated halo (e.g., Tanaka & Haiman 2009) Though still uncertain, progress on this front may come soon from the search for IMBH relics in the low-redshift Universe (see Section 5.2).

Closely related to the problem of SMS evolution and collapse is the viability and nature of quasi-stars (Begelman 2010). These objects form from the collapse to a black hole of a relatively smaller seed mass, surrounded by a massive, giant-like envelope. In this case, the pressure needed in order to support this envelope is provided by the radiation released by accretion at the base of the envelope, in a manner analogous to Thorne–Żytkow objects (Thorne & Żytkow 1975). Ball et al. (2011) and Ball, Tout, and Żytkow (2012) were the first to produce detailed numerical models for the evolution of quasi-stars; however, these early efforts neglected hydrodynamics and nuclear-burning. Including these effects will be vital in any complete model for the evolution of quasi-stars. Fiacconi and Rossi (2016, 2017) investigated the structure of rotating quasi-stars and found no stable solutions for truly supermassive objects. At the moment, it is also unclear what accretion history would naturally produce a quasi-star. These questions should be investigated further.

Another outstanding issue emerging from numerical investigations of SMS evolution is the possibility of ‘supermassive supernovae’. While a relativistic collapse may be reversed for $Z\mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}$0.005 (Fuller et al. 1986), an SMS cannot be formed through the atomically cooled halo scenario at such high metallicity (though recall the discussion in Section 2 of other possibilities). Conversely, for primordial metallicities, only an extremely narrow range of masses were found to explode in the monolithic case (Chen et al. 2014), and further studies will be necessary in order to verify whether this result is robust against small changes to the input physics. If SMSs can indeed undergo thermonuclear explosions, they would provide a remarkable nucleosynthetic signature in the early Universe. Note, however, that the elemental abundances of extremely metal-poor stars in the Milky Way Galaxy (Umeda & Nomoto 2003; Ishigaki et al. 2014) and damped Lyman α systems (DLAs) at high redshifts (Kobayashi, Tominaga, & Nomoto 2011) suggest that chemical enrichment there was driven more predominantly by core-collapse supernovae from ${\sim}$ 10-40 M_⊙ stars (Ishigaki et al. 2018). Observational prospects for investigating the deaths of SMSs and their possible mass return are discussed further in Section 5.

4. Formation rates of black hole seeds

New pathways that can lead to the formation of massive seeds are now emerging. If DCBHs can indeed form in the early Universe, the obvious question becomes how often? And where? Modelling the statistics of early quasar formation necessarily requires a careful understanding of the physics of star formation in primordial halos, with the first stars’ contribution to the build-up of Lyman–Werner radiation and metals in the early Universe setting bounds on the formation of atomically cooled halos (recall Section 2). For those black hole seeds which do form (‘massive’ or not), an immense supply of infalling matter must also then be available within a very short time frame in order to produce ${\sim} 10^{9}\,\mathrm{M}_{\odot}$ quasars by $z\sim6$, particularly as even under the most optimistic assumptions the growth of black hole seeds via mergers appears to be insufficient to account for the observations (see, e.g., Tanaka & Haiman 2009).

Increasingly, the need for available gas to feed the growth of the first quasars points to massive seeds as the likely progenitors of the most massive high-redshift quasars (e.g., Agarwal et al. 2016a; Valiante et al. 2017a; Gallerani et al. 2017; Pacucci et al. 2017c, 2018). Massive seeds have the advantage of boosting further growth, and being born in more gas-rich environments (Latif et al. 2013b; Pacucci, Volonteri, & Ferrara 2015a; Volonteri et al. 2016). Low-mass seeds may still play a role; however, many questions remain. Howmay gas be brought in from cosmological scales to the accreting black hole? And how can feedback during bursts of accretion and the response of the gas to that feedback be mitigated, allowing the seed to continue to grow?

Concurrent with the growth in our understanding of their formation, there has been an explosion of new efforts to construct semi-analytic models (SAMs) that can predict the abundance and growth of DCBH seeds (e.g., Dijkstra et al. 2008; Agarwal et al. 2012; Dijkstra et al. 2014; Agarwal et al. 2014; Habouzit et al. 2016b; Valiante et al. 2016), with halo properties derived either from Monte Carlo algorithms based on the extended Press–Schechter formalism (Bond et al. 1991) or from cosmological N-body simulations. Theoretical models now agree on the expected number density of DCBHs/SMSs, so long as the same recipes (i.e., $J_{\rm{crit}}$ threshold and metal pollution) are implemented. At the same time, however, there has been increasing pressure to move away from simplified recipes, particularly the idea of a single value of $J_{\rm{crit}}$ existing in nature for atomically cooled halos. Rather, the shape of the incident spectra and the importance of H⁻ photo-detachment, along with ${\rm{H}}_{2}$ photo-dissociation, must be accounted for in any realistic model (see discussion in Section 2). The importance of other key criteria, such as avoiding external metal pollution, has also likely been overestimated in the past (cf. Section 4.2). Going forward, further progress may require new observational constraints, such as the mass function and abundance of IMBHs (see Section 5). There is now a consensus that the number of seeds in a theoretical model should be larger than the number of $z=6$ quasars, as simulations indicate that not all seeds formed will later grow into $10^9\,\mathrm{M}_{\odot}$ BHs. If mergers between seeds are unimportant for their growth, then IMBHs found in the local Universe could trace the evolution of host galaxies that underwent only a limited number of mergers. In the following, we focus on addressing these and other questions regarding the likely abundance of black hole seeds in the Universe, whether massive or not, and summarise a number of recent results.

4.1. Light and heavy seeds

The very recent possible detection of first star formation occurring ${\sim} 180 \, \mathrm{Myr}$ after the Big Bang (Bowman et al. 2018), combined with the detection of massive quasars already formed and shining by $z \sim 7$ (Fan et al. 2006; Mortlock et al. 2011; Wu et al. 2015), indicate that the growth process to form the first super-massive black holes (SMBHs) occurred in $\mathbin{\lower.3ex\hbox{$\buildrel\lt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 600 \, \rm{Myr}$. As discussed above, these time constraints profoundly challenge our view of how black holes grow: how massive must the ‘seeds’ of the first SMBHs have been? In this context, the word ‘seed’ refers to the original black hole that, growing in mass by accretion and mergers, generates an SMBH of order ${\sim} 10^{9-10} \, \mathrm{M}_{\odot}$. In the following, we expand upon the previous discussion of the ‘seeding’ problem in order to place the ‘heavy’ seed scenario in context, by contrasting it with the viability of ‘light’ seeds.

Assuming that accretion is the main driver of growth, the standard scenario predicts that the black hole mass $M_{\bullet}$ increases exponentially, with an e-folding time scale $t_\mathrm{S} \sim 0.045\,\epsilon_{0.1} \, \mathrm{Gyr}$, named the Salpeter time [see also Eq. (1)]. Here, $\epsilon_{0.1}$ is the matter-radiation conversion factor normalised to the standard value of 10%, valid if the growth occurs at the Eddington rate. Another implicit assumption of this treatment is that there is a constant gas inflow to feed the growth. If this is not the case, a duty cycle $\mathcal{D}$ can be used to express the fraction of cosmic time during which the black hole is actually accreting. The relevant equation to describe the time evolution of the mass of the seed is thus:

(17)\begin{equation} M_{\bullet}(t) = M_{\bullet ,0} \exp\left({\frac{t}{t_\mathrm{S}}}\right)\!, \label{eqn17} \end{equation}

where $M_{\bullet ,0}$ is the mass of the black hole seed at $z \sim 20$, corresponding to a cosmic age of ${\sim} 180 \, \mathrm{Myr}$. A common terminology in the field differentiates ‘light seeds’ ($\mathbin{\lower.3ex\hbox{$\buildrel\lt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 10^{2} \, \mathrm{M}_{\odot}$) from ‘heavy seeds’ ($\mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 10^{4} \, \mathrm{M}_\odot$) (e.g., see the recent review Valiante et al. 2017a). Light seeds are formed at the endpoint of the evolution of Pop III stars, while heavy seeds are formed by a number of different processes active in the early Universe and described below.

It is very challenging to envision that the $z \sim 7$ SMBHs were formed from light seeds, if their accretion rates are indeed restricted by the Eddington limit. In fact, this growth process would require constant accretion at the Eddington rate until $z \sim 7$ to allow the formation of ${\sim} 10^9 \, \mathrm{M}_\odot$ quasars. This process, while not strictly impossible, is very unlikely, due to the stringent requirement of a steady gas reservoir with low angular momentum. In order to solve the problem, two pathways can be devised directly from the evolution equation [Eq. (17)]: either decreasing the time-scale $t_S$ or increasing the initial mass $M_{\bullet ,0}$.

To decrease the e-folding time scale, while keeping $M_{\bullet ,0}$ in the light seed regime, it is necessary to reach super-Eddington accretion rates (Volonteri & Rees 2005), often indicated with $\lambda_{\rm Edd} \equiv \dot{M}/\dot{M}_{\rm Edd} {\gt} 1$. Note that, at least in a spherically symmetric accretion scenario, this is equivalent to assuming that the radiative efficiency $\epsilon$ is lower than the standard value adopted for Eddington-limited accretion. Several works (Alexander & Natarajan 2014; Madau et al. 2014; Volonteri, Silk, & Dubus 2015) have predicted the occurrence of largely super-Eddington accretion episodes at high redshift.

The second possibility assumes that the environmental conditions (i.e., pristine or very low-metallicity gas) of the early Universe ($z \mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}$10–15) allowed alternative pathways to form massive black hole seeds: (i) the direct collapse of un-enriched and self-gravitating pre-galactic disks (Begelman et al. 2006; Lodato & Natarajan 2006, 2007), (ii) the collapse of a primordial atomiccooling halo into a DCBH (Bromm & Loeb 2003; Shang et al. 2010; Johnson et al. 2012), or (iii) the formation of a very massive star from runaway stellar mergers in a dense cluster (Omukai, Schneider, & Haiman 2008; Devecchi & Volonteri 2009; Davies, Miller, & Bellovary 2011). Alternatives to these two scenarios include the possibility that at least a fraction of $z{\gt}6$ quasars might be lensed (e.g., Wyithe & Loeb 2002; Fan et al. 2019; Pacucci & Loeb 2019). A lensing effect would in fact decrease the mass of the SMBH powering the quasar.

Currently, the general consensus is that massive seeds are the most likely progenitors of early SMBHs (e.g., Agarwal et al. 2016a; Valiante et al. 2017a; Gallerani et al. 2017; Pacucci et al. 2017c; Pacucci, Natarajan, & Ferrara 2017b; Pacucci et al. 2018), due to a combination of two factors. First, as can be clearly seen in Eq. (17), larger initial masses $M_{\bullet ,0}$ boost the growth early on (Pacucci et al. 2015a; Volonteri et al. 2016); since the process follows an exponential growth law (assuming that the black hole grows constantly), this results in a very relevant boost. Second, the environment where heavy seeds formed tend to be characterised by higher gas densities (e.g., Latif et al. 2013b; Pacucci et al. 2015a), at the cusp of the density distribution of primordial halos; this, in turn, increases the chances of feeding the black hole seed with large gas inflows early on. Several works (e.g., Pelupessy, Di Matteo, & Ciardi 2007; Alvarez, Wise, & Abel 2009) have shown that light seeds, instead, are formed in low-density regions, mainly due to the effect of stellar feedback that sweeps the surrounding gas away, causing a delay before the central density is rebuilt.

Thus, the combination of large seed mass and high density leads to a very efficient and rapid growth of heavy seeds. On the contrary, light seeds will struggle to grow, at least at very early times, given the low seed mass and the low-density environments in which they are born.

From a purely analytical point of view, it is possible to show that even light seeds can grow to the SMBH regime early on in the history of the Universe (Volonteri & Rees 2005). The main concern for this scenario, however, is the availability of a steady inflow of gas at super-Eddington rates. Mass supply rates available in high-z galaxies certainly allowed for episodes of super-Eddington accretion. Assuming a velocity dispersion σ ($\sigma_{100}$ is expressed in units of $100 \, \mathrm{km \, s^{-1}}$) for the host, the characteristic free-fall rate of self-gravitating gas exceeds the Eddington rate by a factor of (Begelman & Volonteri 2017; Pacucci et al. 2017c):

(18)\begin{equation} \dot{M}_\mathrm{ff} \sim 10^5 \sigma_{100}^3 \dot{M}_{\rm Edd} \left( \frac{M_{\bullet }}{10^6 \, \mathrm{M_{\odot}}}\right)^{-1}. \label{eqn18} \end{equation}

Indeed, simulations (e.g., Dubois, Volonteri, & Silk 2014) and analytical estimates (e.g., Wyithe & Loeb 2012) suggest that super-Eddington infall rates were quite common in the early Universe. Begelman and Volonteri (2017), for example, show that the fraction of AGN accreting at super-Eddington rates could be as high as ${\sim} 10^{-3}$ at $z=1$ and ${\sim} 10^{-2}$ at $z=2$. Directly comparing these considerations with the growth of massive seeds, however, requires more detailed calculations of their rates of formation and growth, which we explore in the following subsection.

4.2 Narrowing in on the number of DCBHs

Recent efforts to determine the abundance of DCBHs in the $z{\gt}6$ Universe have found number densities spanning several orders of magnitude, from ${\sim} 10^{-1}$ to $10^{-9}$ cMpc${-3}$. This is because these estimates have employed a large variety of models using various methods and simulation set-ups, which we outline below (see also Valiante et al. 2017b). It is necessary for models to include at least three main ingredients which lead to the physical conditions needed for DCBH formation: self-consistent star formation in a halo right from when it first appeared as a minihalo, treating metal-pollution in said halo, and local build-up of the LW radiation field emanating from the first galaxies.

As discussed earlier, exposure to strong LW radiation is one of the possible ways to achieve the low ${\rm{H}}_{2}$ fraction needed for DCBH formation (Omukai 2001; Ciardi et al. 2000; Haiman et al. 2000). Ahn et al. (2009) demonstrated that the spatial variation of LW radiation intensity was key for this scenario. Most of the models exploring DCBH abundances now include a spatially varying LW radiation flux from local irradiating sources. This radiation intensity is computed either directly from stellar particles depending on their age, distance, and redshift (Agarwal et al. 2014; Habouzit, Volonteri, & Dubois 2017), or from the stellar mass painted on DM halos (Agarwal et al. 2012; Dijkstra et al. 2014; Habouzit et al. 2016b; Chon et al. 2016). Overall, studies agree on the fact that the maximum local specific intensity from Pop III stars appears to be almost always below the critical intensity, whereas a majority of pristine halos exposed to radiation from Pop II stars see a level of $J_{\rm{crit}}$ above the critical radiation threshold. The spatial distribution of radiation intensity is in good agreement between various studies, using either position-based LW radiation modelling (Agarwal et al. 2012; Chon et al. 2016) or probabilistic determination of the radiation field (Dijkstra et al. 2008).

Another important requirement for the DCBH scenario is that the gas should be free of metals. This requirement can decrease the number of DCBH candidates by orders of magnitude (see, e.g., Figure 5 in Ferrara et al. 2014). There is a general consensus that both prior episodes of star formation in halo progenitors and galactic winds from nearby halos may play a fundamental role in governing the fate of the collapsing gas. In principle, stellar winds from neighbouring galaxies could prevent DCBH formation on scales of $ \le 10{\mkern 1mu} {\rm{kpc}}$, thereby reducing their abundance (Dijkstra et al. 2014). However, typical separation between direct collapse candidate halos and LW sources is of the same order (${\sim}10\, \rm{kpc}$, Agarwal et al. 2012; Habouzit et al. 2016b). Hydrodynamical studies of metal-mixing in halos due to external SNe winds have found, however, that the extent of the impact is often over-estimated in large simulations. In practice, metal-pollution in primordial gas is complicated by a number of considerations, including the rate at which winds expel metals from their original galaxies, the propagation of metals through the Intergalactic medium (IGM), and subsequent mixing in gas of varying density (Cen & Riquelme 2008; Maio et al. 2019; Smith et al. 2015). Recently, Agarwal et al. (2017c) devised a semi-analytical model which explores the worst case scenario for DCBH formation, given metal-pollution driven by the same galaxies which provide the necessary LW flux. They report that even at an extremely short distance, ${\sim}300 \rm{pc}$, and with instantaneous metal mixing, the irradiating galaxy is not able to sufficiently pollute the DCBH candidate halo (Omukai et al. 2008; Latif et al. 2016c) to prohibit DCBH formation.

In general, there are two approaches employed by the community in order to follow the evolution of any proto-galaxy, using analytic prescriptions depending on their dark matter halo:

1. – pure semi-analytic models (pSAMs): which employ analytic algorithms (e.g., Monte Carlo) based on the extended Press-Schechter formalism (EPS Press & Schechter 1974; Lacey & Cole 1993). This approach allows one to statistically capture the abundance and mass assembly history of halos in the Universe; however, the trade-off is a lack of halo properties such as positions and spin, and thus such an approach relies on a probabilistic determination and implementation of physical effects.

2. – hybrid semi-analytic models (hSAMs): that employ cosmological N-body simulations (e.g., Springel 2005) to extract DM halo properties (e.g., mass and spatial distribution) and superimpose their analytic models. While such an approach provides physical properties such as positions, spins, and virial radius of halos at any given redshift, resolution limitations do not allow one to create a robust statistical sample of $z\sim 6$ quasar hosts while resolving the smallest Pop III host halos at the same time.

Differences between pSAMs employing PS merger trees and hSAMs using N-body dark matter only simulations (or hydrodynamical simulations) arise mainly due to the physical processes that can be included in the models. pSAMs assign probabilities to physical properties, such as the probability for dark matter halos to have a given stellar mass, of being star-forming at a given time or being metal-free. This is because pSAMs are designed to present a statistical view of the first galaxies and cannot capture the inter-dependence of halo assembly histories. On the other hand, hSAMs are able to follow the individual histories of halos and galaxies, and capture the underlying interdependence. Therefore whether a halo is star-forming is naturally captured, and the local radiation intensity or metal pollution can be consistently computed from the position-based distribution of halos produced by the simulations.

Dijkstra et al. (2008, 2014) and Valiante et al. (2016) are studies that use their own pSAMS. Dijkstra et al. (2008) compute the probability distribution function of DM halos that are exposed to the expected LW flux on the basis of their clustering properties. They found that only a small fraction, ${{\lt}}10^{-6}$, of all atomic cooling halos are exposed to a sufficient LW flux ($J_{\rm LW}{\gt}10^3$), producing a number density of ${{\lt}}10^{-6}$ $\rm{cMpc}^{-3}$ potential DCBHs hosts. In a later study, Dijkstra et al. (2014) estimate the abundance to be $n_{\rm DCBH}\sim 10^{-9}$–$10^{-6}$ cMpc${-3}$ between $z=20$ and 7, and explore the effects of model assumptions on their results, including the LW photon escape fraction and $J_{\rm{crit}}$, underlying the role of galactic winds in decreasing the formation rate of DCBH sites. In their pSAM designed to explore the formation of a $z\sim 6$ SMBH J1148, Valiante et al. (2016) predict a mean number density of ${\sim} 10^{-7}$ cMpc${-3}$ DCBHs at $z{\gt}6$. Although several DCBH seeds form in the progenitors of the $10^{13}\,\mathrm{M}_\odot$ DM host halo of J1148, only a fraction of these seeds eventually end up contributing to the mass budget of the final SMBH.

Agarwal et al. (2012, 2014) and Habouzit et al. (2016b) are some examples of hSAMs. Agarwal et al. (2012) report a higher number density in the range $10^{-2}$–$10^{-1}$ cMpc${-3}$ for $J_{\rm{crit}}=30$. In their recent work Habouzit et al. (2016b) found that the number density of DCBH sites can lie anywhere in the range $10^{-7}$–$10^{-2}$ cMpc${-3}$ depending on the value of $J_{\rm{crit}}$ imposed. If a higher critical flux is imposed for DCBH formation ($J_{\rm{crit}} {\gt} 100$), then Dijkstra et al. (2014) and Habouzit et al. (2017) find a lower DCBH abundance that is just sufficient to reproduce the population of quasars. In other words, the abundance of DCBH sites is inversely proportional to the value of $J_{\rm{crit}}$ imposed in the models. The hSAMs self-consistently track each halo’s mass assembly history (and thus the star formation) and compute the build-up of the global and local LW radiation flux on the basis of the positions of halos at any given redshift. Studies done using smaller simulation boxes that resolve minihalos and allow for a more developed chemical network have led to higher DCBH site abundances. This can be particularly attributed to the lower $J_{\rm{crit}}=30$ (Agarwal et al. 2012, 2014) which in turn suggests that the DCBH scenario may be able to seed the SMBHs of more normal galaxies. A recent study by Habouzit et al. (2017) argues that this strongly depends on the implementation of Supernova (SN) feedback; weak SN feedback may better explain BHs in normal galaxies, however, simulations with strong (SN) feedback produce galaxies and BHs in better agreement with observed statistics.

Alternative direct collapse scenarios, such as the galaxy merger-driven model, have also been studied with the aid of hSAMs applied to the Millenium simulation (Bonoli, Mayer & Callegari 2014; MB18). These models, that apply several constraints on the properties of host galaxies of direct collapse seeds, still end up with a predicted number density of bright QSOs at $z \sim 6$ that is a factor of few higher than the observed one (MB18).

Habouzit et al. (2017) perform a comparison between the pSAM study of Dijkstra et al. (2014) and a variety of hSAMs (Agarwal et al. 2012; Agarwal et al. 2014; Habouzit et al. 2016b), and find that Dijkstra et al. (2014) overestimates the stellar mass in halos. Dijkstra et al. (2014), however, underestimate both the number of galaxies contributing to the LW flux and the size of metal-polluted regions. Interestingly, the effect of these varying assumptions can sometimes be to compensate for each other, providing similar results for the number density of viable DCBH sites (Habouzit et al. 2017). The studies discussed in this section provide upper limits on the abundance of DCBHs, with a consensus emerging that the number density of DCBHs is greater than the observed number density of quasars at $z\sim 10 \ \rm{cGpc^{-3}}$ (e.g., Bañados et al. 2016). However, in their current formulation, they may yet be unable to explain the presence of central black holes in galaxies like the Milky Way and dwarf galaxies, which have much higher number densities, ~0.1 cMpc${-3}$.

5. Observations past and future

A significant motivation for the many theoretical developments outlined in the preceding sections has been the remarkable discovery of ${\sim} 10^9\,\mathrm{M}_{\odot}$ quasars in the early Universe, more of which continue to be found (e.g., Mortlock et al. 2011; Wu et al. 2015; Matsuoka et al. 2018; Bañados et al. 2018). The challenge now is to find unique, unambiguous signatures for the formation of these objects which may be realistically sought after immediately or in the near future. In particular, the need now is to distinguish between the formation of massive black holes from low and high mass seeds, either by finding distinct signatures of the formation history of particular observed quasars at $z\sim$7, directly tracing the evolution of related progenitors at even higher redshift with next generation telescopes, or searching for evidence of either formation scenario in the local Universe.

The Luminosity Function (LF) of quasars at high redshift can provide an essential probe of the growth of SMBH (e.g., Tanaka & Haiman 2009). Ongoing efforts such as the VISTA surveys (Sutherland 2009), as well as future observations with the Large Synoptic Survey Telescope (Ivezić et al. 2008) and Euclid (Refregier et al. 2010) will probe the evolution of the quasar LF over cosmic timescales, eventually reaching out to $z\sim$ 7–8 (Manti et al. 2017; LSST Science Collaboration et al. 2009). Constraining the number density and accretion efficiencies of very high-z quasars in this way may soon make it possible to assess the viability of light and heavy seed models.

New diagnostics have also been proposed. For the massive seed scenario, X-ray ionisation by very early massive black holes may leave a distinct imprint on the 21 cm background at $z \mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 20$ (Tanaka, O'Leary, & Perna 2016), and the sites of DCBH formation may reveal themselves by a unique 3 cm (restframe) emission line (Dijkstra, Sethi, & Loeb 2016a). The first direct, unambiguous evidence of any Pop III stars may come from observations of their explosive deaths as supernovae at high-redshift (Hummel et al. 2012; Magg et al. 2016; Hartwig, Bromm, & Loeb 2018a). Owing to the high effective temperatures expected for Pop III stars, it has long been argued that strong He ii recombination line emission can provide a clear signature for their presence within a stellar population (e.g., Schaerer 2003; Visbal, Haiman, & Bryan 2015; Cai et al. 2015). He ii emission alone, however, is not a reliable diagnostic. Without additional line measurements, a He ii nebula’s ionising source cannot be unambiguously distinguished from, e.g., an AGN (Bowler et al. 2017b), and conversely, even if a stellar population produces a strong He II-ionising flux, this may not produce a correspondingly high He ii line luminosity if the ionisation parameter is sufficiently low (see, e.g., discussion in Woods & Gilfanov 2013, Appendix 2).

The difficulty in interpreting high-z emission line spectra is clearly demonstrated by the case of CR7 (Sobral et al. 2015). Initially thought to be the first definitely Pop III-like galaxy ever observed (though for earlier claims see, e.g., Malhotra & Rhoads 2002; di Serego Alighieri et al. 2008), this bright $z = 6.6$ Lyman-α emitter has provided a critical test of our ability to discriminate between theoretical interpretations of high-z galaxy observations. In the following, we begin with a focused overview of the rapid developments and lessons learned in studying CR7, before discussing alternative means to either directly or indirectly probe the formation and growth of black holes in the early Universe. We then discuss how, in the near future, the search for IMBHs may yield new insight into the ‘initial mass function’ (IMF) of black hole seeds, while the JWST and other upcoming telescopes may reveal how they grow with time, or even reveal SMSs themselves (or their explosions!), while future gravitational wave observatories may reveal their growth via mergers.

5.1. What have we learned from CR7?

The CR7 system was first identified as a Lyman-α emitter by Matthee et al. (2015) and spectroscopically confirmed by Sobral et al. (2015). It is the brightest Lyman-α emitter at $z=6.6$, with $L_{\mathrm{Ly}\alpha} \simeq 10^{44}\,\mathrm{erg}\,\mathrm{s}^{-1}$ and $L_\mathrm{He{ii}}=10^{43.3 \pm 0.4}\,\mathrm{erg}\,\mathrm{s}^{-1}$ (Sobral et al. 2015). More remarkable than the sheer brightness of this source was the absence of any metal line detection and the remarkable equivalent width (EW) of Heii (EW$_{\mathrm{HeII}} = 80 \pm 20$ Å, Sobral et al. 2015), indicating a hard ionising spectrum, and suggesting the presence of Pop III stars (Tumlinson & Shull 2000; Oh, Haiman, & Rees 2001). Sobral et al. (2015) indeed found that the SED of the complete system could only be reproduced when including a young (${{\lt}}5\,\mathrm{Myr}$) population of massive metal-free stars, with a total stellar mass of at least $10^7\,\mathrm{M}_\odot$ for one of the three clumps. This was one of the most compelling observations of a Pop III-like galaxy and has therefore triggered many theoretical interpretations and further observations.

By re-analysing the VLT/X-Shooter data of Sobral et al. (2015), Shibuya et al. (2018) found only upper limits on the Heii luminosity with EW$_{\mathrm{He}\,{{ii}}} {\lt} 60$ Å and no other detected emission lines in the UV. In a more recent re-analysis of the X-Shooter data, Sobral et al. (2017) confirmed the detection of He ii at lower significance than before (${\sim} 3\sigma$) with $L_{\mathrm{He}\,{\rm{II}}} \simeq 10^{43.4 \pm 0.5}\,\mathrm{erg}\,\mathrm{s}^{-1}$. Using HST/WFC3 grism observations and a re-analysis of the VLT data, Sobral et al. (2017) find a EW$_{\mathrm{He}\,{\rm{II}}} {\lt} 10$ Å, and a gas metallicity of $0.05 {\lt} Z/{\rm Z}_\odot {\lt} 0.2$. They also report that there is no convincing evidence for strong variability (${{\gt}}0.3$ mag) of any of the components of CR7, which would be a possible smoking gun for an accreting BH.

Bowler et al. (2017b) use Spitzer/Infrared Array Camera observations and interpreted the blue IRAC colour as a [O iii] and H$\beta$ detection (see also Bowler et al. 2017a; Pacucci et al. 2016). They advocate that CR7 could be a low-mass, narrow-line AGN or a young, low-metallicity starburst with a metallicity of $Z \sim 1/200\,{\rm Z}_{\rm{\odot}}$. As discussed in more details below, Pacucci et al. (2017a) showed that this updated photometry can still be explained by a DCBH model. In fact, while a zero-metallicity environment is required at formation, subsequent star formation can certainly provide low levels of metals even in the early phases of DCBH accretion.

Based on recent ALMA observations, Matthee et al. (2017) reveal a clear [C ii] signature within multiple components of CR7. They do not find any FIR continuum emission due to dust, which may also disfavour an AGN. They argue that SFRs of 5–28 $\mathrm{M}_\odot\,\mathrm{yr}^{-1}$ in the individual clumps, together with a metallicity in the range $0.1 {\lt} Z/{\rm Z}_\odot {\lt} 0.2$, can account for the observations.

In addition to the original interpretation of CR7 as a Pop III-like galaxy, Sobral et al. (2015) also proposed the possibility of a DCBH. This scenario is favourable over other AGN formation channels because of its natural absence of metals during the formation process. Soon after the initial discovery, theoretical models employing hydrodynamical simulations (Pallottini et al. 2015) and a semi-analytical model (Agarwal et al. 2016b) indicated that CR7 could indeed host a DCBH. In particular, clumps B and C could provide the required LW flux for clump A to have formed a DCBH, which would have the required observational features about 100 Myr after its formation.

Separately, Hartwig et al. (2016) showed that CR7 cannot host the necessary $10^7\,\mathrm{M}_\odot$ of Pop III stars. Chemical feedback from previous SNe prevents the accumulation of sufficient metal-free gas. An alternative scenario by Visbal, Haiman, & Bryan (2016) demonstrates how photoionisation feedback could delay the collapse of a pristine halo so that it can accumulate enough gas for a massive Pop III burst, in agreement with the original observations of CR7.

The velocity offset between Lyman-α and Heii, and the extent of the Lyman-α emitting region imply a source lifetime of at least 10 Myr, in contradiction with a very young population of massive stars (Smith, Bromm, & Loeb 2016). Smidt et al. (2016) demonstrated, using Enzo (Bryan et al. 2014) and radiative transfer in post-processing, that an AGN powered by a black hole of mass $3 \times10^7\,\mathrm{M}_\odot$ accreting at $f_\mathrm{Edd} = 0.25 - 0.9$ can account for the observations. Additional radiative processes have been suggested to prevent pristine gas from collapsing for long enough to accumulate the required high masses, but it is challenging to reproduce the required $10^7\,\mathrm{M}_\odot$ of Pop III stars (Xu et al. 2016; Visbal et al. 2016; Visbal, Bryan, & Haiman 2017; Yajima & Khochfar 2017).

In the light of metal line detections in CR7, Agarwal et al. (2017b) present an improved analytical model. If a galaxy is externally polluted with metals after the formation of a central DCBH, such a system is within 3σ of the new observations. Also Pacucci et al. (2017a) present a new radiation-hydrodynamic model that can reproduce the observations by Bowler et al. (2017b) with a DCBH of mass $M_\mathrm{BH} = 7 \times 10^6\,\mathrm{M}_\odot$, a gas metallicity of $Z=5\times 10^{-3}$, and a line-of-sight column density of neutral hydrogen of $3 \times 10^{24}\,\mathrm{cm}^{-2}$.

Although CR7 has turned out to be less outstanding than initially thought, it has triggered many important discussions with results that will remain valid, regardless of the nature of CR7. Various mechanisms can delay the collapse of pristine gas and only strong photoionisation feedback can result in a Pop III stellar population more massive than $10^6\,\mathrm{M}_\odot$ (Visbal et al. 2016). Thanks in part to the intense scrutiny given to CR7, we now have a much greater understanding of the conditions under which a DCBH can form, such as the window of opportunity to provide enough LW flux without polluting the halos (Hartwig et al. 2016; Agarwal et al. 2016b, 2017b), and what observational signatures can distinguish a DCBH from other energetic sources at high redshift, such as the 3 cm fine structure line (Dijkstra, Gronke, & Sobral 2016b). A number of alternative diagnostics have now been proposed which may confirm or refute the DCBH scenario, both for individual high-z objects and as a fraction of all SMBHs formed, which we address in the following sections.

5.2 What do IMBHs tell us about the seeds of the first SMBHs?

Intermediate-mass BHs (IMBHs) of $10^{2}$–$10^{5}$ $\mathrm{M}_{\odot}$ are thought to be the local relics of those high-z seed BHs that did not become supermassive. They should be found in nearby dwarf galaxies with low mass, low metallicity, and which have not significantly grown through accretion and mergers, hence resembling the first galaxies formed in the early Universe (Volonteri 2010). As a direct trace of high-z black hole formation, then, which kind of IMBH masses should we expect for different SMBH seed scenarios, and how frequent are IMBHs in dwarf galaxies?

Figure 11. From Ferrara et al. (2014). Schematic view of the scenario for intermediate mass black hole formation and growth accounted for in the cosmological model by Ferrara et al. (2014). The IMF of IMBHs derived by the authors model all of these different physical processes, thus taking into account both IMBH formation via truly direct collapse and via an SMS phase.

Ferrara et al. (2014) investigated the IMF of direct collapse BHs along with the mass distribution of their hosting halos. Figure 11 illustrates the possible paths for direct collapse formation and growth, which have been accounted in their model. To study the final outcome and frequency of such distinct paths the authors used a statistical merger tree approach, which includes metal-pollution from both the inhomogeneously enriched IGM and the star-forming progenitor halos (Salvadori et al. 2014). Their results show that the IMF of DCBHs extends over a large range of masses, $M_{\rm{BH}}\approx (0.2-20) \times 10^5\,\mathrm{M}_{\odot}$, and it is bi-modal. Further, they find that the hosting halos of DCBHs have low masses, $M_\mathrm{h} \approx 10^{7.5}$–$10^8\,\mathrm{M}_\odot$ at $z\approx (7-15)$. In the local Universe, therefore, un-evolved metal-poor dwarf galaxies that have masses $M_\mathrm{h} {\lt} 10^9\,\,\mathrm{M}_{\odot}$ are thus the key objects to search for the living relics of such IMBHs. Yet, given the strict conditions required for IMBHs to form from direct collapse (see Figure 11 and Section 2), a low black hole occupation fraction is expected in today’s dwarf galaxies if direct collapse was the dominant seeding mechanism.

Alternatively, SMBH seeds can form as remnants of massive Pop III stars. The observed properties of ancient dwarf galaxies in the Local Group are all consistent with the idea that these small systems hosted Pop III stars and that were likely among the first star-forming galaxies (Salvadori & Ferrara 2009; Salvadori, Skúladóttir, & Tolstoy 2015). These theoretical expectations has been reinforced over the years by chemical abundance studies of ancient stars in the newly discovered ultra-faint dwarf galaxies, which revealed the chemical imprint of Pop III stars with masses $M_*\approx$ (10–60) $\mathrm{M}_{\odot}$ (e.g., more recently Spite et al. 2018). Cosmological simulations for the formation of the first stars show that Pop III stars can also be very massive, $M_*\approx 10$–$1\,000\,\mathrm{M}_{\odot}$ (Hosokawa et al. 2011; Hirano et al. 2014, 2015), a prediction that is consistent with the Pop III IMF constraints inferred with Milky Way halo stars (de Bennassuti et al. 2017). In conclusion, black hole seeds from Pop III stars likely have masses $M_{\rm{BH}}\approx10$–$1\,000\,\mathrm{M}_{\odot}$ and their living relics may be found in present-day dwarf galaxies (Salvadori & Ferrara 2009; Salvadori et al. 2015; Hartwig et al. 2015; Magg et al. 2018, 2019). Since Pop III stars are more common than direct collapse BHs, these lighter seeds should yield a higher black hole occupation fraction in present-day dwarf galaxies (Volonteri, Lodato, & Natarajan 2008; van Wassenhove et al. 2010).

Several studies have therefore focused on finding IMBHs actively accreting as AGN in dwarf galaxies. A few hundred of such low-mass AGN, with black hole masses $\mathbin{\lower.3ex\hbox{$\buildrel\lt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 10^{6}$ $\mathrm{M}_{\odot}$, have been found using either emission line diagnostics or the width of broad emission lines from the gas, assumed to be virialised, around them (e.g., Greene & Ho 2004, 2007; Dong, Greene, & Ho 2012b; Moran et al. 2014; Reines, Greene, & Geha 2013; Sartori et al. 2015; Chilingarian et al. 2018; Nguyen et al. 2019). The AGN nature of some of these sources has been reinforced by the detection of hard X-ray emission (Dong et al. 2012a; Baldassare et al. 2015, 2017). The use of deep X-ray surveys such as Chandra COSMOS Legacy (Civano et al. 2016) has provided a few more candidates (e.g., Schramm et al. 2013; Pardo et al. 2016; Mezcua et al. 2018a) as well as observational evidence that a population of IMBHs must exist in dwarf galaxies, at least out to $z\sim$1.5 (Mezcua et al. 2016). Their detection is challenging, however, due to their low luminosity and possibly mild obscuration (Mezcua et al. 2016). The redshift record-holder of a low-mass AGN in a dwarf galaxy is cid_1192, at $z=2.39$ (Mezcua et al. 2018a), which reinforces the scenario in which dwarf galaxies host the relic seed BHs of the early Universe. Such deep X-ray surveys have additionally allowed us to derive the AGN fraction in dwarf galaxies, taken as a lower limit to the local black hole occupation fraction (Pardo et al. 2016; Mezcua et al. 2018a).

Even when correcting for completeness, the low AGN fraction found (i.e., 0.4% for $$z\leq$$ 0.3; Mezcua et al. 2018a) and its decrease with stellar mass (Aird, Coil, & Georgakakis 2018; Mezcua et al. 2018a) seem to favour the direct collapse formation scenario for seed BHs (Volonteri et al. 2008; van Wassenhove et al. 2010; Bellovary et al. 2018). This scenario is also favoured by the finding that the black hole mass scaling relations flatten in the low-mass regime (Greene, Ho, & Barth 2008; Jiang et al. 2011; Graham & Scott 2013, 2015; Mezcua 2017; Martín-Navarro & Mezcua 2018; Shankar et al. 2019), as expected from simulations (Volonteri et al. 2008; van Wassenhove et al. 2010). Whether this flattening is an observational bias (Mezcua 2017), the result of dwarf galaxies being dominated by supernova feedback (Martín-Navarro & Mezcua 2018), or the existence of a bimodality in the accretion efficiency of BHs (Pacucci et al. 2018) remains unclear.

In addition to dwarf galaxies, some IMBH candidates have also been found in ultraluminous X-ray sources (Farrell et al. 2009; Mezcua et al. 2011; Webb et al. 2012; Mezcua et al. 2013, 2015, 2018b), possibly as the remnant core of a dwarf galaxy that was stripped in the course of a minor merger. We note that mergers of dwarf galaxies are found to be very frequent and could lead to BH coalescence and rapid BH growth of the primordial seed (e.g., Deason, Wetzel, & Garrison-Kimmel 2014; Mezcua 2019). The effects of AGN feedback could also significantly impact BH growth, by either shutting down or triggering star formation and thus the BH fodder (e.g., Dashyan et al. 2018; Mezcua et al. 2018b; Penny et al. 2018; Regan et al. 2018; Mezcua 2019). If IMBHs in dwarf galaxies have significantly grown through mergers and feedback, then they should not be considered as the fossils of the seed BHs of the early Universe (Mezcua 2019). This has important implications for seed BH formation studies that are based on local IMBHs. Evidence for IMBHs could also be found from tidal disruption events (TDEs) (Lin et al. 2017), high-velocity clouds (Oka et al. 2016, 2017), and in globular clusters (Gebhardt, Rich, & Ho 2005; Lützgendorf et al. 2013; Kızıltan, Baumgardt, & Loeb 2017), though in these cases they could have formed locally hence their relation with the high-z seed BHs is dubious (see Mezcua 2017, for a review). Finally, we note that virtually all mass measurements (e.g., Greene & Ho 2004, 2007; Dong et al. 2012b; Moran et al. 2014; Reines et al. 2013; Chilingarian et al. 2018) rely on virial mass estimates based on empirical line/continuum relations for AGN in the optical part of the spectrum, under the assumption that the emission lines of interest (e.g., Hα or Hβ lines) originate in the vicinity of an active nucleus (i.e., broad or narrow line AGN). Validating and improving these calibrations are areas of ongoing study (Shankar et al. 2019), which we discuss further in the following subsection.

5.3. Searching for IMBHs in low luminosity AGN in the multimessenger era

The discovery of hundreds of SMBHs in the $$\mathbin{\lower.3ex\hbox{$\buildrel\lt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}}10^{6}\rm\,{M}_{\odot}$$ range in the centres of nearby low-mass AGN (see previous section) has set low-luminosity AGN (LLAGN) as prime targets in the search for IMBHs. With the vast majority of candidate IMBHs in nearby low-mass AGN lying in the $$10^{5}-10^{6}\rm\,{M}_{\odot}$$ range (see recent review by Mezcua 2017), these measurements tentatively suggest a prevalence of the direct collapse scenario, indicated by the presence of primarily heavy seeds. However, it is very likely that these findings reflect the limitations of our observational capabilities, as well as a bias towards more luminous and therefore more massive candidates in the LLAGN sample. In addition to the selection bias of more luminous host galaxies targeted for spectroscopy in the SDSS sample (see, e.g., discussions by Greene & Ho 2004 and Barth, Greene, & Ho 2008), the spectral and spatial resolution of telescopes (like the 5 m telescope used in the Palomar survey) lies right at the threshold of detecting the signatures of a slowly accreting black hole with a mass on the order of a few $$10^{5}\rm\,{M}_{\odot}$$ . This is in terms of both the size of the broad line region and the broadening of emission lines (see also Koliopanos 2017 and references therein).

To investigate the fainter end of LLAGN and dwarf galaxies, one can draw from the correlations between large-scale properties of galaxies and the mass of their central SMBH. One of the most well-known scaling relations is that between black hole mass and stellar bulge velocity dispersion ( $$M_\mathrm{BH}{-}{\sigma}$$ : e.g., Ferrarese & Merritt 2000; Gebhardt et al. 2000). While this relation exhibits a small scatter and is very well established, it is luminosity biased (as more luminous and massive galaxies are selected to populate the relation) and requires the spectroscopic analysis of the host bulge emission. Furthermore, there are considerable doubts on the applicability of the M _BH-σ relation in the low mass regime. Most models that aim to reproduce the M _BH-σ relation, through the interaction of massive outflows (during episodes of super-Eddington accretion) with the host galaxy, predict a size for the region of influence that for $$M_\mathrm{BH}\,{{\lt}}\,10^{5}\rm\,{M_{\odot}}~$$ is substantially smaller than what can be resolved by current or near-future optical telescopes (e.g., Silk & Rees 1998; King 2003).

To bypass some of these drawbacks and also access fainter sources, the scaling relation between $$M_\mathrm{BH}$$ and luminosity of the host bulge (L) can be considered. Actually preceding the discovery of the $$M{\text{-}}{\sigma}$$ relation, the $$M{\text -}L$$ relation (e.g. Dressler 1989; Kormendy 1993) can be applied to fainter sources and spectroscopy is not required. While the $$M{\text -}L$$ relation is also biased towards higher mass/luminosity, recent works have attempted to better constrain the low mass regime, discovering strong indications of a steeper slope in the $$M{\text -}L$$ relation, at $$M_\mathrm{BH}\,{{\lt}}\,10^{6}\,\mathrm{M}_\odot$$ (Graham & Scott 2013). Using these improved scaling relations, Graham and Scott (2013) presented 40 lower luminosity spheroids that contain AGN which appear to have a central black hole mass that falls in the IMBH range. Other scaling relations between host galaxy and central black hole include a correlation between $$M_{\rm{BH}}$$ and the galactic Sérsic index ( $$n_{\rm{sph}}$$ , a measure of bulge concentration see, e.g., Graham & Driver 2007; Savorgnan et al. 2016), and a correlation between the morphology of the spiral arms and the mass of the central black hole in disk galaxies (Seigar et al. 2008; Davis et al. 2012) formulated as a relation between $$M_\mathrm{BH}$$ and the spiral pitch angle (PA, a measure of the tightness of spiral arms).

Figure 12. Left panel: spectra for an SMS accreting at 1.0 M_⊙ yr⁻¹ at 100 000 yr (Haemmerlé et al. 2018a). Blue: spectrum of the star itself; red: spectrum of the star in the dense accretion envelope that creates it. Right panel: AB magnitudes for the SMS in its envelope in the JWST F250 (2.5 μm, black), JWST F356 (3.56 μm, green), JWST F444 (4.44 μm, blue), and JWST F460 (4.60 μm, red) filters.

The different scaling relations each have their benefits and drawbacks, which may be more notable in the under-explored low mass regime (Shankar et al. 2019). However, combining the different methods provides a more robust average prediction. Recently, Koliopanos et al. (2017) combined IR, optical, and FUV observation of LLAGN and the M-L, M-PA, M-n _sph relations, along with X-ray and radio observations of the accreting central black hole and employing a well-known relation between the radio luminosity ( $$L_{r}$$ , 5 GHz), the X-ray luminosity ( $$L_{x}$$ , 0.5–10 keV), and the black hole mass (the ‘fundamental plane of black hole activity’: Merloni, Heinz, & di Matteo 2003), to estimate the SMBH mass of candidate LLAGN and probe the consistency between the relations in the low mass regime. The authors demonstrated that the different methods yielded consistent results in the low mass regime. Prompted by the promising results of this pilot study, an endeavour has begun to extend it to a complete census of the mass of the central SMBHs in nearby AGN (within 150 Mpc). The census cross-correlates the Third Reference Catalogue of Bright Galaxies (RC3; de Vaucouleurs et al. 1991), supplemented with the Catalogue of Neighbouring Galaxies (Karachentsev et al. 2004) and the AGN sample with measured central velocity dispersions (Ho et al. 2009; Xiao et al. 2011), from the Palomar survey (Ho et al. 1997) and the 4th Data Release of the Sloan Digital Sky Survey (Greene & Ho 2007), with the XMM-Newton and Chandra catalog and the VLA and e-Merlin radio surveys of nearby AGN, in search of candidates with multi-wavelength signatures of the accreting central black hole. So far, more than 50 candidates—that have available X-ray, radio, FUV, optical, and IR observations from a slew of major telescopes (e.g., XMM-Newton, Chandra, VLA, HST, Spitzer, GALEX, SDSS, KPNO, etc.)—have been identified (Koliopanos et al., in preparation). When completed, this endeavour promises to yield one of the most robust and comprehensive mass measurements of IMBHs in LLAGN to date.

In the era of multi-messenger astronomy, this type of multi-wavelength study is a necessary step towards the ongoing search for IMBHs and the seeds of SMBHs. By scrutinising the remarkable wealth of observational data that are already available in numerous astronomical catalogs, we will also provide the target sample for the next generation of observatories such as the JWST, Athena X-ray telescope, and the Square Kilometre Array (SKA) project. The identification of the seeds to SMBHs is already within reach and promises to be one of the great discoveries of the next decade.

5.4 Direct detection of SMSs with future observatories

5.4.1 Can we observe SMSs themselves?

As discussed above, rapidly accreting SMSs may be the precursors to the first quasars in the Universe. Could this initial stage of early SMBH formation be detected in the Near infra-red (NIR) by upcoming wide-field surveys such as Euclid and the Wide-Field Infrared Survey Telescope (WFIRST), or in narrower but deeper surveys by JWST and the coming generation of 30–40 m telescopes on the ground? SMSs typically have luminosities of $$\sim10^6\rm\,{L_{\odot}}$$ at the onset of central hydrogen burning that later rise to $$\sim10^{10}\rm\,{L_{\odot}}$$ by the time the star collapses to a black hole. There are two main challenges to their detection today.

First, although they are extremely luminous they are also relatively cool, with surface temperatures of 5 000−10 000 K, set by H⁻ opacity in their outer layers. At these temperatures, they are not intrinsically bright sources of UV radiation that would be redshifted into the NIR today. However, a few SMSs accreting at rates below $${\sim}0.005\, \mathrm{M}_\odot\,\mathrm{yr}^{-1}$$ (Haemmerlé et al. 2018a) or intermittently, in bursts (Sakurai et al. 2015), may evolve along hotter (bluer) tracks at times and be brighter today.

Second, these stars are only bright for relatively short durations of 10⁵−10⁶ yr, and do not appear as transients. However, although current evolution models of SMSs cannot resolve pulsations, they are in fact expected to exhibit them, and this might enhance prospects for their detection in future surveys. Depending on the frequency of pulsations in the rest frame, variations in the luminosity of the star might be manifested over the cadences of surveys currently planned for Euclid, WFIRST, JWST, as well as ground-based ELTs in the near future.

Source frame spectra for an SMS accreting at 1 $$\rm{M_{\odot}}\,yr^{-1}$$ before and after reprocessing by its dense accretion envelope are shown in the left panel of Figure 12 (Surace et al. 2018). Absorption by the envelope at the H and He ionisation edges is evident, as is Lyα emission. Continuum absorption due to H⁻ opacity is also visible at wavelengths shorter than 1.65 μm. The reprocessing of shorter wavelength photons by the envelope of the star into Lyα enhances prospects for their detection, since at z ~ 10 these photons are redshifted in the NIR.

We show NIR AB magnitudes for the SMS in the right panel of Figure 12. Unfortunately, predicted SMS H band AB magnitudes are well below the sensitivities of both Euclid and WFIRST at z ~ 10–20. SMS could, however, still appear in JWST and ELT surveys out to z ~ 20 in the 3.56 and 4.44 μm NIRCam filters. The likelihood for detection of SMSs is likely greatest at z ~ 10, as their numbers are thought to peak at this epoch and they are more easily seen at lower redshifts.

5.4.2 Do SMSs explode?

To date, no model of rapidly accreting, zero-metallicity SMSs has produced an SN; however, it has recently been suggested that monolithically formed SMSs within a relatively narrow (5 000 M_⊙) window around 55 500 M_⊙ (Chen et al. 2014) do end their lives in extremely energetic (~10⁵⁵ erg) thermonuclear explosions, capable of destroying the protogalaxies in which they occur (Whalen et al. 2013a, 2013b, Johnson et al. 2013). Numerical studies have found that such explosions would be visible to Euclid and WFIRST at the epochs at which they would likely occur (Figure 13; Whalen et al. 2013c). Furthermore, some models of rotating stars above $$10^6\rm\,{M_{\odot}}$$ (Montero, Janka, & Müller 2012) have been found to explode, as have metal-enriched (0.005 Z_⊙) non-accreting ~500 000M_⊙ stars (Fuller et al. 1986). While it is unclear at the moment how such massive (55 000–500 000 M_⊙) stars may form monolithically, the abundance patterns predicted for their explosions (e.g., Chen et al. 2014) suggest that these must be rare events; in particular, they are inconsistent with observations of the Milky Way and DLAs (Kobayashi 2017, see also discussion in Section 3.4). Interestingly, however, mass return and chemical enrichment from an SMS undergoing putative eruptive mass loss during its lifetime on the hydrogen-burning main sequence may help to explain the unusual abundances observed in globular clusters (Denissenkov & Hartwick 2014; Gieles et al. 2018), although a number of other scenarios exist.

Figure 13. Lower limits on the NIR light curves of the thermonuclear explosion of a 55 500 M_⊙ star. Each light curve is shown in the filter in which the explosion is brightest at the given redshift. Note that the wavelength of the filter limits the redshift at which a transient can be detected in it because the IGM absorbs all flux blueward of the Lyman limit in the rest frame of the event prior to the end of cosmological reionisation. Thus, Euclid and WFIRST can only detect transients in the H band out to z ~ 15 while JWST and the ELTs could detect events at z ≳ 20.

Matsumoto et al. (2015, 2016) examined the observational signatures of highly energetic but directional gamma-ray bursts (GRBs) driven by the collapse of SMSs. Such events would be characterised by central engines with unusually long durations of 10⁴ - 10⁵ s or more because of the times required for the outer layers of the star to collapse into the black hole accretion disk system (BHAD) at its centre, and for its jet to pierce these massive, infalling layers. Matsumoto et al. (2016) find that the gamma-ray signal of SMS GRBs would be detected by Swift and its afterglow would be visible to Euclid, WFIRST, and JWST in the NIR. However, to drive a GRB the core of the SMS would have to be spun up to very high angular momenta prior to collapse, in order to produce a BHAD capable of driving a jet for the times required to break out of the star. The high angular velocities required to form a BHAD in the core approach the breakup limit in stars, and are thought to be why the rate of GRBs is only 10⁻³–10⁻⁴ of all core collapse events today. The SMS could not be born with such angular velocities because the $$\Omega\Gamma$$ instability would have broken up the star before it could become supermassive, and no known mechanism could spin it up at late times prior to its collapse (Haemmerlé et al. 2018b). In any event, the SMS GRB rate is expected to be only a small fraction of the rate of DCBH formation, and the detection rate would be smaller by roughly a factor of $$1/4\pi$$ because they are beamed events, making their detection even by wide-field NIR surveys problematic.

If the earliest quasars also formed from low mass black hole seeds via super- or hyper-Eddington accretion, their birth might have been signalled by less massive but still highly energetic Pop III GRBs. The gamma-ray signal from Pop III star GRBs could in principle be found by Swift even at $$z\mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 20$$ (e.g., Natarajan et al. 2005; Bromm & Loeb 2006; Nomoto et al. 2010; Mészáros & Rees 2010; Nagakura, Suwa, & Ioka 2012), but the time dilation of the burst might not immediately identify it to be a high redshift event because of the large range of central engine times that have now been found for GRBs in the local Universe. Prompt spectroscopy of the NIR afterglows of these events in future missions such as EXIST (Gou et al. 2004; Whalen et al. 2008; Toma, Sakamoto, & Mészáros 2011; Wang et al. 2012; Mesler et al. 2012, 2014) or by radio observations from the ground (Ioka & Mészáros 2005) could pinpoint their redshifts. Their event rate, while low, would still be expected to be much higher than those of SMS GRBs (e.g., Bromm & Loeb 2002; de Souza, Yoshida, & Ioka 2011).

Finally, we note that fragmentation in the vicinity of a DCBH deep in the centre of a massive accretion disk of an atomically cooled halo could lead to energetic transients known as TDEs. TDEs occur when stars forming from such fragments are torn apart by tidal forces just before crossing the event horizon of the DCBH. In this picture, X-rays from the nascent black hole may also trigger the formation of Pop III stars in its vicinity within 5–10 Myr of birth (e.g., Machacek, Bryan, & Abel 2003; Aykutalp et al. 2014). Statistically, some of these stars would then be scattered onto the loss cone of the black hole within ~1 Myr. Numerical studies indicate that the jet emitted in this event could be detected by Swift and eROSITA even at z ~ 20 (Kashiyama & Inayoshi 2016). The X-ray signal from TDEs could be distinguished from those of DCBH GRBs by their much longer durations, $$10^{5 - 6} (1+z)$$ s. Detections of the X-ray signal could then be followed up by radio or NIR observations of the afterglow with the extended Very Large Array and the SKA or JWST to pinpoint its redshift of origin. The rate of TDEs would be on par with that of DCBH formation itself, much higher than the rate of SMS SNe or GRBs, and would therefore be the likeliest electromagnetic transient associated with DCBH birth.

5.5 Observability of DCBHs

To begin this section, it is important to point out that DCBHs are still black holes: as such, we can expect them to have the typical spectral signatures of black holes. If this is true, how can we possibly distinguish them from normal high-z AGN? In fact, once formed, a black hole does not carry a label indicating which formation channel it followed. It is through the interaction of the black hole with its galactic environment that we might be able to understand if an accreting object possibly formed as a DCBH. Differently from other formation mechanisms, DCBHs are predicted to form in a pristine environment (e.g., Bromm & Loeb 2003; Agarwal et al. 2014; Volonteri et al. 2016). The smoking gun for a DCBH would, thus, be the observation of a black hole spectrum with the absence of any metal line. As the case of CR7 (described in Section 5.1) taught us, however, the detection of very feeble metal lines does not rule out the DCBH interpretation (Agarwal et al. 2016b; Pacucci et al. 2017a). In fact, while the conventional model for the formation of a DCBH requires a metal-free gas, once the black hole is in place the gas can certainly become polluted with metals resulting from star formation. This would influence the resulting spectrum of the DCBH. As pointed out in Pacucci et al. (2017a), the DCBH would still be the preferred formation model to explain the spectrum of a high-z black hole seed if the metal content of the host galaxy is $$\mathbin{\lower.3ex\hbox{$\buildrel\lt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 10^{-4} \, \mathrm{Z_{\odot}}$$ , where Z_⊙ indicates the solar metallicity.

The first studies, both analytical (Yue et al. 2013) and numerical (Pacucci et al. 2015b), to predict the spectrum of a DCBH indicated predominant emission in the observed IR (≳1 μm) and X-ray ( $${\sim} 1 \, \mathrm{keV}$$ ) part of the electromagnetic spectrum. Photons with frequencies higher than Lyman- $$\alpha$$ are absorbed by photo-electric processes and re-emitted in the IR. Figure 14 shows a classic DCBH spectrum as computed in Pacucci et al. (2015b). Accretion on a typical DCBH (initial mass $$10^5 \, \mathrm{M}_{\odot}$$ in a zero-metallicity environment) lasts $${\sim} 100$$ Myr. The figure shows that most of the accretion process on a typical DCBH should be easily observable by the JWST (see also Natarajan et al. 2017; Volonteri et al. 2017), with peak magnitudes even reaching $$m_{\rm AB} \sim 20$$ for a very short period of time ( $${\ll} 1 \, \mathrm{Myr}$$ ). The HST should be able to observe their emission up to $$z \sim 10$$ close to the Lyman- $$\alpha$$ line. In the X-ray part of the spectrum, the CDF-S is predicted to observe the emission of a typical DCBH close to the peak of the high-energy emission, for a short period of time (~10–20 Myr). The same result should be obtainable by ATHENA with a much shorter ( $${\sim}300$$ ks) integration time, while Lynx will be able to easily observe most of the accretion process (Vito et al. 2016; Vikhlinin 2018).

Figure 14. Time evolution of the spectrum of a classic DCBH with initial mass 10⁵ M_⊙ in a pristine (z=0) environment. The infrared, optical, and X-ray bands are highlighted with shaded regions, while the unprocessed spectrum is reported, at peak luminosity (t=115 Myr), with a dashed line. The flux limits for future (JWST, ATHENA, and Lynx) and current (HST and CDF-S) observatories are also shown. Image from Pacucci et al. (2015b).

The two main environmental factors that can change the shape of the spectrum for DCBHs are: (i) the presence of metals, and (ii) very large (Compton thick, $$N_H \gg 10^{24} \, \mathrm{cm^{-2}}$$ ) absorbing column densities. The presence of a non-negligible ( $$Z \mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} {10^{ - 5}}{\mkern 1mu} {{\rm{Z}}_ \odot }$$ ) amount of metals in the galactic environment of a DCBH produces a shift of photons from the X-ray to the IR part of the spectrum, due to Auger’s processes on high-energy electron shells (see, e.g., Pacucci et al. 2017a). Increasing the absorbing column density of the gas (even at $$Z=0$$ ) has a similar effect, with a decrease of the X-ray emission followed by an increase of IR re-emission. The effects of an increase in metallicity/column density are very similar, but not degenerate since the presence of metals also produces metal lines in the spectrum (Agarwal et al. 2016b; Pacucci et al. 2017a). The effect of an increase of the column density above Compton-thick levels is particularly interesting for DCBH studies, since most of the high-z sources are predicted to be heavily obscured (e.g., Comastri et al. 2015). The spectral predictions for heavily obscured DCBHs (Pacucci et al. 2015b; Valiante et al. 2018) indicate that the detection of these sources could be challenging in high-energy bands, while in the IR JWST will be able to observe them. Thus, the JWST will certainly be the main observatory to unravel highly obscured sources in the high-z Universe (Natarajan et al. 2017; Valiante et al. 2018). Very large values of the absorbing column density may also lead to photon-trapping effects, with a reduced radiative efficiency of the accretion disk and, consequently, very large accretion rates, possibly above the Eddington limit (Sa¸dowski 2009; Begelman & Volonteri 2017). As pointed out in Pacucci et al. (2015b), this would reduce the accretion time scale on to the DCBH by a factor $$\sim 10$$ , with typical evolutionary times as low as $${\sim}10 \, \mathrm{Myr}$$ . This effect would certainly add to the challenge of detecting highly obscured sources in the early Universe.

Even if DCBHs are observable in the early Universe, possibly up to z ~ 20–25 with the JWST (Natarajan et al. 2017), another interesting question is how we can separate them from other high-z sources. The availability of high-resolution spectra, fundamental to confirm the detection of a DCBH (see the discussion in Gallerani et al. 2017; Pacucci et al. 2017a), is limited only to a small fraction of the sources observed, mostly the brighter ones. DCBHs, instead, are expected to be found in the high-z Universe, where the number of sources with spectroscopic information is of course limited. Deep, wide-field surveys (e.g., the CANDELS survey, see Koekemoer et al. 2011; Grogin et al. 2011) contain photometric information about hundreds of thousands of sources, thus offering an invaluable amount of information to find sources whose photometry is at least compatible with the one theoretically predicted for DCBHs. Once a source is selected as a DCBH candidate, follow-up spectroscopic observations (if possible) would eventually determine their actual nature. The first photometric method to identify black hole seed candidates in deep multi-wavelength surveys was developed in Pacucci et al. (2016). Supported by numerical simulations, this work predicted that DCBHs are characterised by a steep spectrum in the IR (1.6–4.5μm), i.e., by very red colours. A comparison with other high-z categories of sources clearly shows that the IR spectrum of DCBHs is predicted to be significantly steeper. The method selected the only 2, X-ray detected, objects found in the CANDELS/GOODS-S survey with a photometric redshift $$\mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 6$$ . To date, the selected objects represent the most promising black hole seed candidates, possibly formed via the DCBH scenario, with predicted mass $${\gt}10^5 \, \mathrm{M}_{\odot}$$ .

Several works have complemented these selection criteria (Natarajan et al. 2017; Volonteri et al. 2017; Valiante et al. 2018) adding, for instance, constraints on the ratio of X-ray flux to rest-frame optical flux. In particular, upon their birth, DCBHs have BH mass to halo mass ratios much higher than expected for BH remnants of Pop III stars which have grown to the same mass ( $${\sim}10^6\,\mathrm{M}_\odot$$ ). A future X-ray mission, such as Lynx, combined with IR observations should be able to distinguish high-redshift DCBHs from smaller black hole seeds, due to the correspondingly high BH mass to stellar mass ratios of the faintest observed quasars, with inferred BH masses below $${\sim} 10^6\,\mathrm{M}_\odot$$ (Visbal & Haiman 2018).

While the development of improved selection criteria will definitely help in developing a reliable catalogue of DCBH candidates, only the advent of the next generation of telescopes (e.g., JWST, ATHENA, Lynx) will likely lead to the detection of the first, confirmed, DCBH.

5.6 Future insights from gravitational waves

The gravitational wave detector LIGO has proven its capabilities by detecting binary black hole mergers with total masses of a few times $$10\,\mathrm{M}_\odot$$ (Abbott et al. 2017). In the near term, future observations hold substantial promise for furthering our understanding of not only stellar-mass black holes, but their more massive cousins as well. At design sensitivity, LIGO may yet detect IMBH binaries (either IMBH+stellar-mass black hole or of IMBH-IMBH) of $${\leq} 2 \times 10^3\,\mathrm{M}_\odot$$ (e.g., Fragione, Ginsburg & Kocsis 2018). Future interferometers are designed to detect even more massive black hole mergers, such as KAGRA (Aso et al. 2013) or LISA (Babak et al. 2017) at lower frequencies. The expected design of LISA has the highest sensitivity in the mass range 10⁴–10⁶ M_⊙, for which black hole mergers are detectable out to $$z{\gt}20$$ . This is the same mass range for remnants of SMS that may form in binaries, as recent high-resolution simulations suggest (Latif et al. 2013c; Becerra et al. 2015; Latif et al. 2016b; Chon et al. 2018; Regan & Downes 2018a,b). Depending on the stellar evolution of supermassive binaries and the initial separation of their compact remnants, these objects will merge and can be observed using LISA.

Recently, Hartwig et al. (2018b) have calculated the detection rates of such black hole binary (BHB) mergers and the resulting constraints on the nature of SMS formation. Since there are many unknowns, such as the binarity of SMS, their stellar binary evolution, and the corresponding delay times until coalescence, this ignorance can at present only be parameterised using the scaling factor $$f_\mathrm{bin}$$ , which quantifies the number of SMS remnant BHBs per pristine atomic cooling halo that will merge on a sufficiently short time scale. In their most optimistic model with $$f_\mathrm{bin}=1$$ and $$J_{\rm{crit}}=30$$ , Hartwig et al. (2018b) show that we can expect up to $${\sim} 80$$ mergers per year of the remnants of supermassive binary stars. However, a major challenge will be to discriminate these events from other formation channels, such as black hole mergers as a consequence of galaxy mergers. To break this degeneracy, Hartwig et al. (2018b) compare their anticipated detection rates to the results of Sesana et al. (2007), who derive the rates of black hole mergers in the same mass range as a consequence of galaxy mergers. Hartwig et al. (2018b) demonstrate that black hole mergers from supermassive stellar binaries dominate the LISA detection rates at $$z \mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 15$$ , with about $$0.6$$ detections per year in their optimistic scenario (see Figure 15).

These results also allow one to determine what constraints on the abundance of supermassive stellar binaries can be drawn from a non-detection of such black hole mergers. In Hartwig et al. (2018b), their optimistic model ( $$f_\mathrm{bin}=1$$ , $$J_{\rm{crit}}=30$$ ) can already be excluded with 95% certainty after 5 yr of non-detections at $$z{\gt}15$$ . After 10 yr of non-detections, an upper limit on the binarity of SMSs of $$f_\mathrm{bin}{\lt}0.5$$ ( $$J_{\rm{crit}}=30$$ ) can be derived at the $$2\sigma$$ level.

Further refinements of this work, as well as predictions for other observables (e.g., X-ray and UV emission), await a more refined theory of the formation and evolution of SMS multiples. Given the expected prevalence of fragmentation based on recent numerical simulations (e.g., Chon et al. 2018; Regan & Downes 2018a, 2018b), understanding binary interactions in the supermassive regime is a priority.

Figure 15. Rates for the merging of black hole binaries (BHB) from supermassive stellar binaries as a function of redshift compared to the models by Sesana et al. (2007). Only the optimistic scenario with $$f_\mathrm{bin}=1$$ and $$J_{\rm{crit}}=30$$ can produce a population of BHB mergers at $$z \mathbin{\lower.3ex\hbox{$\buildrel\gt\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 15$$ that are clearly distinguishable from other channels of BHB formation. The total rate of such uniquely identifiable BHB mergers is $${\sim} 0.6$$ per year, highlighted in grey. Adapted from Hartwig, Agarwal, & Regan (2018b).

6. Conclusions

The discovery of massive high-z quasars has provided invaluable insight, and a great challenge, to the study of SMBH formation. Given the difficulty in reaching such masses so early in cosmic history, their existence has motivated a rapidly growing effort to understand the formation of much more massive ‘seed’ black holes than could be produced by ‘typical’ Pop III stars. A viable mechanism for the formation of truly supermassive ( $${{\gt}}10^5\,\mathrm{M}_{\odot}$$ ) seeds has emerged in the synchronised pairs model for atomically cooled halos, with an increasingly mature theory for the conditions under which they may arise. A number of alternative models, however, continue to challenge this scenario.

All of these models, in turn, provide realistic initial conditions for detailed stellar evolutionary calculations, permitting a clearer picture than ever before of the nature and fate of SMSs. Many questions remain, however, regarding the mechanism by which SMSs may shed angular momentum, and their fate at the highest plausible accretion rates.

The state of the field has advanced, however, to the point that one can begin to consider ‘subgrid’ models for SMSs. The statistics of SMS formation can also now be addressed; however, the essential conditions for atomically cooled halos to form massive black holes must be robustly determined before the viability of the heavy seed model can be clearly evaluated in this way. In the meantime, observational efforts to determine the seeds of SMBHs continue, from the search for IMBHs, to deep surveys of the distant Universe. Recent efforts have demonstrated that discriminating between light and massive seed models is feasible using multi-wavelength electromagnetic data from the planned and upcoming JWST, WFIRST, and Athena missions, as well as anticipated gravitational wave observations of high-z black hole mergers to be detected by the proposed LISA mission (Ricarte & Natarajan 2018; Pezzulli et al. 2017; Volonteri et al. 2017). In the near future, the next generation of electromagnetic and gravitational wave observatories promises to revolutionise our understanding of the first stars and quasars, and will likely reveal the origin of these titans of the early universe.