The cross-power spectrum between 21 cm emission and galaxies in hierarchical galaxy formation models

Abstract

INTRODUCTION

The prospect of measuring the 21-cm power spectrum from the epoch of reionization is a focus of modern theoretical cosmology (e.g. Morales & Wyithe 2010). A very successful technique has been to employ an N-body code to generate a distribution of haloes, and then apply radiative transfer methods in post-processing to model the generation of ionized structure on large scales using various models for the ionizing sources (e.g. Ciardi, Stoehr & White 2003; Sokasian et al. 2003; Iliev et al. 2007, 2008; Zahn et al. 2007; Trac & Cen 2007; Shin, Trac & Cen 2008; Trac, Cen & Loeb 2008). However, when constructing models to assign ionizing luminosities to dark matter haloes, most studies have used a constant mass-to-luminosity relation. On the other hand, the degree to which the important astrophysics governing formation and evolution of high-redshift galaxies will influence observations of the 21-cm power spectrum is not well known. To improve on the source modelling for calculation of the ionizing photon budget in reionization simulations, several studies (Benson et al. 2006; Raičević, Theuns & Lacey 2011; Lacey et al. 2011) have used galform (Cole et al. 2000; Baugh et al. 2005; Bower et al. 2006) combined with Monte Carlo merger trees. However, these studies calculated only the global evolution of reionization, and are not able to address the reionization structure. Most recently, Kim et al. (2013) have combined galform implemented within the Millennium-II dark matter simulation (Boylan-Kolchin et al. 2009), with a semi-numerical scheme to describe the resulting ionization structure. Kim et al. (2013) demonstrated the sensitivity of the ionization structure to the astrophysics of galaxy formation, and found that the strength of supernovae (SNe) feedback is the most important quantity.

In addition to the 21-cm power spectrum, several studies have previously analysed the cross-power spectrum (correlation) between redshifted 21-cm observations and galaxy surveys (Furlanetto & Lidz 2007; Lidz et al. 2009, 2011; Wiersma et al. 2013). These models showed that the cross-power spectrum should be observable, but do not provide a self-consistent link between the astrophysics of galaxy properties and the reionization structure. For example, Furlanetto & Lidz (2007) and Lidz et al. (2009, 2011) used a simple one-to-one relation between luminosity and dark matter halo mass. Conversely, in Wiersma et al. (2013), the cross-power spectrum was predicted using a semi-numerical code for 21-cm emission based on dark matter overdensity cross-correlated with a semi-analytic model for galaxies. As a result, the calculation did not include the direct relation between galaxies and ionization structure. In this paper our aim is to determine whether the cross-power spectrum can be used to infer the properties of high-redshift galaxy formation. We present predictions for the cross-power spectrum between 21-cm emission and galaxies using the model of Kim et al. (2013) who directly combined detailed models of high-redshift galaxy formation using galform with a semi-numerical description, and predict the resulting redshifted 21-cm power spectrum of different reionization histories. This model provides self-consistent results because the ionizing sources and observed galaxies are the same. These galaxies include both the observed luminous galaxies and the low-mass (∼10⁸ M_⊙) galaxies thought to drive reionization.

We begin in Sections 2 and 3 by describing the implementation of galform, our method for modelling the ionization structure, the cross-power spectrum and cross-correlation function, and the cross-correlation coefficient. The cross-power spectra from our method, and the effect of feedback processes on the cross-power spectra, are presented in Section 4. In Section 5 we describe the observational uncertainty. We finish with some conclusions in Section 6.

THE galform GALAXY FORMATION MODEL

In this section we summarize the theoretical galaxy formation modelling based on Kim et al. (2013) that is used in our analysis in order to describe the new features for this paper.

We implement the galform (Cole et al. 2000) model, within the Millennium-II dark matter simulation (Boylan-Kolchin et al. 2009). In this study, we specifically use the Lagos implementation of galform (Lagos et al. 2012) model described in Kim et al. (2013). The simulation has a cosmology including fractional mass and dark energy densities with values of Ω_m = 0.25, Ω_b = 0.045 and Ω_Λ = 0.75, a dimensionless Hubble constant of h = 0.73, and a power spectrum normalization of σ₈ = 0.9. The particle mass of the simulation is 6.89 × 10⁶ h⁻¹ M_⊙ and we detect haloes down to 20 particles (the minimum halo mass corresponds to ∼1.4 × 10⁸ h⁻¹ M_⊙) in the simulation box of side length L = 100 h⁻¹ Mpc.

