On the Detection of Supermassive Primordial Stars. II. Blue Supergiants

TL;DR

This paper models the spectra and brightness of hot, blue supermassive primordial stars, suggesting they could be detected by upcoming telescopes like JWST, Euclid, and WFIRST at redshifts up to 10-12.

Contribution

It provides the first detailed spectral and photometric predictions for blue supermassive primordial stars, highlighting their potential observability.

Findings

JWST can detect these stars up to z ~ 10-12

Euclid and WFIRST could detect them with gravitational lensing

Blue supermassive stars are brighter and hotter than previously modeled red hypergiants

Abstract

Supermassive primordial stars in hot, atomically-cooling haloes at $z \sim$ 15 - 20 may have given birth to the first quasars in the universe. Most simulations of these rapidly accreting stars suggest that they are red, cool hypergiants, but more recent models indicate that some may have been bluer and hotter, with surface temperatures of 20,000 - 40,000 K. These stars have spectral features that are quite distinct from those of cooler stars and may have different detection limits in the near infrared (NIR) today. Here, we present spectra and AB magnitudes for hot, blue supermassive primordial stars calculated with the TLUSTY and CLOUDY codes. We find that photometric detections of these stars by the James Webb Space Telescope (JWST) will be limited to $z \lesssim$ 10 - 12, lower redshifts than those at which red stars can be found, because of quenching by their accretion envelopes. With moderate gravitational lensing, Euclid and the Wide-Field Infrared Space Telescope (WFIRST) could detect blue supermassive stars out to similar redshifts in wide-field surveys.

Equations

$$F(\lambda )=\frac{L\left(\frac{\lambda }{1+z}\right)}{\frac{\lambda }{1+z}(1+z)4\pi d_{\mathrm{L}}^2(z)}.$$

$$d_{\mathrm{L}}(z)=(1+z)c/{\mathrm{H}}_0{\int }_0^z\frac{1}{\sqrt{{\Omega }_{\mathrm{M}}(1+z{)}^3+{\Omega }_{\lambda }}}dz.$$

$$\sum _{\lambda }\frac{L(\lambda )}{\lambda }\Delta \lambda =L_{\mathrm{bol}}.$$

$$m_{\mathrm{AB}}=-2.5\mathrm{lo}{\mathrm{g}}_{10}\frac{{\int }_0^{\infty }F(\lambda )T(\lambda )d\lambda }{{\int }_0^{\infty }{F}_0(\lambda )T(\lambda )d\lambda }.$$

$$\frac{\mathrm{d}N}{\mathrm{d}z\mathrm{d}\Omega }={\dot{n}}_{\mathrm{SMS}}{t}_{\mathrm{SMS}}{r}^2\frac{\mathrm{d}r}{\mathrm{d}z},$$

$$r(z)=\frac{c}{H_0}{\int }_0^z\frac{\mathrm{d}{z}^{\prime }}{\sqrt{{\Omega }_{\mathrm{m}}(1+z^{\prime }{)}^3+{\Omega }_{\Lambda }}}.$$

Full text

††pubyear: 2019††pagerange: On the detection of supermassive primordial stars. II. Blue supergiants.–On the detection of supermassive primordial stars. II. Blue supergiants.

1 Introduction

Supermassive stars (SMSs) have long been the subject of analytical studies (e.g., Iben, 1963; Chandrasekhar, 1964; Fowler, 1964, 1966) and numerical simulations (e.g., Appenzeller & Fricke, 1972; Shapiro & Teukolsky, 1979; Fuller et al., 1986; Baumgarte & Shapiro, 1999; Sun et al., 2018; Butler et al., 2018). But credible scenarios for their formation have only recently been found: supermassive primordial star (SMS) formation in atomically-cooling primordial haloes at high redshifts exposed to either unusually strong Lyman-Werner (LW) UV fluxes (Latif et al., 2014b; Agarwal et al., 2016; Chon et al., 2017; Wise et al., 2019) or highly supersonic baryon streaming motions (Latif et al., 2014a; Hirano et al., 2017; Schauer et al., 2017) or the formation of stars powered by self-annihilation of dark matter (’dark stars’; Spolyar et al., 2008; Freese et al., 2008). Strong UV backgrounds or streaming motions can suppress star formation in a halo until it reaches masses of $\sim$ 107 M*⊙* and virial temperatures of $\sim$ 104 K that trigger rapid atomic cooling that leads to catastrophic baryon collapse that can build up a star at rates of up to $\sim$ 1 M*⊙* yr*-1* (Lodato & Natarajan, 2006; Wise et al., 2008; Regan & Haehnelt, 2009; Inayoshi et al., 2014; Latif & Volonteri, 2015). Such stars may have been the origin of the first quasars, a few of which have now been discovered at $z>$ 7 (Mortlock et al., 2011; Bañados et al., 2018; Smidt et al., 2018).

