EDGES and JWST with 21cm global signal emulator

The 21cm line is observed as emission/absorption against the radio background such as the CMB, and then the 21cm line brightness temperature is given as (e.g. Furlanetto et al., 2006a)

where $T_{\rm S}$ is the spin temperature, $x_{\rm HI}$ is the neutral fraction of hydrogen gases, $\delta$ is matter over-density and $T_{\gamma}=2.725/(1+z)$ is the temperature of CMB. Note that we do not consider the enhanced radio backgrounds (e.g. Feng & Holder, 2018; Fialkov & Barkana, 2019; Reis et al., 2020; Ewall-Wice et al., 2020) proposed to reproduce the very strong absorption signal reported in Bowman et al. (2018).

We used a modified version of 21cmFASTv2²²2There are several recent updates on the 21cmFAST such as the inhomogeneous Lyman-Werner and relative-velocity feedback (Muñoz et al., 2022). (Mesinger et al., 2011; Park et al., 2019), which is similar to the methodology used in Yoshiura et al. (2018); Minoda et al. (2023), for evaluating the 21cm line brightness temperature. Specifically, we added the running $\alpha_{\rm s}$ and the running of running $\beta_{\rm s}$ to the primordial power spectrum. The detailed method used in the 21cmFAST for solving the ionization and evolution of spin temperature is well described in Mesinger et al. (2011); Park et al. (2019). We here describe the method briefly. The initial matter power spectrum will be used for creating an initial high-resolution matter density map and estimation of the halo mass function³³3 The choice of halo mass function, the stellar population synthesis, and cosmology can affect the 21cm global signal with a few - tens of mK (Mirocha et al., 2021). . The three-dimensional matter density map of each redshift bin is evaluated using 2LPT from the initial matter distribution with high resolution. The ionization is solved by comparing the number of neutral baryons and recombination with ionizing photons (Sobacchi & Mesinger, 2014). The heating of IGM is caused by the X-ray emission⁴⁴4 Note that additional heating processes such as the Lyman-$\alpha$ heating produced by the Lyman-$\alpha$ photons can change the shape of global signal (Villanueva-Domingo et al., 2020; Ghara & Mellema, 2020; Mittal & Kulkarni, 2021). , calculated as angle-averaged specific intensity, with the optical depth being taken into account. The Lyman-$\alpha$ background is originated from the X-ray excitation of the neutral hydrogen atom and direct stellar emission of the Lyman series. The 21cm signal is evaluated using the maps of neutral fraction, spin temperature, and matter density. The global signal (i.e. sky averaged) 21cm signal is calculated using the resultant 21cm brightness temperature maps.

The emission from luminous objects is controlled by several astrophysical parameters, which we list following Park et al. (2019); the fraction of gas converted into stars in a halo of mass $10^{10}M_{\odot}$, $f_{\rm star}$, the power-law index of $f_{\rm star}$ as a function of halo mass, $\alpha_{\rm star}$, a fraction of time scale of the star formation against Hubble time, $t_{\rm star}$, the escape fraction of ionizing photons of a halo of mass $10^{10}M_{\odot}$, $f_{\rm esc}$, the power-law index of $f_{\rm esc}$ as a function of halo mass, $\alpha_{\rm esc}$, turn over mass, $M_{\rm turn}$, which is the lowest halo mass to host the efficient star formation, X-ray luminosity per star formation, $L_{\rm X}$, and minimum X-ray energy, $E_{0}$.

In addition to astrophysical parameters, the small-scale primordial power spectrum can also play an important role in the evolution of the 21cm global signal. The primordial power spectrum can be parametrized as

where $n_{s}$ is the spectral index, $\alpha_{s}$ and $\beta_{s}$ are the running and the running of the running parameters. $k_{\rm ref}=0.05{\rm Mpc^{-1}}$ is the reference scale. In Yoshiura et al. (2018, 2020a), the effects of $\alpha_{\rm s}$ and $\beta_{\rm s}$ on the 21cm global signal are investigated. The enhancement of the matter power spectrum owing to the positively large value of $\alpha_{\rm s}$ and $\beta_{\rm s}$ makes the structure formation faster. As a result, the redshift of the 21cm absorption signal is shifted to earlier. Thus the primordial power spectrum is also important in predicting the 21cm global signal and we vary $n_{s}$, $\alpha_{s}$, and $\beta_{s}$ as well as the astrophysical parameters in the following analysis.