Fig. 1 shows the relation between the UV magnitude (the rest-frame 1500 Å AB magnitude) including the effects of dust extinction of galaxies and the host halo mass (top), and between the total Lyman continuum luminosity (⁠|$\dot{N}_{\rm Lyc}$|⁠) of each galaxies and the host halo mass (bottom) from the galform model. Of particular note is that the luminosity of an ionizing source is not simply proportional to the host halo mass as is often assumed in reionization models (Lidz et al. 2009, 2011; Iliev et al. 2011). In part this is because of the distribution of satellite galaxies. The broad scatter of the relation indicates that physically motivated modelling for ionizing sources during the reionization should be included to understand the epoch of reionization. We note that this magnitude is not the same as ionizing luminosity. However, as shown in Fig. 1, the UV magnitude (the rest-frame 1500 Å AB magnitude) is closely related to the ionizing luminosity. Furthermore, Fig. 1 shows that the predicted ionizing luminosity to mass ratio from the model is not a simple one-to-one relation between luminosity and dark matter halo mass.

Figure 1.

The relation between the UV magnitude (the rest-frame 1500 Å AB magnitude) and the host halo mass (top), and between the total Lyman continuum luminosity (⁠|$\dot{N}_{\rm Lyc}$|⁠) of each galaxies (bottom) for galaxies and the host halo mass at z = 7.272 from galform. In each panel, black and blue dots represent central and satellite galaxies. Red, orange and yellow colours represent 1 (68.3 per cent), 2 (95.4 per cent) and 3-sigma (99.7 per cent) levels.

Open in new tabDownload slide

THE IONIZATION MODEL

In this section we summarize the calculation of the ionized structure (Section 3.1) and describe calculation of cross-power spectrum and cross-correlation function (Sections 3.2 and 3.3).

Semi-numerical scheme to calculate the evolution of ionizationed structure

Mesinger & Furlanetto (2007) introduced an approximate but efficient method for simulating the reionization process, referred to as a semi-numerical technique. In this paper we apply a semi-numerical technique to find the ionization structure resulting from galform galaxies within the Millennium-II dark matter simulation.

Based on equation (2), individual cells can have Q_cell ≥ 1. On the other hand, cells with Q_cell < 1 may be ionized by photons produced in a neighbouring cell. In order to find the extent of ionized regions we therefore filter the Q_cell field using a sequence of real space top hat filters of radius R (with 0.3906 < R < 100 h⁻¹ Mpc), producing one smoothed ionization field Q_R per radius. At each point in the simulation box we find the largest R for which the filtered ionization field is greater than unity (i.e. ionized with Q_R ≥ 1). All points within the radius R around this point are considered ionized. Ionization cells with 0 < Q_cell < 1 which are not part of an ionized Q_R ≥ 1 region retain their values.

The cross-power spectrum

The cross-correlation function

THE CORRELATION BETWEEN 21-CM EMISSION AND GALAXIES

In this section we present predictions for the cross-power spectrum, cross-correlation function and cross-correlation coefficient between 21-cm emission and galaxies as a function of redshift, luminosity and host halo mass (Section 4.1). We also discuss the effect of feedback processes on the cross-power spectrum, cross-correlation function and cross-correlation coefficient (Section 4.2).

Predictions for the correlation between 21-cm emission and galaxies