Stellar evolution models indicate that primordial (Pop III) stars growing at these rates can reach masses of a few 105 M*⊙* before, in most cases, collapsing to black holes (direct collapse black holes, or DCBHs; Umeda et al., 2016; Woods et al., 2017; Haemmerlé et al., 2018b; Haemmerlé & Meynet, 2019). A few non-accreting Pop III SMSs may explode as thermonuclear transients (Montero et al., 2012; Whalen et al., 2013b; Johnson et al., 2013b; Whalen et al., 2013c, 2014; Chen et al., 2014). Pop III SMSs are the leading candidates for the origin of the earliest supermassive black holes (SMBHs) because the environments of ordinary Pop III star BHs are hostile to rapid growth (Whalen et al., 2004; Alvarez et al., 2009; Whalen & Fryer, 2012; Smith et al., 2018). In contrast, DCBHs are born with much larger masses and in much higher densities in host galaxies capable of retaining their fuel supply even when it is heated by X-rays (Johnson et al., 2013a). But Pop III star BHs, in principle, could reach large masses by super- or hyper-Eddington growth if there is enough gas to fuel their rapid growth (Madau et al. 2014; Volonteri et al. 2015; Pezzulli et al. 2016; Inayoshi et al. 2016 – see Mayer et al. 2015; Mayer & Bonoli 2019 for other pathways to the formation of these quasars and Valiante et al. 2017; Woods et al. 2018 for recent reviews on the first quasars).

Most studies have found that rapidly accreting Pop III stars evolve as cool, red hypergiants along the Hayashi limit, with surface temperatures of 5,000 - 10,000 K due to H- opacity in their atmospheres, at least until they reach $\sim$ 105 M*⊙* (Hosokawa et al., 2013). Haemmerlé et al. (2018a) found that SMSs can remain cool even above these masses and reach luminosities $\gtrsim 10^{10}$ L*⊙. But Woods et al. (2017) found that SMSs evolving from similar initial conditions quickly settle onto hotter, bluer tracks with temperatures of 20,000 - 40,000 K. Haemmerlé et al. (2018a) found that Pop III SMSs accreting at low rates ($\lesssim$ 0.005 M⊙* yr;-1) also evolve along blue tracks, as may stars with clumpy accretion due to fragmentation or turbulence in the accretion disk (Sakurai et al. 2015 – but see Sakurai et al. 2016). Whether these differences are due to opacities, code physics (such as the numerical treatment of convection), accretion physics and boundary conditions, or numerical resolution remains unknown.

What are the prospects for observing blue Pop III SMSs today? Johnson et al. (2012) semi-analytically examined the spectral features of similar stars and predicted that they would be characterized by strong Balmer emission and the conspicuous absence of Ly$\alpha$ lines due to absorption by their envelopes. The source of this flux was the hypercompact H II region of the star, whose ionising UV was trapped close to its surface by the density and ram pressure of the inflow (which has also found to be true in cosmological simulations of highly-resolved atomically-cooled haloes; Becerra et al., 2018). Freese et al. (2010), Zackrisson et al. (2010a), Zackrisson et al. (2010b) and Ilie et al. (2012) modeled the spectra of hot, blue Pop III dark stars. They found that these objects could be observed today even by $8-10\,\mathrm{m}$ telescopes on the ground, primarily because of their high surface temperatures (20,000 - 30,000 K), larger masses (up to 107 M*⊙*) and longer lives (up to 107 yr; see also recent reviews by Freese et al., 2016; Banik et al., 2019). Most recently, Surace et al. (2018) calculated spectra for cool, red SMSs and found that some will be visible to the James Webb Space Telescope (JWST; Gardner et al., 2006; Kalirai, 2018) at $z\lesssim$ 20 and at $z\sim$ 10 - 12 to Euclid (Laureijs et al., 2011) and the Wide-Field Infrared Space Telescope (WFIRST; Spergel et al., 2015) if they are gravitationally lensed. Hartwig et al. (2018) also found that the relics of such stars would be uniquely identifiable with the gravitational wave detector LISA at $z>$ 15 if they form in binaries.

