Interpreting LOFAR 21-cm signal upper limits at $z\approx 9.1$ in the context of high-$z$ galaxy and reionisation observations

It is assumed that the typical stellar mass of a galaxy, $M_{\ast}$, can be directly related to its host halo mass, $M_{\rm h}$ (e.g. Kuhlen & Faucher-Giguère, 2012; Dayal et al., 2014; Behroozi & Silk, 2015; Mitra et al., 2015; Mutch et al., 2016; Sun & Furlanetto, 2016; Yue et al., 2016) through a power-law relation⁴⁴4At $z\gtrsim 5$ this power-law dependence between stellar mass and halo mass is consistent with the mean relation from semi-analytic model predictions (e.g Mutch et al., 2016; Yung et al., 2019) and semi-empirical fits to observations (e.g. Harikane et al., 2016; Tacchella et al., 2018; Behroozi et al., 2019). normalised to a dark matter halo of mass $10^{10}$ $M_{\odot}$:

where $f_{\ast}$ is the fraction of galactic gas in stars,

with $f_{\ast,10}$ being the normalisation of the relation and $\alpha_{\ast}$ its power-law index.

The star-formation rate (SFR) can then be estimated by dividing the stellar mass by a characteristic time-scale,

where $H^{-1}(z)$ is the Hubble time and $t_{\ast}$ is a free parameter allowed to vary between zero and unity.

Equivalently, the UV ionising escape fraction, $f_{\rm esc}$, is also parameterised to vary with halo mass,

with $f_{\rm esc,10}$ again normalised at a halo of mass $10^{10}$ $M_{\odot}$ with the power-law index, $\alpha_{\rm esc}$.

Finally, to mimic the inability of small mass halos to host active, star-forming galaxies (i.e. because of inefficient cooling and/or feedback processes), a duty-cycle is included to suppress their contribution:

In effect, some fraction, $(1-f_{\rm duty})$, of dark matter haloes with mass $M_{\rm h}$ are unable to host star-forming galaxies, with the characteristic scale for this suppression being set by $M_{\rm turn}$ (e.g. Shapiro et al., 1994; Giroux et al., 1994; Hui & Gnedin, 1997; Barkana & Loeb, 2001; Springel & Hernquist, 2003; Mesinger & Dijkstra, 2008; Okamoto et al., 2008; Sobacchi & Mesinger, 2013a, b).

In summary, this results in six free parameters describing the UV galaxy properties in our model, $f_{\ast,10}$, $f_{\rm esc,10}$, $\alpha_{\ast}$, $\alpha_{\rm esc}$, $M_{\rm turn}$ and $t_{\ast}$.

Prior to reionisation, it is thought that the IGM undergoes heating in the early Universe due to X-rays. The likely, dominant source of these X-rays are stellar remnants within the first galaxies. In 21CMFAST, X-ray heating is included by calculating the cell-by-cell angle-averaged specific X-ray intensity, $J(\boldsymbol{x},E,z)$, (in erg s^-1 keV^-1 cm^-2 sr^-1), by integrating the co-moving X-ray specific emissivity, $\epsilon_{\rm X}(\boldsymbol{x},E_{e},z^{\prime})$ back along the light-cone:

where $e^{-\tau}$ accounts for attenuation by the IGM. The co-moving specific emissivity, evaluated in the emitted frame, $E_{\rm e}=E(1+z^{\prime})/(1+z)$, is,

where $\bar{\delta}_{\rm nl}$ is the mean, non-linear density in a shell around $(\boldsymbol{x},z)$ and the quantity in square brackets is the SFR density along the light-cone.

The normalisation, $L_{\rm X}/{\rm SFR}$ (erg s^-1 keV^-1 $M^{-1}_{\odot}$ yr), is the specific X-ray luminosity per unit star formation escaping the host galaxies. This is assumed to be a power-law with respect to photon energy, $L_{\rm X}\propto E^{-\alpha_{\rm X}}$, which is attenuated below a threshold energy, $E_{0}$, where photons are absorbed inside the host galaxy. The specific luminosity is then normalised to an integrated soft-band ($<2$ keV) luminosity per SFR (in erg s^-1 $M^{-1}_{\odot}$ yr), which is taken to be a free parameter:

This limit of $2\,{\rm keV}$ approximates to an X-ray mean-free path of roughly the Hubble length at high redshifts, implying harder photons do not contribute to heating the IGM (e.g. McQuinn, 2012).

In summary, we have three free parameters describing the X-ray properties, $L_{{\rm X}<2\,{\rm keV}}/{\rm SFR}$ , $E_{0}$ and $\alpha_{\rm X}$.