Fig. 2 shows the redshift evolution of the cross-power spectrum (top-left) and cross-correlation coefficient (bottom-left panel), and of the cross-correlation function (right panel) between redshifted 21-cm emission and galaxies. We show three examples which have UV magnitude limits, M_AB(1500 Å) − 5 log(h) < −18, in the model. This magnitude threshold corresponds to the deepest ‘wide’ area survey with Wide Field Camera 3/infrared and the Cosmic Assembly Near-Infrared Deep Extragalactic Legacy Survey on Hubble Space Telescope (Bouwens et al. 2011; Finkelstein et al. 2012). At each redshift, we calculate a mass-averaged ionization fraction, 〈x_i〉. From the correlation function, galaxies and 21-cm emission are anti-correlated at small separations while at large separations we find a weak correlation. These regions are separated by a transition wavenumber at which the cross-correlation coefficient and cross-correlation function change from negative to positive. Galaxies are correlated with 21-cm emission on scales larger than the ionized regions, but anti-correlated on smaller scales. The size of ionization regions therefore corresponds to this transition wavenumber. We find that the transition wavenumber from negative to positive cross-correlation coefficient increases as redshift decreases since the size of ionized regions generated by galaxies increases as the Universe evolves. We note that the cross-correlation function at z = 6.197 (〈x_i〉 = 0.95) has a shape that is different at small scales. We interpret this as being due to noise because the 21-cm emission regions are rare. In addition to the limitation of the integration range mentioned in Section 3.3, this noise indicates that a larger simulation box is necessary if one wants to measure the shape of the cross-correlation function at small k. Fig. 3 shows the comparison of cross-power spectra and cross-correlation functions between different host halo mass thresholds at z = 7.272 (〈x_i〉 = 0.55). We find that more massive haloes exhibit stronger anti-correlation as expected (Lidz et al. 2009; Wiersma et al. 2013). The same trend is also shown in Fig. 4 where we compare the results from calculations with different UV magnitude [M_AB(1500 Å) − 5 log(h)] thresholds. Fig. 4 shows that the transition wavenumber is similar for galaxy samples selected at different luminosity thresholds, since this scale is primarily set by the size of H ii regions.

Figure 2.

Redshift evolution of the cross-power spectrum and cross-correlation function between 21-cm fluctuations and the galaxies which have the UV magnitude less than −18 in the model. Left panel: the absolute value of the cross-power spectrum (top) and cross-correlation coefficient (bottom). Right panel: the corresponding cross-correlation function. In each panel, dotted (red), dash–three dotted (orange), dashed (yellow), dash–dotted (green), long-dashed (blue), and solid (purple) lines represent results from at z (〈x_i〉) = 9.278 (0.056), 8.550 (0.16), 7.883 (0.36), 7.272 (0.55), 6.712 (0.75) and 6.197(0.95), respectively.

Open in new tabDownload slide

Figure 3.

Comparison of the cross-power spectrum and cross-correlation function for different host halo mass thresholds at z = 7.272 (〈x_i〉 = 0.55). Left panel: the absolute value of the cross-power spectrum (top) and cross-correlation coefficient (bottom). Right panel: the corresponding cross-correlation function. In each panel, the dotted (dark brown), dot–dashed (brown), dashed (orange), long-dashed (yellow) lines show the cross-correlation using galaxies which are included in 10⁹, 10¹⁰, 10¹¹ and 10¹²h⁻¹ M_⊙, respectively.

Open in new tabDownload slide

Figure 4.

The same as Fig. 3 but results are computed based on different UV magnitude thresholds. In each panel, the dotted (black), dot–dashed (blue) and long dashed (sky-blue) lines show the cross-correlation using galaxy samples which are, respectively, more luminous than magnitude limits of −10, −15 and −20.

Open in new tabDownload slide

The effect of feedback processes

In order to investigate the effect on the power spectrum of different feedback processes in galaxy formation, we follow a similar method to Kim et al. (2013). We use the Lagos et al. (2012) galaxy formation model as our fiducial case, and then consider two variants of this (hereafter called NOSN models) which have SNe feedback turned off. We use two variants of the NOSN model. First, we consider the inclusion of photoionization feedback using V_cut = 30 km s⁻¹, where V_cut is a threshold value of the host halo's circular velocity (Kim et al. 2013). Secondly, we removed both SNe feedback and photoionization feedback by setting V_cut = 0 km s⁻¹. We refer to this second model as NOSN (no suppression) in this paper. Since turning off SNe feedback in the Lagos et al. (2012) model changes the bright end of the UV luminosity function, we have changed some other parameters so that the NOSN models still match the observed UV luminosity functions at z = 7.272. Specifically, we introduce a stellar initial mass function dominated by brown dwarfs, with ϒ = 4, and also reduce the star formation time-scale in bursts by setting f_dyn = 2 and τ_{★burst, min} = 0.005 Gyr [see Cole et al. (2000) and Lacey et al. (2011) for more details of these parameters]. In Fig. 5 we show the resulting comparison of cross-power spectra and cross-correlation functions at z = 7.272 (〈x_i〉 = 0.55). The locations of transition wavenumbers between the Lagos et al. (2012) model and the two NOSN models are significantly different (see also ionization structure for these models in Kim et al. 2013). In particular, the Lagos et al. (2012) model has a larger transition scale. On small scales, the cross-correlation function of the Lagos et al. (2012) model shows stronger anti-correlation than the two NOSN models between 21-cm emissions and galaxies.