There are two challenges to modeling spectra for blue SMSs. First, unlike the cool, red stars in Surace et al. (2018), blue SMSs cannot be approximated as blackbodies (BBs) because they have much higher ionising fluxes due to their higher surface temperatures, and much of this flux is absorbed by their own atmospheres. Second, these stars are deeply embedded in hot, dense, accretion shrouds that reprocess flux from the star into longer wavelengths. Accurate spectra for blue SMSs require both stellar atmosphere models and radiative transfer through the envelope of the star. Such spectra are crucial to predicting detections of blue SMSs at high redshift, which would capture primordial quasars at the earliest stages of their development. Here, we calculate spectra and NIR magnitudes for hot, blue Pop III SMSs at high redshift with the TLUSTY and Cloudy codes. Our models are described in Section 2, and we discuss spectra, NIR magnitudes and detection rates for blue SMSs in Section 3. We conclude in Section 4.

2 Numerical Method

We model the atmospheres and spectra of blue SMSs with the TLUSTY code (Hubeny & Lanz, 1995) and how their accretion envelopes reprocess these spectra with the Cloudy code (Ferland et al., 2017). The emergent spectra are then cosmologically redshifted and convolved with filter functions to obtain their NIR AB magnitudes today.

2.1 TLUSTY Atmosphere Models

We consider SMSs accreting at 1.0 M*⊙* yr*-1* and 0.1 M*⊙* yr*-1* as fiducial cases. These stars were evolved in the Kepler stellar evolution code and discussed in detail in Woods et al. (2017). Surface temperatures, $T_{\mathrm{eff}}$, and bolometric luminosities for both stars over their lifetimes are shown in the upper row of Figure 1. The spectrum of the 1.0 M*⊙* yr*-1* star was calculated at 1.51 $\times$ 105 yr, about halfway through its lifetime of 2.51 $\times$ 105 yr, when it has a mass of 1.51 $\times$ 105 M*⊙* and a surface temperature $T_{\mathrm{eff}}=$ 36,963 K. The spectrum of the 0.1 M*⊙* yr*-1* star was calculated at 8.01 $\times$ 106 yr, about halfway through its lifetime of 1.63 $\times$ 106 yr, when it has a mass of 8.01 $\times$ 104 M*⊙* and a $T_{\mathrm{eff}}=$ 22,093 K. The bolometric luminosities of the two SMSs are 1.89 $\times$ 1043 erg s*-1* and 1.13 $\times$ 1043 erg s*-1*, respectively. As in Surace et al. (2018), we neglect the luminosity of the accretion shock at the surface of the stars because it is negligible at the velocities and densities of the infall onto the star (at most $\sim$ 104 L*⊙*).

The surface gravities of these stars are $\mathrm{log}(g)\approx 3.148$ and $\mathrm{log}(g)\approx 2.203$ for the 1.0 M*⊙* yr*-1* and 0.1 M*⊙* yr*-1* SMSs, respectively. TLUSTY has great difficulties converging for surface gravities as low as these, and we have therefore settled for spectra generated using TLUSTY v.205 with $\mathrm{log}(g)=3.25$ and 2.35 (i.e., offsets by $\Delta\mathrm{log}(g)\approx 0.1$ and 0.15). The stellar atmospheres are based on non-LTE, zero metallicity and primordial abundances of H and He. The resulting TLUSTY spectra have then been rescaled to match the actual bolometric luminosites of the two stars. Comparisons to zero-metallicity models with similar temperatures (but somewhat lower surface gravities) in the publicly available TLUSTY grids of Lanz & Hubeny (2003) and Lanz & Hubeny (2007) do not reveal significant problems due to these $\mathrm{log}(g)$ discrepancies, although we cannot rule out the possibility that we are slightly underestimating the ionizing flux in the case of the 0.1 M*⊙* yr*-1* model.