We run the 21cmFAST for roughly 49000 different parameter sets in the PC cluster of CfCA and store the global 21cm line signal. The simulation was performed in a cosmological volume of $128{\rm Mpc^{3}}$, $192^{3}$ cells for initial condition and $64^{3}$ cells for ionized and temperature maps⁵⁵5 We tested various sets of resolutions and find that the resolution does not significantly affect the evolution of the 21cm global signal. . We divide the redshift into 53 bins over the range of $6<z<30$. Parameters are randomly selected from a range motivated by Park et al. (2019). The ranges of parameters of the primordial power spectrum ($n_{\rm s}$, $\alpha_{\rm s}$, and $\beta_{\rm s}$) are consistent with the 5 $\sigma$ confidence level of the constraints from Planck Collaboration et al. (2020). The ranges of parameters are listed in Table 1. The set of parameters and the resultant 21cm global signal is used for training our emulator as described below. We will use roughly 44000 data as training data and the rest as test data.

As the emulator of a 21cm global signal, we employ an ANN-based architecture which is similar to some previous works (Cohen et al., 2020; Bevins et al., 2021b; Bye et al., 2022). We perform the optimization of ANN models using 12 different architectures⁶⁶6 The 12 models are the combinations of the number of hidden layers (2, 3, and 4) and the number of neurons at the hidden layers (4,8,16, and 32). . The hidden layers have a Tanh activation function except for the last layer which has a linear activation. The input layer of the network for the global signal has 11 inputs (8 astrophysical parameters, 3 primordial power spectrum parameters), and 1 redshift index. The redshift indices (from 1 to 53) correspond to the output redshift bins of 21cmFAST. The network yields one output corresponding to the brightness temperature at the redshift index. The loss function used to optimize the network is Mean Squared Error. We use the batch size of 53, ADAM optimizer, and learning rate of 0.002. We train the network for 100 epochs using the all training data set. As quantitative evaluators of the accuracy of our network, we employ root mean squared error (RMSE) as in Bevins et al. (2021b); Cohen et al. (2020). The RMSE for a test data set is given as

where $N$ is number of redshift bins, $z_{j}$ is redshift at $j$-th bin, $\delta T_{\rm sim}$ is simulated 21cm brightness temperature from 21cmFAST and $\delta T_{\rm pred}$ is 21cm signal predicted using the emulator. As a result, the architecture consists of 1 input layer, 4 hidden layers with 32 neurons for each layer, and 1 output layer, achieving the mean RMSE of 2.45 mK which is the smallest among the other architectures. The figure 1 shows examples of the 21cm global signal predicted by our emulator. The accuracy might be improved further by optimization of training parameters and the architecture. However, this RMSE is more or less consistent with the previous work and lower than the expected thermal noise error of the EDGES data at a frequency bin. Thus, we consider that the current emulator is accurate enough for this work.

Note that the redshift bins of our training data sets do not correspond to the frequency bins of the EDGES data one by one. Thus, in the Bayesian analysis, we predict the 21cm brightness temperature by linearly interpolating the values to evaluate the 21cm signal at the frequency bins of the EDGES data.

Using the parameters used in the 21cmFAST, we can evaluate the high-$z$ UVLF following Park et al. (2019). The star formation rate (SFR) is described as

where $H^{-1}_{z}$ is Hubble time at the redshift $z$ and the stellar mass $M_{\rm*}$ is given as

where $M_{\rm h}$ is halo mass and the $f_{\rm star,m}$ is less than 1 and given as $f_{\rm star,m}$ = $f_{\rm star}M^{\alpha_{\rm star}}_{\rm h,10}$. $M_{\rm h,10}$ is halo mass normalized by $10^{10}M_{\odot}$. The SFR is converted to rest-frame UV luminosity by the conversion factor $K_{\rm UV}$ = $1.15\times 10^{-28}M_{\odot}~{}{\rm yr^{-1}/ergs~{}s^{-1}Hz^{-1}}$ (assuming the Salpeter IMF, Sun & Furlanetto, 2016). The UV magnitude $M_{\rm UV}$ is derived from the luminosity by following AB magnitude relation (Oke & Gunn, 1983). Finally, the UVLF is given as