Figure 5.

The same as Fig. 4 but results are computed based on different feedback processes. In each panel, solid (red), dot–dashed (light grey) and long dashed (dark grey) lines represent our model, NOSN(V_cut = 30 km ;s⁻¹), and NOSN(no suppression) models, respectively.

Open in new tabDownload slide

We note that for the three different models, we have used the same mass-averaged ionization fraction, 〈x_i〉, at each redshift (listed in Section 3.1). While these models are forced to have the same ionization history, Fig. 5 still shows different cross-power spectrum shapes. This is because of the effect of feedback processes. SN feedback suppresses the formation of galaxies within small dark matter haloes, and consequently this process enhances the galaxy bias of the ionizing emission. Photoionization feedback also suppresses the formation of low-luminosity galaxies, but the effect is not significant compared to the effect of SNe feedback (see also Kim et al. 2013). This weak effect of photoionization feedback is revealed by small difference between two NOSN models.

DETECTABILITY

In this section we describe the error estimation of the cross-correlation coefficient (Section 5.1) and discuss observational requirements for future galaxy surveys (Section 5.2). Our examples are based on the MWA-like observations of the 21-cm signal combined with various hypothetical galaxy redshift surveys.

Error estimate in the cross-correlation coefficient

The first term in equation (15) comes from a sample variance within the finite volume of the survey and the second term comes from the thermal noise of the 21-cm telescope. We have assumed specifications of the MWA for the calculation of thermal noise. In the thermal noise term, T_sys ∼ 250[(1 + z)/7]^2.6K denotes the system temperature of the telescope; B = 8 MHz is the survey bandpass; t_int is the integration observing time. We use 1000 h total observing time in this calculation; D and ΔD are the comoving distance to the survey volume and the comoving survey depth, |$\Delta D = 1.7(\frac{B}{0.1 {\rm MHz}})\sqrt{\frac{1+z}{10}}(\frac{\Omega _{\rm m}\,h^2}{0.15})^{-1/2}$| (Furlanetto, Oh & Briggs 2006), respectively; n(k_⊥) denotes the number density of baselines in observing the transverse component of the wave vector, where |$k_{\perp } = \sqrt{1-\mu ^2}k$|⁠. Observing the signal of k_⊥ in each Fourier cell is related to the length of baseline and the antenna configuration. Here, we follow the method of Morales (2005), Bowman, Morales & Hewitt (2006) and Datta, Bharadwaj & Choudhury (2007) for calculation of n(k_⊥). The maximum value of the transverse component of the wave vector is k_{⊥, max} = 2πL_max/(Dλ), where L_max = 750 m is the maximum baseline distance in the antenna array. This limit is due to the maximum angular resolution of the telescope related to L_max. On the other hand, the minimum line-of-sight wavenumber is set by the bandpass k_min = 2π/ΔD; The observed wavelength is λ = 0.21 m × (1 + z), and A_e is the effective collecting area of each antenna. We use A_e ∼ N_dipλ²/4 (Bowman et al. 2006), where N_dip = 16 is the number of dipoles. We have assumed 500 antenna elements.¹

Fig. 6 shows the 21-cm power spectrum with errors estimated based on equation (15) for cases including different feedback processes. The 21-cm power spectra show obvious differences between the models for SNe feedback, especially at large scales. Fig. 6 reinforces the importance of detailed modelling of galaxy formation during reionization (Kim et al. 2013).

Figure 6.

The 21-cm power spectrum with estimated errors, based on an 800 deg² survey area, at z = 7.272. We assume 1000 h total observing time, and, based on the assumption of 8 MHz bandwidth, the survey depth is about 0.2 redshift units. Red represents the power spectrum from our model including SNe feedback with V_cut = 30 km s⁻¹. The light grey and dark grey lines represent power spectrum from the NOSN models with V_cut = 30 km s⁻¹ and no suppression, respectively.

Open in new tabDownload slide