The evolved density and velocity fields in 21CMFAST are obtained following second-order Lagrange perturbation theory (e.g Scoccimarro, 1998). Reionisation is then determined from the evolved density field by tracking the balance between the cumulative number of ionising photons against the number of neutral hydrogen atoms plus cumulative recombinations in spheres of decreasing radii. Each cell is flagged as ionised when,

where $\bar{n}_{\rm rec}$ is the cumulative number of recombinations (e.g. Sobacchi & Mesinger, 2014) and $n_{\rm ion}$ is the cumulative number of IGM ionising photons per baryon inside a spherical region of size, $R$ and corresponding overdensity, $\delta_{R}$,

where $\bar{\rho}_{b}$ is the mean baryon density and $N_{\gamma/b}$ is the number of ionising photons per stellar baryon⁵⁵5We adopt a value of 5000, corresponding to a Salpeter initial mass function (Salpeter, 1955); however note that this is highly degenerate with $f_{\ast}$. The final term of Equation 9, $(1-\bar{x}_{e})$, corresponds to ionisations by X-rays, which are expected to contribute at the $\sim 10$ per cent level (e.g. Ricotti & Ostriker, 2004; Mesinger et al., 2013; Madau & Fragos, 2017; Ross et al., 2017; Eide et al., 2018).

The thermal state of the neutral IGM (and its partial ionisations) are tracked in each cell by following adiabatic heating/cooling, Compton heating/cooling, heating through partial ionisations, as well as the heating/ionisations from X-rays (discussed in the previous section). Combining all these, the IGM spin temperature, $T_{\rm S}$, is then computed as a weighted mean between the gas and CMB temperatures, depending on the density and local Lyman-$\alpha$ (Ly$\alpha$) intensity impinging on each cell (Wouthuysen, 1952; Field, 1958).

This Ly$\alpha$ background (see e.g. Mesinger et al. 2011 for further details) is estimated by summing the contribution from: (i) excitations of neutral hydrogen by X-rays ($J_{\alpha,{\rm X}}$) and (ii) direct stellar emission of photons between Ly$\alpha$ and the Lyman limit ($J_{\alpha,\ast}$). For (i), this is set by the X-ray heating rate assuming that energy injection is balanced by photons redshifting out of Ly$\alpha$ resonance (Pritchard & Furlanetto, 2007). For (ii), any Lyman-$n$ resonance photons absorbed by the neutral IGM cascade with a recycling fraction passing through Ly$\alpha$ producing a background that is the sum over all Lyman resonances (e.g. Barkana & Loeb, 2005). Note, the soft UV spectra of the stellar emission component of the first sources is currently held fixed and we do not include other possible sources of soft UV, such as quasars (see e.g. Qin et al., 2017; Ricci et al., 2017; Mitra et al., 2018).

As 21CMFAST employes an excursion-set formalism for determining which cells are ionised, it is impacted by the resultant non-conservation of ionising photons. This arises within cells that exceed the ionisation criteria (Eq. 9), where the excess ionising photons that are left over after ionising the neutral hydrogen are not propagated onward. This behaviour acts as an effective bias on the ionising escape fraction, $f_{\rm esc}$. Roughly speaking, this corresponds to a loss of $\sim 10-20$ per cent of the ionising photons (e.g. McQuinn et al., 2005; Zahn et al., 2007; Paranjape & Choudhury, 2014; Paranjape et al., 2016; Hassan et al., 2017; Choudhury & Paranjape, 2018; Hutter, 2018; Molaro et al., 2019).

Both Choudhury & Paranjape (2018) and Molaro et al. (2019) have introduced new, explicitly photon conserving algorithms for semi-numerical simulations. While being orders of magnitude faster than full radiative-transfer simulations, they are still considerably slower than conventional semi-numerical simulations. Under Bayesian parameter estimation approaches such as 21CMMC, these schemes are still intractable when forward modelling the 21-cm signal in the high-dimensional parameter spaces required to characterise the ionising, soft UV, and X-ray properties of the first galaxies.

Instead, Park et al., in prep, introduce an approximate correction to the ionising photon non-conservation issue, correcting for the effective bias on $f_{\rm esc}$ by analytically solving for the correct evolution of the ionisation fraction for a given source model, under the assumption of no correlations between the sources and sinks (which should only impact the final $\sim 10$ per cent of the EoR; Sobacchi & Mesinger 2014).

This analytic solution is compared against a calibration curve from 21CMFAST considering only ionisations (i.e. no recombinations or spin temperature evolution). If photons were conserved, these two curves would be identical. However, in practice the 21CMFAST history is delayed due to the loss of ionising photons. In order to recalibrate 21CMFAST, we decrease the redshift used for determining the ionisation field⁶⁶6This effectively amounts to boosting the number of ionisations to compensate for the ionising photon non-conservation. However, we do not know $a$-priori how much to boost ionisations by, thus we modify the redshift to allow our calculation to be performed on-the-fly.. Note, this modified redshift is only used for modifying the ionisation field, all other quantities are calculated at the original redshift. We find that for the bulk of the EoR this corresponds roughly to a shift in redshift of $\Delta z\sim 0.3\pm 0.1$.