2.2 Cloudy Models

We use the TLUSTY spectra of both stars as the input spectra in our Cloudy models of the flux that emerges from the accretion envelopes of the stars, whose spherically-averaged density and temperature profiles are shown in the bottom row of Figure 1. They are taken from an Enzo cosmology code (Bryan et al., 2014) simulation of the collapse of an atomically cooled halo after the formation of the accretion disk that creates the star (see also Figure 2 of Surace et al., 2018). These envelope models do not account for feedback from the SMS perturbing the structure of the infalling gas, but it is not expected to be important because ionizing radiation from the star is trapped close to its surface, as we discuss in the next section. We surround the 1.0 M*⊙* yr*-1* SMS with the profile at 0.238 Myr after the formation of the disk and the 0.1 M*⊙* yr*-1* SMS with the profile at 1.738 Myr because the envelope has time to build up to higher central densities with the more slowly accreting star. These profiles are tabulated in Cloudy with 70 bins that are uniformly partitioned in log radius, with inner and outer boundaries at 0.015 pc and 927 pc. The temperatures in the envelope are set by the virialization of cosmic flows well above it rather than by radiation from the star because ionising UV from the star is confined to very small radii deep in the cloud. Since these temperatures determine to what degree the envelope is collisionally excited, and therefore how it reprocesses photons from the star, we require Cloudy to use the Enzo temperatures for the envelope instead of inferring them from the spectrum of the star.

Cloudy then solves the equations of radiative transfer, statistical and thermal equilibrium, ionisation and recombination, and heating and cooling to determine the excitation and ionisation state of the gas surrounding the star and calculate its emergent spectrum. These calculations use tables of recombination coefficients from Dere et al. (1997) and Landi et al. (2012) and ionic emission data from Badnell et al. (2003) and Badnell (2006). Each spectrum has 8228 bins that are uniformly partitioned in log $\lambda$. We convert the luminosity in each bin, $L(\lambda)=\lambda L_{\lambda}$ in erg s*-1*, to the flux density, $F_{\lambda}$ in erg s*-1* cm*-2* $\mu$m*-1*, needed to compute AB magnitudes (equations 1 - 3 in Rydberg et al., 2018) by

[TABLE]

Here, $\lambda$ is the wavelength in the observer frame and $d_{\mathrm{L}}(z)$ is the luminosity distance:

[TABLE]

This is done to conform to the Cloudy convention that

[TABLE]

AB magnitudes, $m_{\mathrm{AB}}$, in specific filters are then calculated from

[TABLE]

Here, $T(\lambda)$ is the filter transmission function and $F_{0}(\lambda)=3.630781\times 10^{-20}c\lambda^{-2}~{}\mathrm{ergs}~{}\mathrm{cm}^{-2}~{}\mathrm{s}^{-1}~{}\mathrm{\mu m}^{-1}$, the reference spectrum for AB magnitudes. We assume cosmological parameters from the second-year Planck release: $\Omega_{\mathrm{M}}=0.308$, $\Omega_{\Lambda}=0.691$, $\Omega_{\mathrm{b}}=0.0223$, $h=0.677$, $\sigma_{8}=0.816$, and $n=0.968$ (Planck Collaboration et al., 2016). Flux blueward of Ly$\alpha$ in the rest frame of the star is set to zero in the AB magnitude calculation because of absorption by the neutral IGM at $z>$ 6.

3 Blue Supermassive Stars

3.1 Stellar Spectra