The error on the galaxy power spectrum is expressed in equation (16). The galaxy shot-noise is dependent on the number density of galaxies observable (n_gal), k_∥ = μk, and σ_χ = cσ_z/H(z), where σ_z is the galaxy redshift error. Here, we assume a Gaussian distribution of redshift errors.

Observational requirements for future galaxy surveys

Figure 7.

The S/N for the cross-correlation coefficient as a function of survey area and relative redshift error at z = 7.272. Left panel: plots of S/N as a function of survey area (A_survey) for different models. We assume σ_z = 0.05. Central panel: a plot of S/N as a function of redshift error, σ_z, with A_survey = 5 deg² for the default model. In the left and central panels, solid (brown), long-dashed (light grey) and dotted (dark grey) lines represent results from our model including SNe feedback with V_cut = 30 km s⁻¹ (the default model), the NOSN models with V_cut = 30 km s⁻¹ and no suppression, respectively. Right panel: plots of S/N as a function of σ_z at different wavenumbers for default model. Dotted, dashed, dot–dashed, three dot–dashed, long dashed and solid lines represent k = 0.1, 0.2, 0.4, 0.6, 0.8 and 1.0 h Mpc⁻¹, respectively. In each panel, we assume 1000 h total observing time.

Open in new tabDownload slide

As a specific example, we also calculate the total S/N as a function of σ_z by assuming the survey area of A_survey = 5 deg². The total S/N in the central panel of Fig. 7 shows that measurements will require redshift uncertainties less than 0.1. The NOSN models show a similar shape to the default model, but have lower S/N. Lower accuracy redshifts (σ_z > 0.1) wash out the cross-correlation signal. An error of σ_z ∼ 0.1 provides measurement only on larger scales (k[h Mpc⁻¹] < 0.2) (right-hand panel in Fig. 7). To measure the cross-correlation over a broad range of k, redshift uncertainties, σ_z, less than 0.01 will be required.

Fig. 7 illustrates the conditions required for measurement of the 21-cm-galaxy cross-correlation. Before concluding we discuss these requirements with respect to real galaxy surveys. In this study, we have used UV magnitude cuts to select galaxy samples which relate to observed Lyman-break galaxies. However, Lyman-break galaxies, which are photometrically selected, have σ_z ≳ 0.5 at z ∼ 6.5 (Beckwith et al. 2006) much longer than the σ_z < 0.1 requirement. As a result, Lyman-break surveys will not be sufficient to detect the cross-correlation (Wiersma et al. 2013). On the other hand, Lyα emitters selected from narrow-band surveys have σ_z ∼ 0.05–0.1 (Ouchi et al. 2008, 2010). Thus, a detection could be made based on the precision and volume of current Lyα surveys. Our semi-numerical model does not predict Lyα luminosity (see Orsi et al. 2008). However, the difference between simulated populations of Lyα emitters and Lyman-break galaxies is found not to be significant (Dayal & Ferrara 2012). For the purpose of our calculation, we therefore use star-forming galaxies with UV magnitudes [M_AB(1500 Å) − 5 log (h)] corresponding to the number density of n_gal = 1.6 × 10⁻⁴ h³ Mpc⁻³, which is seen in the Subaru Deep Field at z ∼ 6.6 (Kashikawa et al. 2006).

The largest Lyα survey (Ouchi et al. 2010) covered 1 deg² at z ∼ 6.6, and used a narrow-band filter with a central wavelength of 9196 Å and a full width at half-maximum of 132 Å. These values correspond to a survey depth of Δz ∼ 0.11 at z = 6.6. This is smaller than, but comparable to, the survey depth of MWA observation which is Δz ∼ 0.3 corresponding to the bandwidth of 8 MHz assumed for this paper. For the survey at z ∼ 7.3, Shibuya et al. (2012) have a survey depth of Δz ∼ 0.18, which used the central wavelength of 10052 Å and a full width at half-maximum of 214 Å. This is also smaller than, but comparable to, the MWA observation of the survey depth of Δz ∼ 0.38.

While 1 deg² represents the largest high-redshift survey at the current time, future surveys will be larger. For example, in the next 5 years Hyper Suprime-Cam on the Subaru telescope will observe 10⁵ galaxies at z ∼ 5.7 and 6.5 in a survey area of ∼30 deg², and 100 s of galaxies at z ∼ 7 in a survey area of 3 deg² (M. Ouchi, private communication). As shown in Fig. 7, this increased survey area will improve the S/N, so that the cross-power spectra signal could be detected with high significance.²