Combining all the cosmological fields discussed in the previous section, we compute the observable cosmic 21-cm signal as a brightness temperature contrast relative to the CMB temperature, $T_{\rm CMB}$ (e.g. Furlanetto et al., 2006):

where $\tau_{\nu_{0}}$ is the optical depth of the 21-cm line,

Here, $x_{\mathrm{H{\scriptscriptstyle I}}{}}$ is the neutral hydrogen fraction, $\delta_{\rm nl}\equiv\rho/\bar{\rho}-1$ is the gas over-density, $H(z)$ is the Hubble parameter, ${\rm d}v_{\rm r}/{\rm d}r$ is the gradient of the line-of-sight component of the velocity and $T_{\rm S}$ is the gas spin temperature. All quantities are evaluated at redshift $z=\nu_{0}/\nu-1$, where $\nu_{0}$ is the 21-cm frequency and we drop the spatial dependence for brevity. Additionally, we include the impact of redshift space distortions along the line-of-sight as outlined in Mao et al. (2012); Jensen et al. (2013); Greig & Mesinger (2018).

In the middle right panel of Figure 1 we present a census of all existing constraints on the IGM neutral fraction against the full range of reionisation histories disfavoured by LOFAR (shaded region). We compare against the dark pixel statistics of high-$z$ quasars (QSOs; McGreer et al. 2015), Ly$\alpha$ fraction (Mesinger et al., 2015), the clustering of Ly$\alpha$ emitters (LAEs; Sobacchi & Mesinger 2015), the Ly$\alpha$ equivalent width distribution of Lyman-break galaxies (LBGs; Mason et al. 2018; Hoag et al. 2019; Mason et al. 2019), the neutral IGM damping wing imprint from high-$z$ QSOs (Greig et al., 2017; Davies et al., 2018; Greig & Mesinger, 2018) and the midpoint of reionisation ($z_{\rm Re}$) from Planck (Planck Collaboration et al., 2018).

As discussed in the previous section, the models disfavoured by the LOFAR upper limits are driven to produce a midpoint of reionisation at $z\approx 9.1$. The full range of IGM neutral fractions disfavoured by LOFAR at $z\approx 9.1$ correspond to $0.15\lesssim\bar{x}_{\mathrm{H{\scriptscriptstyle I}}{}}\lesssim 0.6$. Compared to the existing observational constraints, this places these disfavoured models at roughly $\gtrsim 2\sigma$. Unfortunately, owing to the still large amplitude of the LOFAR upper limits, and the correspondingly large uncertainties, we have verified that a joint analysis is entirely constrained by the existing limits. However, prospects for this will improve in the near future with the further processing of existing observational data and multiple frequency (redshift) limits.

Our constraints on $\bar{x}_{\mathrm{H{\scriptscriptstyle I}}{}}$ are broadly consistent with Ghara et al. (2020), with any differences being explained by different choices in adopted priors between the two works. The Ghara et al. (2020) limits on $\bar{x}_{\mathrm{H{\scriptscriptstyle I}}{}}$ appear from two regions of parameter space (c.f their Table 5): (i) ‘cold’ reionisation ($0.40<\bar{x}_{\mathrm{H{\scriptscriptstyle I}}{}}<0.55$) and (ii) patchy X-ray heating ($\bar{x}_{\mathrm{H{\scriptscriptstyle I}}{}}>0.92$). For (i), reionisation must be ongoing and the IGM must be cold. Our results are consistent with (i), given that Ghara et al. (2020) assume a fixed $M_{\rm turn}=10^{9}~{}M_{\odot}$ while also including a hard prior on the neutral fraction excluding values below 0.19, while we have a flat prior down to zero. For (ii), Ghara et al. (2020) require reionisation to be in its infancy and X-ray heating to be dominated by luminous, highly biased ($M_{h}\geq 10^{10}~{}M_{\odot}$) sources with very soft SEDs ($E_{0}=0.2$ keV). We do not recover their (ii) because our prior on $M_{\rm turn}$ has a hard upper limit of $M_{h}\leq 10^{10}~{}M_{\odot}$.

Unlike the reionisation history which is constrained to be within a narrow range around $\bar{x}_{\mathrm{H{\scriptscriptstyle I}}{}}\approx 0.5$, the UV LFs are not strongly constrained. Instead, they show extremely large scatter in the shapes and amplitudes, especially at the faint end. As such, we cannot simply provide a constrained range of plausible UV LFs as we did previously for the reionisation history (shaded region of middle right panel of Figure 1). Therefore, in the six panels in the top right of Figure 1 we plot the resultant UV luminosity functions at $z=6$, 7, 8, 10, 12 and 15 for 500 randomly sampled models drawn from the posterior which exceed the LOFAR upper limit. This is one of the advantages of using 21CMFAST in that the built-in parameterisation is capable of directly outputting UV LFs.