We compare TLUSTY spectra for the blue 0.1 M*⊙* yr*-1* and 1.0 M*⊙* yr*-1* stars to those of blackbodies at the same temperatures and luminosities in the upper panels of Figure 2. In both cases the atmosphere of the star has little effect on its spectrum at wavelengths redward of its blackbody peak except for some relatively weak emission and absorption lines, but the picture is different at shorter wavelengths. The sharp drops in luminosity at 504 Å in the 0.1 M*⊙* yr*-1* star and at 227 Å in the 1.0 M*⊙* yr*-1* star are due to the ionisation limits of He I and He II, respectively. There is virtually no absorption due to hydrogen just blueward of its ionisation limit except for small features at 912 Å because most of it has been been ionised by the star. There is a weak Ly$\alpha$ line at 1216 Å and H$\alpha$ and weak Paschen lines at 6560 Å, 12800 Å, and 18,800 Å.

3.2 Reprocessed Spectra

We show Cloudy spectra for the two stars after absorption and re-emission by their accretion envelopes in the bottom row of Figure 2. The ionizing UV fluxes of the 0.1 and 1.0 M*⊙* yr*-1* stars are trapped at radii of 0.033 pc and 0.02 pc in our Cloudy models. These are the resolution limits of the Enzo model at 0.238 and 1.738 Myr, so as expected the strong inflows quench the ionizing UV of both stars. Strong continuum absorption due to ionisation of neutral H is evident below 912 Å, with additional steps in absorption at 504 Å and 227 Å due to the ionisation of He I and He II, respectively. These features are stronger with the 0.1 M*⊙* yr*-1* star than the M*⊙* yr*-1* star because its envelope has collapsed to higher central densities by 1.786 Myr. Strong H$\alpha$ and Paschen absorption lines are visible at 6560 Å, 12800 Å, and 18800 Å. Most absorption blueward of the Lyman limit is re-emitted as the continuum and numerous lines at wavelengths above 5000 Å.

In contrast to the red stars in Surace et al. (2018), absorption and re-emission by the accretion envelopes of blue stars do not enhance their spectra in most of the bands that would be redshifted into the NIR today. A potential exception is Ly$\alpha$: in contrast to red SMSs, both spectra exhibit very strong Ly$\alpha$ emission lines that are pumped by the much higher UV fluxes of the blue stars. Although these lines are strong it is not clear how much of this Ly$\alpha$ flux would be observed in the NIR today, for two reasons. First, many of the Ly$\alpha$ photons would be resonantly scattered into a halo of large radius but low surface brightness, so the star might not appear to be a strong point source of this flux. Second, repeated resonant scatterings broaden the line over time so some of the flux in principle could fall outside a given filter after being redshifted into the NIR today. This is not expected to be a large effect because the Ly$\alpha$ photons are only scattered at most $\sim$ 3% from line center before their optical depth in the wings falls below unity and they stream freely through the universe (Smith et al., 2015). If the maximum displacement of the photon from line center is 0.03$\lambda_{0}$ $=$ 36 Å in the rest frame it would be $\sim$ 0.04 $\mu$m for a $z=$ 10 SMS, or about an order of magnitude smaller than the typical width of JWST wide band NIR filters. A detailed treatment of Ly$\alpha$ radiative transfer in the primordial IGM is beyond the scope of this paper so we calculate AB magnitudes for the stars with and without the Ly$\alpha$ line as upper and lower limits.

3.3 NIR Magnitudes

NIR magnitudes for the 1.0 M*⊙* yr*-1* blue SMS in JWST, Euclid and WFIRST bands are plotted in Figure 3. We consider three cases: i) stars with accretion envelopes but no Ly$\alpha$ line; ii) stars with accretion envelopes and their Ly$\alpha$ lines; iii) stars with no envelopes. This latter case is in the event that ionising UV radiation from the star or other dynamical effects such as gravitational torqueing from nearby haloes strip away the envelope of the star. At 2.5 - 4.6 $\mu$m the magnitudes with and without the Ly$\alpha$ line are indistinguishable out to $z=$18, when it begins to be redshifted into the 2.5 $\mu$m JWST NIRCam filter, leading to an increase in brightness of about two magnitudes. A similar effect is visible in this filter with the 0.1 M*⊙* yr*-1* star in Figure 4 but is less prominent because the Ly$\alpha$ line from its envelope is weaker.