where $f_{\rm duty}=\exp(-M_{\rm turn}/M_{\rm h})$ is a suppression factor for star formation in small halos. The halo mass function ${dn}/{dM_{\rm h}}$ is calculated following (Murray et al., 2013) with the parametrization for the primordial power spectrum⁷⁷7We here assume the amplitude of the primordial power spectrum at the reference scale is $2.207\times 10^{-9}$ (Planck Collaboration et al., 2016). (2) as in e.g. Yoshiura et al. (2020b), using the top-hat window function.

Previous works (Park et al., 2019; Yoshiura et al., 2020b) have shown that the UVLF at $z\leq 10$ is a powerful quantity to constrain the astrophysical parameters and primordial power spectrum. However, when we perform the Bayesian analysis using the EDGES data corresponding to $z>13$, we do not include the UVLF in the likelihood. We calculate the UVLF at $z=12$ and $z=16$ only to discuss the implication of EDGES to the recent JWST results.

Our emulator is constructed with 11 parameters. We demonstrate the emulator with the EDGES low-band data which covers the frequency range of the 21cm line before the cosmic reionization. Thus, we fix the escape fraction as the $\log_{10}f_{\rm esc}=-1.5$ and $\alpha_{\rm esc}=-0.25$ since these parameters are mainly relevant to the reionization. We also consider a model with fixing the parameters for the primordial power spectrum as $n_{\rm s}=0.9665$, $\alpha_{\rm s}=0.0$, and $\beta_{\rm s}=0.0$. It would be worth mentioning that the model with $n_{\rm s}=0.9625$, $\alpha_{\rm s}=0.002$, and $\beta_{\rm s}=0.01$, which corresponds to the mean values of Planck constraint (Planck Collaboration et al., 2020), has the absorption peak at higher redshifts than that of the model with $\alpha_{\rm s}=\beta_{\rm s}=0$ as shown as the thin solid line in figure 2. However, the influence is minor compared to the effects of astrophysical parameters. We also examine the model without a 21cm signal.

We assess the foreground models with orders of polynomials from 5th to 9th. The prior range of the parameters used in the fitting is listed in Table 1 and 2. We perform the analysis for multiple combinations of these models as listed in Table 3.

Regarding the systematics, we found that the model without the sinusoidal systematics, described as Eq. (11), has significantly larger RMS and lower $\ln\mathcal{Z}$. This might suggest the existence of sinusoidal systematics. Therefore we decided to use sinusoidal systematics in all models.

We perform the Bayesian analysis using all frequency bands (51-99MHz) or only lower bands (51-87 MHz) and use 2 different noise models. For clear discussion, we define 4 cases for the combination of the frequency band and noise model, namely FullBW(full bandwidth), HcutBW (high-frequency cut bandwidth), TN (theoretical noise, $A_{\rm n}=1$, $T_{\rm wn}=0$ K), and EN (enhanced noise, $A_{\rm n}=2.25$, $T_{\rm wn}=0.015$ K). In addition, depending on the model framework as to how we fix the parameters, we separate the models into 4 labelled as M1, M2, M3 and M4 (see Table 3). We also define 11 IDs for different combinations of 21cm signal and foregrounds (see Table 3). We refer, for example, to the model of M1-FG8 and FullBW-TN as FullBW-TN-M1-FG8.

The Bayesian evidence values of all models are summarized in Table 3. We compare the evidence of 44 models¹²¹²12 The Bayesian evidence indicates how well the model describes the data. The result of 51¡$\nu$¡99 MHz should not be directly compared with the result of 51¡$\nu$¡87 MHz. . The evidence values become significantly high if we assume an enhanced noise model (i.e. $A_{\rm n}$=2.25, $T_{\rm wn}$=0.015). This indicates the presence of additional noise sources such as calibration error and polarized foregrounds as discussed in Sims & Pober (2020) and/or the presence of unknown systematic error which is not properly modeled with the sinusoidal function. It should be mentioned that the values assumed for the enhanced noise ($A_{\rm n}$=2.25, $T_{\rm wn}$=0.015) are not physically motivated ones. Thus, additional systematic spectrum features may improve the evidence values without the enhancement of noise. We also find that the models with 5th-order polynomials have low evidence values except for HCutBW-EN. Thus the polynomial function of higher than 6th-orders is required to remove the foregrounds.