We compare these UV LFs against existing observations from unlensed UV LFs at $z=6-8$ (orange circles; Bouwens et al. 2015; Bouwens et al. 2017) and at $z=10$ (pink squares; (Oesch et al., 2018)). The vast majority of the UV LFs disfavoured by LOFAR are strongly disfavoured by the existing UV galaxy limits (up to two orders of magnitude larger in amplitude). Again, this is driven by models producing excessive amounts of ionising photons to have reionisation $\sim 50$ per cent complete by $z\approx 9.1$. While the UV LFs are independent of the number of ionising photons, to achieve cold reionisation at $z\approx 9.1$ we prefer to have both a high $f_{\rm esc}$ and high stellar-to-halo mass relation (see the $f_{\ast}-f_{\rm esc}$ panel in Figure 1). Interestingly, unlike the existing constraints on the IGM neutral fraction, there appear to be several models in excess of the LOFAR limits which appear to be capable of producing UV LFs consistent with the observed ones. However, the vast majority are strongly inconsistent, owing to the relatively small statistical uncertainties on the observed UV LFs.

Once again, this highlights the fundamental value of observing the cosmic 21-cm signal. Constraints on the amplitude of the 21-cm PS are capable of providing limits on the underlying UV galaxy LFs, beyond what is capable from existing space telescopes. Our galaxy model, which has an average star-formation time-scale evolving with redshift as the Hubble time (or analogously the halo dynamical time; Equation 3), and a redshift-independent stellar-to-halo mass relation is capable of reproducing current observations of high-$z$ ($z$ 6-10) UV LFs (e.g. figure 4 of Park et al. 2019). Moreover, these scaling relations are consistent with more sophisticated high-$z$ SAMs (e.g. Mutch et al., 2016; Yung et al., 2019), as well as empirical scaling relations (e.g. Tacchella et al., 2018). Therefore the fact that the UV LFs corresponding to the models inconsistent with the LOFAR observations are generally above the observed UV LFs is a significant result, highlighting that such models are already ruled out by current observations.

Next, in Figure 2 we consider the electron scattering optical depth, $\tau_{\rm e}$, using the latest constraints from Planck ($\tau_{\rm e}=0.054\pm 0.007$; Planck Collaboration et al. 2018). Here, the solid vertical line is the mean value from Planck, with shaded regions corresponding to being within $\pm 1\sigma$, $\pm 2\sigma$, $\pm 3\sigma$ of the mean observed value. The red curve is the histogram of $\tau_{\rm e}$ obtained by binning all models disfavoured by the LOFAR upper limits. Clearly, all models in excess of LOFAR are $\gtrsim 2\sigma$ from existing observational constraints. Note again, this arises owing to the assumptions of our astrophysical source parameterisation¹⁴¹⁴14We note that our parametrisation does allow for the population averaged $f_{\rm esc}$ to vary with redshift. This is because the halo mass function evolves with redshift, and we allow $f_{\rm esc}$ to scale with the halo mass through a power law relation whose index, $\alpha_{\rm esc}$, is a free parameter. However, our range of priors on $\alpha_{\rm esc}$ results in a relatively modest redshift evolution of the population averaged galaxy $f_{\rm esc}$.. For example, allowing a strong redshift dependence on $f_{\rm esc}$ could produce models consistent with Planck. However, we note that most hydrodynamical simulations do not find evidence for a strong redshift evolution for $f_{\rm esc}$ (e.g. Kimm & Cen 2014; Paardekooper et al. 2015; Xu et al. 2016; Ma et al. 2020, though see Lewis et al. 2020). The primary driver of this is the requirement to have reionisation at $\sim 50$ per cent by $z\sim 9.1$ for the 21-cm signal to exceed the current LOFAR limits, whereas the observed values from Planck prefer a redshift for the mid-point of $z_{\rm Re}=7.67\pm 0.5$ Planck Collaboration et al. 2018).

In Figure 3 we compare the mean UV background photoionisation rate from 500 randomly sampled models from the posterior in excess of the LOFAR upper limits to observational constrains extracted from the proximity zones of $z>6$ QSOs (e.g. Calverley et al., 2011; Wyithe & Bolton, 2011). Here, the vast majority of the models ruled out by LOFAR are at least two orders of magnitude larger in amplitude than those inferred from the observational constraints ($\gtrsim 3\sigma$). This is consistent with the picture of an excessive amount of ionising photons being required to ensure reionisation is $\sim 50$ per cent complete by $z\approx 9.1$. As a result, these models drastically overproduce the mean background photo-ionisation rate.