From $z\sim$ 14 - 20 the envelope of the 1.0 M*⊙* yr*-1* star somewhat suppresses flux from the star but enhances it at $z>$ 13, especially at $z<$ 7 where reprocessed radiation redward of 5000 Å enhances brightnesses by 3 - 4 magnitudes. Similar enhancements are evident with the 0.1 M*⊙* yr*-1* star at the same redshifts. Absorption in the NIR by the denser envelope of the 0.1 M*⊙* yr*-1* star is more severe, decreasing its brightness down to $z\sim$ 7. In contrast, reprocessing of the spectra redward of 5000 Å by the envelopes of both stars makes them more visible at 7.7 - 25.5 $\mu$m at nearly all redshifts, but their magnitudes remain well below MIRI detection limits. NIRCam AB magnitude limits of $\sim$ 31 will effectively limit detections of the 1.0 and 0.1 M*⊙* yr*-1* stars by JWST to $z\sim$ 12 and 10, respectively.

The 1.0 M*⊙* yr -1 star is brighter by $\sim$ 2 magnitudes in the Euclid and WFIRST filters with the Ly$\alpha$ line than without it at redshifts over which it is shifted into these filters, as shown in the lower panels of Figure 3. The 0.1 M*⊙* yr*-1* star is about one magnitude brighter. Exclusion of this line results in brightnesses that are consistently lower than those for stars without envelopes, and this effect is especially pronounced at lower redshifts where quenching by the envelope is greatest. Quenching at low redshifts is greatest with the 0.1 M*⊙* yr*-1* star because it has the denser envelope. The AB magnitudes of stars with envelopes but no Ly$\alpha$ never rise above 30 and could not be directly detected at $z\gtrsim$ 6 with Euclid or WFIRST, whose practical detection limits are 26 and 28, respectively.

3.4 SMS Detection Rates

The number of SMSs per unit redshift per unit solid angle at a redshift $z$ is

[TABLE]

where $\dot{n}_{\rm SMS}$ is the SMS formation rate per unit comoving volume, $t_{\rm SMS}$ is the average lifetime of an SMS, which we take to be 1 Myr, and $r(z)$ is the comoving distance to redshift $z$,

[TABLE]

Current estimates of $\dot{n}_{\rm SMS}$ vary by up to eight orders of magnitude (Woods et al., 2018), and these models also predict a variety of evolution of $\dot{n}_{\rm SMS}$ with redshift. Habouzit et al. (2016) find that the comoving number density of SMSs rises with decreasing $z$ whereas Valiante et al. (2017) predict that most SMSs form in the narrow range $z\sim$ 16 - 18.

As in Surace et al. (2018) we consider upper and lower limits for $\dot{n}_{\rm SMS}$. The upper limit is the low $J_{\rm crit}$ model of Habouzit et al. (2016), in which most SMSs form at $z\sim$ 10 - 12 and the final comoving SMS number density is $\sim$ $10^{-1}\,{\rm Mpc^{-3}}$. The lower limit is found by assuming that most SMSs form at $z\sim$ 16 - 18, as in Valiante et al. (2017), with a final comoving number density of $\sim$ $10^{-8}\,{\rm Mpc^{-3}}$. The upper limit yields about $4\times 10^{7}$ SMSs per steradian per unit redshift, or around 30 per NIRCam field of view. The lower limit on $\dot{n}_{\rm SMS}$ yields only $\sim 10$ SMSs per steradian per unit redshift, or at most $10^{-5}$ SMSs per NIRCam field of view. There is also some uncertainty in SMS detection rates due to their range of lifetimes, but it is small compared to the uncertainty in $\dot{n}_{\rm SMS}$.

At present, the relative numbers of blue and red SMSs are not known. Although most studies so far have found rapidly accreting SMSs to have extended red envelopes, the codes used to model their evolution lack detailed radiation hydrodynamics and opacities and can only approximate convective mixing, all of which can have profound effects on the structure of the star. Neither blue nor red SMSs have been found in the Hubble Ultra Deep Field to date because its AB mag limit is 29 at 1.38 $\mu$m, well below that expected of either type of star even at $z\sim$ 6. Strategies for the direct detection of SMSs by JWST, Euclid and WFIRST are now under development (Whalen et al., 2019).