In many cases except for FullBW-TN, the evidence is high even without the 21cm signal. Thus, we could not confirm the presence of a 21cm absorption line. In FullBW-TN, we find the statistically significant difference of $\ln Z$ between the FullBW-TN-M4-FG9 (without 21cm signal) and FullBW-TN-M3-FG9 (with 21cm signal). However, the evidence is much lower than that of FullBW-EN. This might caution that the analysis of a 21cm global signal can lead to noise-biased results. However, it should be emphasized that our results do not draw any conclusion regarding the strong absorption reported in Bowman et al. (2018) since our emulator cannot produce the flattened Gaussian shape and absorption of 500 mK reported by EDGES.

We discuss the result of FullBW-EN-M1-FG7 in detail as an example model. Figure 4 shows the posterior distribution of astrophysical parameters. The shape of the posterior distribution of all astrophysical parameters is almost flat over the prior range. Figure 5 shows the 21cm global signal calculated from the posterior samples¹³¹³13 We use fgivenx (Handley, 2018). . The constraint is consistent with zero within 1 $\sigma$ error. The blue region shows the posterior of the 21cm global signal from random samples within priors¹⁴¹⁴14 Specifically, we create 20000 samples of 21cm global signal from random parameters within the prior range. The posterior density distribution, shown as the blue region, is evaluated as 3 $\sigma$ region of the random samples using fgivenx. . The figure indicates that we cannot argue the detection of the 21cm line from this model.

The strong upper limits on the 21cm power spectrum observation measured by HERA (HERA Collaboration et al., 2023) have put the constraints on the spin temperature at $z=10.4$ and $z=7.9$ using various 21cm models including 21cmFAST. Thus, joint analysis with the power spectrum result would give a tighter lower limit on the 21cm global signal at $z<10$. In Bevins et al. (2023), they have performed a joint Bayesian analysis using the SARAS3 and HERA observation data. Our constraints on the 21cm signal are consistent with their result while they have not used the 21cmFAST to obtain the 21cm signal and imposed a wider prior range in the likelihood analysis.

The model of star formation rate used in the 21cmFAST can be applied to evaluate the UVLF following Park et al. (2019). Thus, we can also obtain the posterior distribution of the UVLF. Figure 6 shows the constraint on the UVLF at $z=12$ and $z=16$. The results are well consistent with the recent JWST results (Harikane et al., 2023). The central value of UVLF at $z=16$ slightly offsets the centre of our confidential region.

We first examine the effect of the order of polynomials by comparing M1-FG5 to M1-FG9 of FullBW-EN. The constrained global signal is similar to each other and consistent with zero within 1 $\sigma$ error except the FullBW-EN-M1-FG5. The FullBW-EN-M1-FG5 with the lower polynomial order of the foreground model suggests an absorption at $z\approx 17$ with a peak amplitude of roughly 150 mK but has lower $\ln\mathcal{Z}$ (Figure 7). This result caution that an improper assumption on the order of the polynomial can potentially lead to a biased detection.

Comparing FullBW-EN-M1-FG7 and HCutBW-EN-M1-FG7, we find the constraints on the global signal are similar and consistent with zero at 1 $\sigma$. HCutBW-EN-M1-FG7 (i.e. removing data at $\nu>87\rm MHz$) also shows weak constraints on astrophysical parameters (Figure 4).

Now, we discuss the case of the other noise model, which is theoretically motivated thermal noise. The FullBW-TN (HCutBW-TN) results are consistent with zero within 2 (1) $\sigma$ although we can find complicated structures in the posterior global signal for some cases (e.g. absorption features in Figure 8.). The structure should be caused by clusters of posterior distributions of foreground parameters. For example, the posterior distributions of $d_{0}$ and $d_{1}$ show clusters as shown in Figure 9. Even with a small difference in foreground parameters, the impact can be larger than the 21cm signal because the foreground is significantly bright. The models with theoretical noise model tend to show such clustering posterior distribution of the foreground parameters. This clustering posterior distribution of the foreground parameters might be interpreted as misfitting due to the underestimation of noise. There are other possible noise sources, such as the calibration error or the polarized foreground as described in section 3. We checked that there are two or three peaks in the posterior distributions of astrophysical parameters. Each possible absorption line (Figure 8) might correspond to each cluster (Figure 9).