Based on the requirement from Section 5.2, we assume a 5 deg² galaxy survey field and a redshift error of 0.05 for a future galaxy survey. We also assume 1000 h total observing time. Fig. 8 shows the predicted errors for the cross-correlation coefficient within spherical bins of logarithmic width ϵ = 0.5 at z = 6.712 and 7.272 for such a galaxy survey combined with the MWA. The estimated errors are exponentially increased near the wavenumber of 1 h Mpc⁻¹, because of the limit of k_{⊥, max} for a 21-cm survey. We compare the result with the cross-correlation coefficient from the NOSN models. The result shows that we could observationally distinguish our default model from two different NOSN models.

Figure 8.

The cross-correlation coefficient at z = 6.712 (left panel) and 7.272 (right panel). The error bars are calculated for spherical bins of logarithmic width ϵ ≡ d lnk = 0.5. We assume a 5 deg² galaxy survey field, 1000 h total observing time, the redshift error of 0.05 and galaxy number density of the Subaru Deep Field survey. The solid (orange and red) lines represent the power spectrum from the default model (Lagos et al. 2012) including SNe feedback with V_cut = 30 km s⁻¹. The long dashed (light grey) and dot–dashed (dark grey) lines represent power spectrum from the NOSN models with V_cut = 30 km s⁻¹ and no suppression, respectively.

Open in new tabDownload slide

Lyα observations at z ≳ 7 over a large area are very challenging. The latest Lyα survey at z ∼ 7.3 (Shibuya et al. 2012) has a galaxy number density of ∼6.7 × 10⁻⁶ covering a survey area of 0.48 deg². This value is smaller than the value we assume for our error estimation. Computing the cross-power spectrum corresponding to this number density is not possible owing to the limited box size of our simulation. However, we have checked that the estimated error would approximately increase by a factor of 2, when using this number density.

SUMMARY AND CONCLUSIONS

In this study we have investigated evolution of the cross-power spectrum, cross-correlation function and cross-correlation coefficient between 21-cm emission and galaxies using the model of Kim et al. (2013). This model combines the hierarchical galaxy formation model galform implemented within the Millennium-II dark matter simulation, with a semi-numerical scheme to describe the resulting ionization structure. We find that there is a transition wavenumber, k, at which the cross-correlation coefficient changes from negative to positive (Lidz et al. 2009). This transition wavenumber is associated with the size of the ionized regions generated by galaxies, and increases with decreasing redshift. We also find the same trend in the cross-power spectrum and cross-correlation function. We calculated the cross-power spectrum as a function of UV luminosity [M_AB(1500 Å) − 5 log(h)] and host halo mass. These calculations reveal that bright galaxies and galaxies residing in massive haloes have stronger anti-correlation, but a similar transition wavenumber.

We have studied observational uncertainties in measurement of the cross-correlation coefficient based on the specifications of an upgraded (512 tile) Murchison Widefield Array (MWA) combined with galaxy surveys. The results show that the cross-power spectrum signal could be detected when combined with more than 3 deg² of a galaxy survey at the depth of the future galaxy survey having redshift error <0.1. We have also investigated the dependence on the inclusion of feedback processes in the galaxy modelling. We find that the amplitude of the cross-correlation is larger when SNe feedback is considered and that the cross-correlation coefficient has a different shape compared to a model with no SNe feedback. Thus the cross-correlation could be used to determine the importance of SNe feedback in high-redshift galaxies. Our results imply that detailed modelling of reionization processes and galaxy formation is required to predict an accurate cross-correlation between 21-cm emission and galaxies, and to interpret future observational measurements.

HSK is supported by a Super-Science Fellowship from the Australian Research Council. JSBW acknowledges the support of an Australian Research Council Laureate Fellowship. The Centre for All-sky Astrophysics is an Australian Research Council Centre of Excellence, funded by grant CE110001020. This work was supported in part by the Science and Technology Facilities Council rolling grant ST/I001166/1 to the ICC. The Millennium II Simulation was carried out by the Virgo Consortium at the supercomputer centre of the Max Planck Society in Garching. Calculations for this paper were partly performed on the ICC Cosmology Machine, which is part of the DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and Durham University.

REFERENCES