4 Discussion and Conclusion

In contrast to cooler, redder SMSs that can be found at $z\sim$ 18 - 20, detections of hot, blue SMSs by JWST will be limited to $z\lesssim$ 10 - 12 due to quenching by their accretion envelopes. Likewise, these stars cannot be directly detected by Euclid or WFIRST at $z\gtrsim$ 6. This does not mean that these two missions cannot find blue SMSs because only moderate gravitational lensing is required to boost their fluxes above their detection limits. Their fields of view will enclose thousands of massive galaxies and galaxy clusters, and at $z\sim$ 6 - 12 magnification factors of only 10 - 100 would be required to reveal either star. It is likely that a significant fraction of their wide fields will be magnified by such factors (e.g., Oguri & Marshall, 2010; Rydberg et al., 2018). Higher magnifications may be possible in future surveys of individual cluster lenses by JWST but at the cost of smaller lensing volumes (Whalen et al., 2013a; Windhorst et al., 2018) that enclose fewer objects.

How could blue rapidly accreting SMSs be distinguished from hot blue dark stars of similar mass and redshift? Perhaps the greatest distinction between the two objects is the dense accretion shroud of the SMS, which imprints prominent continuum absoption features on the its spectra redward of Ly$\alpha$ in the rest frame that are absent from those of blue dark stars (compare Figures 2c and 2d to Figure 6 in Freese et al., 2016). In principle, these spectral features could be used to distinguish blue SMSs from hot dark stars of similar mass. Blue SMS spectra also exhibit very prominent Ly$\alpha$ lines due to pumping of the accretion envelope by high-energy UV photons from the star. Dark star spectra lack this feature because they do not have dense envelopes but, as discussed earlier, it is not clear if it could be detected today because of resonant scattering of Ly$\alpha$ by the neutral IGM at $z>$ 6.

The photospheric temperatures of supermassive Pop III stars in atomically cooling haloes (whether they are red, blue, or both in the rest frame) remain an unsolved problem. Although numerical simulations to date broadly agree on the evolution and final fates of rapidly-accreting Pop III SMSs there are key differences between them that remain poorly understood, such as the growth of the convective core mass, the final masses of the most rapidly accreting objects, and the inflation of the photosphere. But there is a general consensus that these discrepancies likely arise from differences in how the models flag the onset of convection (i.e., the Schwarzschild or Ledoux criteria), their ability to follow dynamical instabilities (e.g., KEPLER) or not (e.g., GENEVA and Yorke & Bodenheimer, 2008), and their numerical resolution and boundary conditions at the surface.

Using a code derived from Yorke & Bodenheimer (2008), Hosokawa et al. (2013) found that H*-* opacity in the outermost envelopes of SMSs can greatly expand their photospheres and limit their temperatures to $\sim 0.5$–$1\times 10^{4}$K until becoming blue at masses $\gtrsim 10^{5}M_{\odot}$. Woods et al. (2017), however, find SMSs to be compact and blue from early times in the KEPLER stellar evolution code, without enough H*-* in their atmospheres to expand and cool them. A third study by Haemmerlé et al. (2018a) with the GENEVA code found that rapidly-accreting SMSs are persistently red throughout their lifetimes, although more slowly-accreting ones may be blue (see discussion below). Efforts to benchmark these studies and converge on a solution continue (see, e.g., the recent review by Woods et al., 2018), but the final answer may only come from observations like those proposed here.

Our simulations neglect the effect of radiation pressure due to flux from the star on the flows that create it. Including these effects in cosmological simulations is challenging because they must resolve photospheres, the inner regions of accretion disks, and how the two are connected on very small scales that prevent codes from evolving them for long times. Smith et al. (2017) post processed simulations of highly-resolved atomically cooling haloes with Ly$\alpha$ photon transport and found it could alter flows onto the star. Radiation hydrodynamical calculations by Luo et al. (2018) and Ardaneh et al. (2018) without resonant Ly$\alpha$ scattering found that radiation from the protostar did not significantly alter local flows at early times but did suppress fragmentation, thus promoting the rapid growth of a single, supermassive object. There is somewhat more of a possibility that ionising flux from blue SMSs could blow out gas from the disk and partially expose it to the IGM, but this will simply result in AB magnitudes closer to those of the bare star shown in Figures 3 and 4.