The model with higher $\alpha_{\rm s}$ and higher $\beta_{\rm s}$ can enhance the halo mass function. Therefore, the timing of Lyman-$\alpha$ coupling and X-ray heating is shifted to higher redshift, and hence the position of the 21cm absorption trough moves to a higher redshift as shown in Figure 2. The variation of the primordial power spectrum introduces more degrees of freedom to the shape and position of the global signal profile. Comparing the FullBW-EN-M1-FG9 (fixing $n_{\rm s},\alpha_{\rm s},\beta_{\rm s}$) and FullBW-EN-M2-FG9 (varying $n_{\rm s},\alpha_{\rm s},\beta_{\rm s}$) can give insight into the effect of the primordial power spectrum on the astrophysical parameter constraints (Figure 10). The constraints are weak for both cases while the posterior distributions are not perfectly consistent with each other. For the constraints on parameters of the primordial power spectrum, we also find an almost flat posterior distribution except FullBW-TN-M2-FG9 (Figure 11). The figure indicates that constraints on the primordial power spectrum from 21cm global signal depend on the assumptions on foreground and noise in the analysis. For example, the case of FullBW-TN-M2-FG9, we find an upper limits on the $\beta_{\rm s}$, in which the constraint might be affected biased by unknown systematics and underestimated error. In any case, our analysis indicates that the 21cm global signal has the potential to constrain the running parameters if the shape of the 21cm global signal is well constrained.

As the EDGES data corresponds to the redshift larger than 13, and the ionization would not happen at such high redshift, therefore we fixed parameters related to the escape fraction of ionizing photons in models M1 and M2. On the other hand, in the M3 model, we treat the $f_{\rm esc}$ and $a_{\rm esc}$ as free parameters. In the FullBW-EN-M3-FG9, we find that the posterior distributions of parameters are almost flat and the constraints on the global signal are consistent with zero.

The models of high log ${\cal Z}$ provide the posterior luminosity function consistent with the JWST results as shown in figure. 6. However, the constraints on the luminosity function is weak. On the other hand, in the FullBW-TN-M1-FG9, for example, the constraints on the UVLF become tight and low values of the UVLF are disfavored as shown in Figure 12. Furthermore, in the FullBW-TN-M1-FG9, the presence of a certain absorption line at $z>13$ is preferred in the posterior distribution of the global signal. To create the absorption feature, an effective star formation rate is required to produce a large amount of Lyman-$\alpha$ photons enough to couple the spin temperature with the cold gas temperature. Therefore the recent discovery of high UVLF at high-$z$ might be consistent with the 21cm absorption line at $z\approx 15$. It should be mentioned that the favoured parameters require a high value of $\log_{10}f_{\rm star}>-0.5$ and a low value of $t_{\rm star}<0.3$. We note that Haro et al. (2023) recently derived the follow-up spectroscopic result of JWST and find one of the galaxies candidates at $z=16$ was the galaxy at $z=4.9$ while the galaxy suspected to be at $z>10$ was confirmed as the galaxy at $z=11.4$.

We mention that the required star formation rate to explain the EDGES absorption line has been discussed in Madau (2018). Very recently, Hassan et al. (2023) pointed out the UV luminosity density required to explain the EDGES absorption signal is consistent with the recent JWST observation ($z<15$) and the extrapolation of the galaxy UV observed by HST. They also show that the UVLF at $z=16$ necessary to explain the EDGES absorption is slighly higher than the UVLF extrapolated from UVLF at $z=12$.

In Bera et al. (2022), they have shown such a model is consistent with current observations (e.g. neutral fraction) and pointed out the JWST can probe it. Thus, it should be interesting to perform MCMC analysis simultaneously using the UVLF, neutral fraction, 21cm global signal, and so on.