Pulsations in blue or red SMSs (e.g., Osaki, 1966) could improve their prospects for detection by temporarily boosting their fluxes above the detection limits of surveys. This phenomenon is not captured by the stellar evolution codes used here because their implicit solvers and large time steps do not resolve such oscillations, but it might cause the star to periodically brighten and dim by an order of magnitude on timescales of a few weeks to months in the rest frame. These oscillations might also facilitate their detection because their regularity would differentiate them from dusty, red high-$z$ quasars or low-$z$ interlopers such as cool Milky Way stars. Periodic dimming and brightening in principle could flag high-$z$ SMSs in transient surveys proposed for JWST such as FLARE (Wang et al., 2017).

The original studies on the collapse of pristine, atomically cooled haloes and the subsequent formation of SMSs assumed very high LW backgrounds that completely suppressed H2 formation in their cores, so collapse was nearly isothermal at temperatures of $\sim$ 8000 K and flow rates of 0.1 - 1 M*⊙* yr*-1*. This is why we adopted them as the two fiducial rates in our study (they are also typical of the collapse of atomically cooled haloes due to supersonic baryon streaming motions; Hirano et al., 2017). But in lower LW backgrounds some H2 can form in the core of the halo and enhance cooling there, leading to lower infall rates of a few 10*-3* M*⊙* yr*-1* (e.g., Latif et al., 2015; Regan & Downes, 2018). Such rates result in much less massive stars, perhaps 103 - 104 M*⊙* rather than 105 M*⊙*. It is not clear at this point which of these two populations of SMS was more prevalent in the early universe because average LW background strengths are not well understood and supersonic streaming motions are thought to have produced about as many SMSs as LW backgrounds. Furthermore, it is not clear if these stars evolved along hot blue tracks or cool red ones, although there are indications that some would be blue (Haemmerlé et al., 2018a). The prospects for detection of this second, less massive population of SMSs are unclear because it is not yet known if they were red or blue and they evolved in accretion envelopes with lower densities than those considered here. But they may be more difficult to find because of their lower fluxes. They will be studied in future work.

DCBH birth after the collapse of the SMS is the next stage of primordial quasar evolution, and a number of studies have examined prospects for their detection in future NIR surveys. These objects are also deeply embedded in dense, hot flows and techniques similar to those used here are required to model their spectra. One-dimensional radiation hydrodynamics simulations of DCBH emission post processed with Cloudy have shown that they could be detected by JWST out to $z\sim$ 20 (Pacucci et al., 2015; Natarajan et al., 2017). We will next post process radiation hydrodynamical simulations of the H II regions of DCBHs from $z=$ 10 - 20 with Cloudy to find out to what redshifts they could be found by Euclid, WFIRST and JWST.

Acknowledgements

The authors thank the anonymous referee, whose comments improved the quality of this paper. DJW was supported by the STFC New Applicant Grant ST/P000509/1. EZ acknowledges funding from the Swedish National Space Board. TH was supported by JSPS KAKENHI Grant Number 17F17320 and SCOG was funded by the European Research Council via the ERC Advanced Grant STARLIGHT: Formation of the First Stars (project number 339177) and from the DFG via SFB 881 “The Milky Way System" (sub-projects B1, B2, B8). AH was supported by TDLI though a grant from the Science and Technology Commission of Shanghai Municipality (Grants No.16DZ2260200) and National Natural Science Foundation of China (Grants No.11655002). SCOG was funded by the European Research Council via the ERC Advanced Grant STARLIGHT: Formation of the First Stars (project number 339177) and from the DFG via SFB 881 “The Milky Way System" (sub-projects B1, B2, B8). The simulations were performed on the Sciama High Performance Compute (HPC) cluster which is supported by the Institute of Cosmology and Gravitation (ICG), SEPnet and the University of Portsmouth.