# The thermal Sunyaev–Zel'dovich effect from the epoch of reionization

Ilian T Iliev, Azizah R Hosein, Jens Chluba, Luke Conaboy, David Attard, Rajesh Mondal, Kyungjin Ahn, Stefan Gottlöber, Joseph S W Lewis, Pierre Ocvirk, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Ilian T Iliev, Azizah R Hosein, Jens Chluba, Luke Conaboy, David Attard, et al. The thermal Sunyaev–Zel'dovich effect from the epoch of reionization. Monthly Notices of the Royal Astronomical Society, 2025, 540 (2), pp.1432 - 1448. 10.1093/mnras/staf776. insu-05106141

HAL Id: insu-05106141

https://insu.hal.science/insu-05106141v1

Submitted on 10 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The thermal Sunyaev–Zel’dovich effect from the epoch of reionization

Ilian T. Iliev<fnref n="1" />,<fnref n="1" />* Azizah R. Hosein,<fnref n="1" /> Jens Chluba<fnref n="2" />,<fnref n="2" /> Luke Conaboy<fnref n="2" />,<fnref n="1" />,<fnref n="3" /> David Attard,<fnref n="1" />

Rajesh Mondal,<fnref n="4" /> Kyungjin Ahn,<fnref n="5" /> Stefan Gottlöber,<fnref n="6" /> Joseph S.W. Lewis,<fnref n="7" /> Pierre Ocvirk,<fnref n="8" />

Hyunbae Park,<fnref n="9" /> Paul R. Shapiro,<fnref n="10" /> Jenny G. Sorce<fnref n="11" />,<fnref n="12" /> and Gustavo Yepes<fnref n="13" />,<fnref n="14" />

Accepted 2025 April 15. Received 2025 April 10; in original form 2024 December 6

## ABSTRACT

The thermal Sunyaev–Zel’dovich (tSZ) effect arises from inverse Compton scattering of low energy photons on to thermal electrons, proportional to the integrated electron pressure, and is usually observed from galaxy clusters. We can expect that the epoch of reionization (EoR) also contributes to this signal, but this has not been previously evaluated. In this work, we analyse a suite of fully coupled radiation–hydrodynamics simulations based on RAMSES-CUDATON to study the EoR tSZ signal. We construct lightcones of the electron pressure in the intergalactic medium for $6 \lesssim z$ to calculate the resulting Compton y-parameters. We vary the box sizes, resolutions, and star formation parameters to investigate how these factors affect the tSZ effect. We produce plots of maps and distributions of y, as well as angular temperature power spectra of the tSZ signal obtained from integrating the lightcones constructed for each simulation. We find that the tSZ signal from reionization is generally sub-dominant to the post-reionization one at larger scales ($\ell < 10^4$), but can contribute non-trivially and potentially contaminate the measured signals. At scales probed by current experiments like SPT ($\ell \sim 10^3 - 10^4$), we find that the tSZ signal power spectrum from reionization contributes at roughly a per cent level compared to the current templates, with the quadratic Doppler effect contributing an additional $\sim 10$ per cent to the tSZ signal. At smaller scales the tSZ from reionization peaks and can potentially dominate the total signal and is thus a potentially much more important contribution to take into account in any future, more sensitive experiments.

**Key words:** radiative transfer – methods: numerical – intergalactic medium – cosmic background radiation – dark ages, reionization, first stars – large-scale structure of Universe.

## 1 INTRODUCTION

The epoch of reionization (EoR) marks the last global phase transition of the intergalactic medium, from cold and largely neutral state to the hot and fully ionized one that we observe today. It is the period in cosmic history during which the first luminous objects formed, driving this process. This process lasted for up to a billion years and was likely very patchy, with some volumes around the first sources

ionized very early, and the deep voids far from any sources finishing this much later.

Although ongoing and upcoming surveys (e.g. with the Low Frequency Array (LOFAR), Square Kilometre Array (SKA), and JWSZ) will give us considerably more insight into this epoch, leading into its first direct detection, much still remains unknown. Low-frequency radio surveys will allow us to obtain topographical images of the neutral hydrogen gas from the emitted 21-cm signal, while infrared imaging will shed more light on bright luminous sources at much higher redshifts than are currently available. Additionally, the EoR imprints information into the cosmic microwave background (CMB).

<footnote n="1">¹Department of Physics & Astronomy, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote n="2">²Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester M13 9PL, UK</footnote>

<footnote n="3">³School of Physics and Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD, UK</footnote>

<footnote n="4">⁴Department of Physics, National Institute of Technology Calicut, Calicut 673601, Kerala, India</footnote>

<footnote n="5">⁵Department of Earth Sciences, Chosun University, Gwangju 61452, Republic of Korea</footnote>

<footnote n="6">⁶Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany</footnote>

<footnote n="7">⁷Institut d'Astrophysique de Paris, UMR 7095, CNRS, Sorbonne Université, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="8">⁸Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="9">⁹Center for Computational Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan</footnote>

<footnote n="10">¹⁰Department of Astronomy and Texas Center for Cosmology and Astroparticle Physics, The University of Texas at Austin, Austin, TX 78712-1083, USA</footnote>

<footnote n="11">¹¹Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRISTAL, F-59000 Lille, France</footnote>

<footnote n="12">¹²Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France</footnote>

<footnote n="13">¹³Departamento de Física Teórica, Modulo 8, Facultad de Ciencias, Universidad Autónoma de Madrid, 28049 Madrid, Spain</footnote>

<footnote n="14">¹⁴CIAFF, Facultad de Ciencias, Universidad Autónoma de Madrid, 28049 Madrid, Spain</footnote>

<footnote>*E-mail: I.T.Iliev@sussex.ac.uk</footnote>
While the primary anisotropies of the CMB store information about the primordial density fluctuations at last scattering, the secondary anisotropies are due to interactions between the CMB photons and matter along their paths to us. Therefore, the latter can shed light on the structures and physics in the Universe since the last scattering. One of the contributors to these secondary CMB anisotropies is the Sunyaev–Zel’dovich (SZ) effect.

The SZ effect (Sunyaev & Zeldovich 1972, 1980a,b) (see also reviews by Rephaeli 1995; Birkinshaw 1999; Carlstrom, Holder & Reese 2002; Kitayama 2014; Mroczkowski et al. 2019) arises from the inverse Compton scattering of CMB photons by energetic electrons. This transfers energy to the photons and could distort the CMB spectrum from that of a pure blackbody. There are two types of SZ effect: thermal (tSZ) and kinetic (kSZ). The tSZ effect is caused by thermally moving electrons, while the kSZ effect arises from large-scale bulk motions of electrons. The tSZ effect changes both the energy and the spectrum of the radiation, while kSZ leaves the spectrum as blackbody, but shifted to higher energies. An example of both these effects can be seen in galaxy clusters – when photons travel through a cluster, they are scattered by the electrons in the hot intracluster gas ($T \sim 10^7 - 10^8$ K) which are moving both thermally (tSZ) and with the overall bulk motion of the cluster itself (kSZ).

The kSZ effect has been proposed as a probe of the EoR, both on its own (see e.g. Iliev et al. 2007, 2008; Mesinger, McQuinn & Spergel 2012; Park et al. 2013), and in terms of its cross-correlation with the redshifted 21-cm emission (e.g. Jelić et al. 2010; Tashiro et al. 2010). The reionization process creates ionized patches, whose bulk motions yield fluctuating kSZ secondary anisotropies at the scales corresponding to the patch sizes. These works showed that EoR has appreciable contribution to the small-scale ($\ell > 3000$) CMB anisotropies.

In contrast, to date the tSZ effect has largely been used for detection and studies of galaxy clusters at lower redshift, but has not been not considered in EoR context. Due to its characteristic spectral signature (the tSZ effect lowers the intensity of the photons at frequencies <218 GHz, and raises it at frequencies >218 GHz), the tSZ effect is used as redshift-independent approach to detect clusters. It is also useful for measuring the thermal energy of the intercluster gas and measurements of the Hubble constant. The tSZ effect is also one of the foregrounds which contaminate CMB measurements. Templates of its angular power spectrum are thus used for the purpose of cleaning the primordial CMB signal (e.g. Shaw et al. 2010).

The tSZ effect is a measure of the integrated electron pressure along the line of sight. During the EoR the galaxy groups and clusters have not yet formed, and the typical electron gas temperatures caused by photoionization are of order tens of thousands of K, rather than the million-degree temperatures of intercluster gas. However, due to the effects of Hubble expansion ($\bar{\rho} \propto (1+z)^3$), the matter density of both the IGM and haloes are much higher than it is later on, resulting in correspondingly higher gas pressure. We can therefore reasonably expect that the EoR-produced electron pressure may provide a non-trivial contribution to the total observed effect. Furthermore, since the tSZ effect from EoR has some contribution from the diffuse intergalactic medium (IGM) in the ionized patches, the spatial structure of the temperature fluctuations should be different from the one produced by galaxy clusters, potentially providing a new interesting probe of cosmic reionization. The aim of this paper is to investigate the strength of this effect and its detectability. A significant contribution from this previously neglected contribution to tSZ would have consequences on cluster measurements and would require updating models used to create templates for the tSZ angular

power spectrum. Several upcoming CMB experiments like NICA2¹ and C4² will reach significantly higher resolutions, reaching the small-scale regime of interest, where the reionization contribution to the anisotropies is expected to peak.

We investigate the tSZ signal arising from the EoR by analysing the hydrodynamical data produced by a series of fully coupled radiation–hydrodynamics simulations ran using the RAMSES-CUDATON code (Teyssier 2002; Aubert & Teyssier 2008). Our baseline simulation is Cosmic Dawn II (CoDa II; Ocvirk et al. 2020). The results of further six auxiliary simulations with varying volumes, and resolutions are used to investigate the effects of these variations and of varying star formation parameters on the tSZ signal.

The rest of this paper is organized as follows. We first describe the physics behind the tSZ effect (Section 2) and the simulations used in this project (Section 3). We then outline our methodology (Section 4) and present our results (Section 5) and end with our conclusions (Section 6). The data analysis for this work was performed in part using the SEREN3 PYTHON package.<fnref n="3" />

## 2 THE SUNYAEV–ZEL’DOVICH EFFECT

The SZ effect arises from the inverse Compton scattering of CMB photons by free electrons. When CMB photons travel through a cloud of such electrons, the probability that they will scatter is dictated by the Thomson scattering optical depth,

$$ \tau_e = \sigma_T \int n_e \, dl \sim 2 \times 10^{-3} \left( \frac{n_e}{10^{-3} \, \text{cm}^{-3}} \right) \left( \frac{l}{\text{Mpc}} \right), \quad (1) $$

where the integral is performed along the line of sight, $n_e$ is the electron number density, and the fiducial values denoted are typical for galaxy clusters, where the SZ effect is most commonly observed. Although inverse Compton scattering occurs in a variety of scenarios, the SZ effect usually refers to the scattering of CMB photons in the GHz to THz range of frequencies on non- or mildly relativistic electrons.

### 2.1 The thermal Sunyaev–Zel’dovich effect

When CMB photons pass through the hot intracluster medium (ICM), there is a ~1 per cent chance that it will interact with one of the energetic electrons in the plasma. The scattered photon experiences an energy boost of ~$4 \frac{k_B T_e}{m_e c^2}$ (see e.g. Rybicki & Lightman 1979; Sazonov & Sunyaev 2000), causing a distortion in CMB intensity given by

$$ \Delta I_v \approx I_0 y \frac{x^4 e^x}{(e^x - 1)^2} \left( x \frac{e^x + 1}{e^x - 1} - 4 \right) \equiv I_0 y g(x), \quad (2) $$

where $x = \frac{h_P v}{k_B T_{CMB}} \approx \frac{v}{56.8 \text{ GHz}}$ is the dimensionless frequency, y is the Compton y-parameter (defined below), and

$$ I_0 = \frac{2 (k_B T_{CMB})^3}{(h_P c)^2} = 270.33 \left( \frac{T_{CMB}}{2.73 \text{ K}} \right)^3 \text{MJy sr}^{-1}. \quad (3) $$

If $\frac{\Delta I_v}{I_v} \ll 1$, the signal can be expressed in terms of the CMB temperature, using the derivative of the Planck function with respect

<footnote n="1">¹https://www.epj-conferences.org/articles/epjconf/pdf/2024/03/epjconf_mmUniverse2023_00040.pdf</footnote>

<footnote n="2">²https://cmb-s4.org</footnote>

<footnote n="3">³https://github.com/sully90/seren3</footnote>
to temperature:

$$ \frac{\Delta T_{\text{CMB}}}{T_{\text{CMB}}} \approx y \left( x \frac{e^x + 1}{e^x - 1} - 4 \right) = yf(x), \quad (4) $$

where $f(x)$ is the frequency dependence of the tSZ spectrum in terms of $\Delta T_{\text{CMB}}$. When there are electrons with energies in the relativistic regime, a factor $\delta_{\text{SZ}}(x, T_e)$ is incorporated, and $f(x)$ becomes

$$ f(x) = \left( x \frac{e^x + 1}{e^x - 1} - 4 \right) (1 + \delta_{\text{SZ}}(x, T_e)). \quad (5) $$

For more detailed interpretations of the relativistic corrections to the tSZ effect, see e.g. Challinor & Lasenby (1998), Nozawa et al. (2006), and Chluba et al. (2013). We will neglect such effects here. In the non-relativistic case and in the Rayleigh–Jeans, low-energy limit, we have $f(x \ll 1) \to -2$, while at high frequencies, $f(x \gg 1) \to x - 4$.

The spectral signature of the tSZ can be separated from the kSZ effect around the null of the tSZ effect, and other temperature fluctuations due to its characteristic frequency dependence above. Hence, multifrequency measurements are necessary to distinguish all these effects. Additionally, the tSZ signal dominates over the kSZ effect roughly by an order of magnitude for clusters. This is due to the thermal velocity of electrons ($\sim 10^4$ km s⁻¹) being much higher than the bulk velocity ($\lesssim 10^3$ km s⁻¹).

We can see from equation (4) that the change in the CMB temperature is proportional to the Compton y-parameter, which is thus used to measure the magnitude of the tSZ signal. It is defined as the line of sight integral of the electron pressure,

$$ y \equiv \frac{k_B}{m_e c^2} \int T_e \, d\tau_e = \frac{k_B \sigma_T}{m_e c^2} \int n_e T_e \, dl = \frac{\sigma_T}{m_e c^2} \int p_e \, dl, \quad (6) $$

where $T_e$ is the temperature of the electrons, and $p_e$ is the pressure due to the electrons. It is important to note that the Compton y-parameter is redshift-independent. This can also be explained by considering that $\Delta T_{\text{CMB}}$ (and $\Delta I_v$) is redshifted the same way as $T_{\text{CMB}}$ (and $I_v$). Hence, the tSZ effect does not experience a loss in intensity with redshift (dimming), making it a useful tool for measuring large-scale structure in the universe. The ICM of typical clusters have electron temperatures of 5–10 keV. For massive clusters with central optical depth $\sim 10^{-2}$ we have $y \sim 10^{-4}$. More generally, the y-parameter for clusters is typically $y \gtrsim 10^{-5}$. In practice, all clusters above certain mass cut-off in a given area of sky can be detected with the tSZ effect using CMB frequency maps of the y-parameter.

Another useful observable is the angular power spectrum of the tSZ signal. Komatsu & Seljak (2002) found that the tSZ angular power spectrum is cosmology-dependent, with $D_l \propto \sigma_8^7 \Omega_b^2 h^2$. Later Shaw et al. (2010) used numerical simulations to refine this scaling to $D_l \propto \sigma_8^{8.3} \Omega_b^{2.8} h^{1.7}$. This dependence makes the tSZ power spectrum a powerful tool for constraining these cosmological parameters (see e.g. Barbosa et al. 1996; Sievers et al. 2013; Crawford et al. 2014; Hill et al. 2014; George et al. 2015; Planck Collaboration XXII 2015; Horowitz & Seljak 2017). The tSZ power spectrum is also affected by astrophysical processes, such as active galactic nucleus (AGN) feedback, shock heating, radiative cooling, etc. Komatsu & Seljak (2002), Battaglia et al. (2010), Sehgal et al. (2010), Shaw et al. (2010), and Trac, Bode & Ostriker (2011) have produced templates for the (post-reionization) tSZ power spectrum, fitted for Wilkinson Microwave Anisotropy Probe (WMAP) Lambda cold dark matter (ΛCDM) cosmology, while probing its dependence on cosmological and astrophysical effects. While these studies investigate the tSZ signal in the context of galaxy clusters, this is not the only possible source of such signal. The aim of our current study is to determine

the tSZ effect arising from the EoR, which previously has largely been neglected.

It is instructive to consider the mean values of the y-parameter, $\langle y \rangle$ produced by various large-scale cosmological structures. The tightest constraint on overall $\langle y \rangle$ was done by the COBE Far Infrared Absolute Spectrophotometer (FIRAS) experiment and found $\langle y \rangle < 1.5 \times 10^{-5}$ at the 95 per cent confidence level (Fixsen et al. 1996). Refregier et al. (2000) performed hydrodynamic simulations comprising only gravitational forces (with no star formation) as well as analytic calculations with the Press–Schechter (PS) formalism (Press & Schechter 1974) to compute the tSZ signal. For ΛCDM cosmology, they found $\langle y \rangle = 1.67 \times 10^{-6}$ (simulation) and $2.11 \times 10^{-6}$ (PS), both values being about an order of magnitude below the COBE-FIRAS upper limit. Their projected maps of the y-parameter showed clusters having $y < 10^{-5}$, and groups and filaments having $y \sim 10^{-7}-10^{-5}$. However, they note that projecting a number of simulation boxes along the line of sight on the sky would cause the filamentary objects to be averaged out (da Silva et al. 2001; Seljak, Burwell & Pen 2001). They also find that the majority of the tSZ signal arises from low redshifts ($z < 2$). Their angular power spectra showed the tSZ effect having comparable power to the primary CMB at $l \sim 2, 000$, while groups and filaments contributed $\sim 50$ per cent of the power at $l = 500$ with $\sim 50$ per cent of that power being produced at $z \lesssim 0.1$.

By combining hydrodynamic simulations and analytic models, Zhang, Pen & Trac (2004) estimated

$$ \langle y \rangle = 2.6 \times 10^{-6} \left( \frac{\sigma_8}{0.84} \right)^{4.1-2\Omega_m} \left( \frac{\Omega_m}{0.268} \right)^{1.28-0.2\sigma_8} \quad (7) $$

for a flat ΛCDM WMAP cosmology, with the dominant contribution coming from $z \sim 1$. More recently, Hill et al. (2015) used analytic calculations [including an ICM model (Hill et al. 2014), relativistic corrections (Arnaud, Pointecouteau & Pratt 2005), and a reionisation model (Battaglia et al. 2013)] to compute the total mean Compton parameter of the universe. They found $\langle y \rangle_{\text{ICM}} = 1.58 \times 10^{-6}$, $\langle y \rangle_{\text{IGM}} = 8.9 \times 10^{-8}$ and $\langle y \rangle_{\text{EoR}} = 9.8 \times 10^{-8}$ for the contributions from the ICM, IGM, and EoR, respectively, which gave a total mean value of $\langle y \rangle_{\text{total}} = 1.77 \times 10^{-6}$. The signals from the IGM and EoR are thus expected to be sub-dominant to the ICM signal. However, none of these results is based on full radiative-hydrodynamics EoR simulations, and can thus might not be fully accurate. Furthermore, even if this EoR tSZ signal is sub-dominant, its accurate estimation would be very useful for providing constraints on galaxy formation models, feedback mechanisms, and the thermal history of the universe. It could also be a possible source of errors for cluster tSZ measurements.

## 2.2 Second-order doppler distortions

Zeldovich, Illarionov & Syunyaev (1972) showed that the spectral distortions produced by thermally energetic electrons give rise to a spectrum which can be described as a superposition of blackbodies with a Compton y-parameter, $y = O(v^2)$. This y-parameter depends on the second-order of the electron velocities, $v$, essentially making it the equivalent of a tSZ effect due to the bulk flow. Hu, Scott & Silk (1994) investigated the significance of this quadratic Doppler distortion during the EoR. They showed that when this contribution to the y-parameter is included, equation (6) becomes

$$ y = \frac{\sigma_T}{m_e c^2} \int \left( \frac{1}{3} m_e n_e \langle v^2 \rangle + p_e \right) dl. \quad (8) $$
<table><caption>**Table 1.** Summary of the parameters used in the simulations. A few parameters were not varied between our simulations: massive star lifetime $t_{\star} = 10$ Myr; supernova energy $E_{SN} = 10^{51}$ erg; effective photon energy of 20.28 eV; effective H I cross-section $\sigma_{E} = 2.493 \times 10^{-22}$ m$^2$; and (full) speed of light $c = 299792458$ m s$^{-1}$. CoDa II did not use a temperature threshold criterion for star formation.</caption><thead><tr><th>Simulation set</th><th>CoDa II</th><th>100 h<sup>-1</sup> Mpc</th><th>50 h<sup>-1</sup> Mpc</th><th>25 h<sup>-1</sup> Mpc</th><th>10 h<sup>-1</sup> Mpc</th></tr></thead><tbody><tr><td>Number of nodes (GPUs)</td><td>16 384</td><td>Set-up</td><td></td><td></td><td></td></tr><tr><td>Grid size</td><td>4096<sup>3</sup></td><td>128</td><td>128</td><td>32</td><td>128</td></tr><tr><td>Comoving box size L<sub>box</sub></td><td>64 h<sup>-1</sup> Mpc</td><td>1024<sup>3</sup></td><td>1024<sup>3</sup></td><td>512<sup>3</sup></td><td>1024<sup>3</sup></td></tr><tr><td>Comoving force resolution dx</td><td>23.1 kpc</td><td>145.5 kpc</td><td>72.8 kpc</td><td>72.8 kpc</td><td>14.6 kpc</td></tr><tr><td>Physical force resolution at z = 6</td><td>3.3 kpc</td><td>20.8 kpc</td><td>10.4 kpc</td><td>10.4 kpc</td><td>2.1 kpc</td></tr><tr><td>DM particle number N<sub>DM</sub></td><td>4096<sup>3</sup></td><td>1024<sup>3</sup></td><td>1024<sup>3</sup></td><td>512<sup>3</sup></td><td>1024<sup>3</sup></td></tr><tr><td>DM particle mass M<sub>DM</sub></td><td>4.07 &times; 10<sup>5</sup>M<sub>&odot;</sub></td><td>1.05 &times; 10<sup>8</sup>M<sub>&odot;</sub></td><td>1.31 &times; 10<sup>7</sup>M<sub>&odot;</sub></td><td>1.31 &times; 10<sup>7</sup>M<sub>&odot;</sub></td><td>1.05 &times; 10<sup>5</sup>M<sub>&odot;</sub></td></tr><tr><td>Initial redshift z<sub>start</sub></td><td>150</td><td>80</td><td>80</td><td>80</td><td>80</td></tr><tr><td>End redshift z<sub>end</sub></td><td>5.8</td><td>4.0, 0.0, 4.1, 1.0</td><td>0.0</td><td>5.0</td><td>6.7</td></tr><tr><td>Density threshold &delta;<sub>*</sub></td><td>50 &lang;&rho;&rang;</td><td>Star formation</td><td></td><td></td><td></td></tr><tr><td>Temperature threshold T<sub>*</sub></td><td>off</td><td>20 &lang;&rho;&rang;</td><td>30 &lang;&rho;&rang;</td><td>30 &lang;&rho;&rang;, 50 &lang;&rho;&rang;</td><td>200 &lang;&rho;&rang;</td></tr><tr><td>Efficiency &epsilon;<sub>*</sub></td><td>0.02</td><td>2 &times; 10<sup>4</sup> K</td><td>2 &times; 10<sup>4</sup> K</td><td>2 &times; 10<sup>4</sup> K</td><td>2 &times; 10<sup>4</sup> K</td></tr><tr><td>Mass fraction &eta;<sub>SN</sub></td><td>10 per cent</td><td>Supernova feedback</td><td>10 per cent</td><td>10 per cent, 15 per cent</td><td>10 per cent, 15 per cent</td></tr><tr><td>Stellar particle escape fraction f<sub>esc,*</sub></td><td>0.42</td><td>Radiation</td><td>0.25, 0.3, 0.25, 0.3</td><td>0.25</td><td>0.1, 0.14, 0.2, 0.25, 0.3, 0.5</td></tr></tbody></table>

They also suggest that the bulk flow contribution would be dominated by the tSZ effect. According to their estimations, for a CDM cosmology with $\Omega_{0} = 1$, the quadratic Doppler effect does not yield a measurable average (isotropic) Compton-y parameter. Below we evaluate this second-order Doppler contribution based on our simulations for the currently favoured cosmological model.

## 3 SIMULATIONS

All simulations used in this study were performed using the code RAMSES-CUDATON (Teyssier 2002; Aubert & Teyssier 2008), which is a fully coupled, fixed-grid, hybrid (using Graphics Processing Units, GPU, and Central Processing Units, CPU) code combining N-body dynamics and gas dynamics (ran on the CPUs), and (GPU-based) radiative transfer for simulating large-scale structure and galaxy formation. The N-body dynamics solve for the velocities and positions of collisionless dark matter particles. The gas is treated by solving the Euler equations using a second-order Godunov shock-capturing method (Toro 2009). The parameters and set-up of our simulations are summarized in Tables 1 and 2.

### 3.1 Cosmic Dawn II

The Cosmic Dawn (CoDa) II (Ocvirk et al. 2020) is the largest simulation we use here, which serves as our fiducial case. It has comoving box of 64 h$^{-1}$ Mpc on a side, with 4096$^3$ grid for the gas and radiation dynamics, and 4096$^3$ dark matter particles of mass $4.07 \times 10^5$ M$_{\odot}$, with the gravitational forces solved on a 4096$^3$ grid. The simulation uses the Planck 2014 cosmology (Planck Collaboration XVI 2014): $\Omega_{\Lambda} = 0.693$, $\Omega_{m} = 0.307$, $\Omega_{b} = 0.045$, $H_{0} = 67.77$, and power spectrum normalization $\sigma_{8} = 0.8288$ and slope $n = 0.963$. The simulation starts at redshift $z = 150$ and ends at $z = 5.80$.

For storage efficiency, the data for the full resolution grid was reduced to a coarser grid of 2048$^3$ cells. These full-box, lower

<table><caption>**Table 2.** Summary of our simulation set and the input astrophysical parameters set.</caption><thead><tr><th>Simulation</th><th>δ<sub>*</sub></th><th>f</th><th>ε<sub>*</sub></th><th>η</th></tr></thead><tbody><tr><td>CoDa II</td><td>50</td><td>0.42</td><td>0.02</td><td>0.1</td></tr><tr><td>100-δ20-f0.25-ε3.5</td><td>20</td><td>0.25</td><td>0.035</td><td>0.1</td></tr><tr><td>100-δ20-f0.3-ε5</td><td>20</td><td>0.3</td><td>0.05</td><td>0.1</td></tr><tr><td>100-δ20-f0.25-ε6</td><td>20</td><td>0.25</td><td>0.06</td><td>0.1</td></tr><tr><td>100-δ20-f0.3-ε8</td><td>20</td><td>0.3</td><td>0.08</td><td>0.1</td></tr><tr><td>50-δ30-f0.25-ε2.5</td><td>30</td><td>0.25</td><td>0.025</td><td>0.1</td></tr><tr><td>10-δ200-f0.2-ε3.5</td><td>200</td><td>0.2</td><td>0.035</td><td>0.1</td></tr><tr><td>10-δ200-f0.25-ε3.5-η15</td><td>200</td><td>0.25</td><td>0.035</td><td>0.15</td></tr><tr><td>25-δ30-f0.1-ε4</td><td>30</td><td>0.1</td><td>0.04</td><td>0.1</td></tr><tr><td>25-δ30-f0.14-ε4</td><td>30</td><td>0.14</td><td>0.04</td><td>0.1</td></tr><tr><td>25-δ30-f0.2-ε4</td><td>30</td><td>0.2</td><td>0.04</td><td>0.1</td></tr><tr><td>25-δ30-f0.2-ε2</td><td>30</td><td>0.2</td><td>0.02</td><td>0.1</td></tr><tr><td>25-δ30-f0.2-ε2-η15</td><td>30</td><td>0.2</td><td>0.02</td><td>0.15</td></tr><tr><td>25-δ30-f0.25-ε2.5</td><td>30</td><td>0.25</td><td>0.025</td><td>0.1</td></tr><tr><td>25-δ50-f0.3-ε3</td><td>50</td><td>0.3</td><td>0.03</td><td>0.1</td></tr><tr><td>25-δ50-f0.5-ε1</td><td>50</td><td>0.5</td><td>0.01</td><td>0.1</td></tr></tbody></table>

resolution data of all the gas properties, ionizing flux density, and dark matter density field were kept for all snapshots. Furthermore, a catalogue of the stellar particles and halo catalogues are available. For a more detailed description of CoDa II, we direct the reader to Ocvirk et al. (2020).

### 3.2 10, 25, 50, and 100 h$^{-1}$ Mpc simulations

Our suite of auxiliary simulations use cosmology parameters consistent with the latest constraints from the Planck survey (Planck Collaboration VI 2018): $\Omega_{\Lambda} = 0.682$, $\Omega_{m} = 0.318$, $\Omega_{b} = 0.045$, $H_{0} = 67.1$, and power spectrum normalization $\sigma_{8} = 0.833$ and slope $n = 0.9611$. The 15 simulations used here vary in volume, spatial resolution, and star formation parameters. There are four runs with sides $100 h^{-1}$ Mpc (149.0 Mpc), one run with side $50 h^{-1}$ Mpc (74.5 Mpc).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Calibration of our simulations. For clarity not all simulations are shown as separate lines, but instead the 25 *h*⁻¹Mpc box models are shown as bands, the first band (left one on top panels, top one on bottom panels) indicates the range for early reionization cases (25-δ30-*f*0.14-ε4, 25-δ30-*f*0.2-ε4, 25-δ30-*f*0.25-ε2.5, and 25-δ50-*f*0.3-ε3), while thesecond band (right one on top panels, bottom one in bottom panels) includes the late-reionization ones (25-δ30-*f*0.1-ε4, 25-δ30-*f*0.2-ε2, 25-δ30-*f*0.2-ε2-η15, and 25-δ50-*f*0.5-ε1). The respective observational data points are shown as symbols, as indicated on each panel. We show the evolution of the following globally averaged quantities: (top left) neutral hydrogen fraction (Fan, Carilli & Keating 2006; Ouchi et al. 2010; Greig & Mesinger 2017; Davies et al. 2018; Mason et al. 2018; Hoag et al. 2019); (top right) ionized hydrogen fraction (Greig & Mesinger 2017; Davies et al. 2018; Mason et al. 2018; Hoag et al. 2019).(Middle left) Volume-weighted photoionization rate (Faucher-Giguère et al. 2008; Calverley et al. 2011; Wyithe & Bolton 2011; Becker & Bolton 2013; D'Aloisio et al. 2018). (Middle right) Thomson scattering optical depth due to reionization. The measurement from the Planck Collaboration (Planck Collaboration VI 2018) including its 1σ error (horizontal shaded area) is shown for comparison.</figcaption></figure>

Mpc), and two runs with sides 10 *h*⁻¹ Mpc (14.9 Mpc), all with 1024³ cells and 1024³ dark matter particles. We have also performed a further eight runs with sides 25 *h*⁻¹ Mpc (37.25 Mpc) with 512³ cells and 512³ dark matter particles. All simulations adopted initial power spectrum of density fluctuations based on CAMB code,<fnref n="4" /> except for the 25 *h*⁻¹ volumes which used the Eisenstein & Hu (1999) power spectrum. We do not expect this difference to have a notable effect on our results.

We label our simulations by listing the main parameters being varied, as follows Boxsize-δNum1-fNum2-εNum3(-ηNum4), where Boxsize is the size of the simulation volume in *h*⁻¹ Mpc, Num1 is the value of the density threshold for star formation used, Num2 is the escape fraction per stellar particle, Num3 is the star formation efficiency in percents, and Num4 is the supernova feedback load factor (optional, indicated only if different from the fiducial value of 0.1), for example 25-δ30-*f*0.2-ε2-η15.

The main differences between these simulations and CoDa II are the feedback strength, ionizing photon escape fractions and

star formation parameters. For instance, the cooling floor for the interstellar gas is switched off in CoDa II, i.e. there is no temperature limit below which the cells need to be for stars to form. However, the temperature floor is turned on in the auxiliary simulations. The density threshold for star formation was varied between 200 above mean for the 10 Mpc h⁻¹ boxes to 20 for the 100 Mpc h⁻¹ ones, in order to account for the different grid resolution of these volumes. With the exception of 10-δ200-*f*0.25-ε3.5-η15, all simulations run from *z* = 80 to the end of the EoR or later. Simulation 10-δ200-*f*0.25-ε3.5-η15 was stopped early since its ionization fraction was evolving much faster than expected and reionization was due to be completed too early and sooner than observations suggest (see Fig. 1).

## 3.3 Calibration of the simulations

In Fig. 1, the top 2 panels show the evolution of the mean neutral and ionized hydrogen fractions, respectively, along with the corresponding observational constraints. The observational data suggest that reionization is largely complete by *z* ~ 6, at which point the mean neutral fraction sharply drops to ⟨xH₁⟩ ~ 10⁻<fnref n="4" /> as the ionized

<footnote n="4">⁴https://camb.readthedocs.io/en/latest/</footnote>
patches overlap. The universe becomes largely transparent, and the neutral fraction continues to decrease, but more slowly. This sharp transition marks the end of reionization. Our simulations follow the same general trend.

For clarity, we have grouped the reionization histories of the four 25 Mpc boxes that reionize earlier in the rose-brown band, while the blue band represents the four 25 Mpc simulations that reionize later. The two high-resolution, small-volume simulations (dashed lines) were not run all the way to completion of reionization, but only to $z \sim 6.3$ ($10-\delta200-f0.2-\epsilon3.5$) and $z \sim 6.7$ ($10-\delta200-f0.25-\epsilon3.5-\eta15$), but both follow the same general trend.

The overlap epoch suggested by data is best matched by CoDa II, $100-\delta20-f0.3-\epsilon5$, $100-\delta20-f0.25-\epsilon6$, $50-\delta30-f0.25-\epsilon2.5$, and $25-\delta50-f0.5-\epsilon1$. The additional simulations, which reionize slightly earlier and later than this allow us to probe the effects of the timing of the EoR on the tSZ signal. We note that CoDa II, and to a lesser extent the other early-reionization scenarios tend to overionize after overlap, while the late-reionization models tend to agree with post-reionization data better. This does not have significant effect on the tSZ signals, however.

### 3.3.1 Global photoionization rate

In Fig. 1, bottom left panel we show the evolution of the mean global photoionization rates. The evolution of the mean photoionization rate roughly follows the same pattern for all simulations, mirroring the reionization histories. The early-reionization cases have higher photoionization rates, since their background ultraviolet (UV) flux densities are higher. More ionizing photons in the IGM allows for reionization to proceed more rapidly in these runs. Conversely, the simulations that reionize later have lower rates of photoionization. Around $z \sim 6$, there is a noticeable upturn in the slopes. This corresponds to the end of reionization, where there is a drop in ionizing photon absorption as the neutral hydrogen in the boxes are ionized, increasing the background UV flux. While the slopes match the pattern of the observational data, most of the simulations overshoot the observations post-reionization. This discrepancy is responsible for those simulations having neutral fractions lower than those observed. However, the blue band of 25 Mpc simulations which best match with the neutral fraction observations also agree with the observed data for the photoionization rates. We also note that the photoionization rate evolution in the higher-resolution simulations (CoDa II and the two $10$ Mpc $h^{-1}$ per side volumes) is somewhat different, with initially low values, but a steeper rise around overlap.

The reason for the observed fluctuations of the photoionization rates are early times is a combination of suppression of sources due to radiative feedback and the small number statistics of sources at early times. Since early-on there are relatively few, low-mass and highly clustered sources, they tend to be easily suppressed by photoionization feedback, and their low numbers mean that this is seen as fluctuations in the mean global photoionization rate. Later on there are many more, and on average larger sources, and thus the mean photoionization rate becomes smooth.

### 3.3.2 Optical depth to reionization

Another key observational constraint is the integrated Thomson scattering optical depth to recombination, showing in Fig. 1, bottom left panel, along with the corresponding result from Planck (Planck Collaboration VI 2018) and its $1\sigma$ errors (grey shared area). The helium is not explicitly accounted for in our simulations, thus its

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Evolution of the global star formation rate density. The dust-corrected and dust-uncorrected observations from Bouwens et al. (2015) are indicated by the grey shaded region. The bottom band includes the results from the four 100 Mpc boxes, the top band indicates the range for early-reionization $25~h^{-1}$Mpc box models, and the middle band includes the late-reionization ones.</figcaption></figure>

contribution is included as discussed in (Iliev et al. 2006). Assuming that for $z > 3$ He is singly ionized wherever the ionized fraction of hydrogen is above 0.95, and for $z < 3$ helium is always double-ionized.

All cases are in rough agreement with the Planck results, at the low end, or just outside the $1\sigma$ uncertainty. The late-reionization scenarios naturally have lower optical depths, as there is less free electrons in the IGM along the line of sight between the end of the EoR and the present day; and vice versa. The simulations whose optical depths best agree with observations are $100-\delta20-f0.3-\epsilon8$, $25-\delta30-f0.2-\epsilon4$, and $10-\delta200-f0.25-\epsilon3.5-\eta15$.

### 3.3.3 Star formation rate

Finally, in Fig. 2, we show the evolution of the global star formation rate density. The grey shaded area shows the range between observed dust-corrected and dust-uncorrected star formation rate (SFR) densities from (Bouwens et al. 2015). Here, the rose-brown band represents the set of higher star formation 25 Mpc volumes ($25-\delta30-f0.1-\epsilon4$, $25-\delta30-f0.14-\epsilon4$, $25-\delta30-f0.2-\epsilon4$), while the blue band represents the lower star formation 25 Mpc volumes ($25-\delta30-f0.2-\epsilon2$, $25-\delta30-f0.2-\epsilon2-\eta15$, $25-\delta30-f0.25-\epsilon2.5$, $25-\delta50-f0.3-\epsilon3$), and the yellow band represents the 100 Mpc volumes. Similarly to the observational data, the general trend of the SFR densities for all simulations is a continuous increase, with roughly the same slopes. All simulations are in rough agreement with the SFR data, particularly most 25 Mpc volumes, the 50 Mpc one, and CoDa II. The timing of reionization has little correlation with the SFR.

The variation in the SFR of the simulations is mainly governed by the parameters $\delta_*$ and $\epsilon_*$.$^1$ Since $\delta_*$ indicates the minimum gas overdensity required to trigger star formation, its value is tied to the spatial resolution of the box. Setting $\delta_*$ higher results in fewer stars being formed since the (fixed-grid) simulation is less able to resolve the high-density gas clouds. On the other hand, low $\epsilon_*$ values mean that smaller fraction of the gas is converted into stars, resulting in lower SFR. Hence, the balancing act between these two parameters is needed for obtaining a SFR corresponding to the observational data. Unlike the previous plots, the SFR by itself does not directly indicate the progression of the EoR, i.e. the relationship between them is not as straightforward. The amount of ionizing photons emitted into
the IGM by the stellar population is dictated by the value of $f_{esc,*}$. This parameter can be adjusted to ensure that the timing of the EoR coincides with observations. Thus, it is possible for a simulation to have a low SFR but still complete reionization at a reasonable time. For example 25-δ50-f0.5-ε1 has one of the lowest SFRs, due to its high δ and low ε*, but because its assumed $f_{esc,*}$ is quite high, a greater proportion of photons are emitted into the IGM per stellar particle, increasing the photoionization rate, and allowing reionization to end at $z \sim 6$. Another example of how the SFR history alone does not govern the timing of the EoR is the fact that the simulations which have similar SFR histories (grouped together in the three coloured bands) end reionization at different times. Within each coloured band, there are simulations which finish reionising as early as $z \sim 6.5$ and those which complete reionization as late as $z \sim 5.5$. Again, this is a direct result of the different $\delta_*$, $\epsilon_*$, and $f_{esc,*}$ parameters used for them. We see that increasing $\eta_{SN}$ by 5 per cent has only a small influence on the EoR, and a negligible change in the SFR. 25-δ30-f0.2-ε2 completes reionization slightly sooner ($z \sim 5.8$) than does 25-δ30-f0.2-ε2-η15 ($z \sim 5.7$), and their SFR density histories are almost identical. Hence, this small increase in supernova feedback only minimally influences the reionization process. We also note the high SFR of CoDa II. In addition to the different values used for $\delta_*$, $\epsilon_*$, and $f_{esc,*}$, CoDa II has different simulation settings compared to the others, particularly the switching off of the temperature threshold for star formation. By turning off this setting, all cells, regardless of their temperature, are eligible for star formation if they are above the density threshold. This means that the conditions for star formation are solely based on $\delta_*$ and $\epsilon_*$. This results in there being more regions where star formation is possible. While radiative suppression of star formation still occurs in CoDa II, it is not as intense as it would be if the temperature threshold was turned on.

# 4 METHOD

## 4.1 Electron pressure lightcones

The input data for our calculations are the full-box, reduced-resolution cubes of the volume-weighted gas pressure, and mass-weighted ionization fraction (CoDa II) and the full-resolution cubes of the gas pressure and ionization fraction (the rest of our simulations). These fields are then used to calculate the Compton y-parameter (equation 6). We start by constructing lightcones of the electron pressure, which for pure hydrogen gas is given by

$$p_e = \frac{x_{HII}}{1 + x_{HII}} p_{\text{gas}}.$$

We construct the lightcone as follows. Consider a photon travelling along a line of sight, cell by cell, from initial redshift, $z_i$, to a final one, $z_f$. The age of the universe, initial time, and redshift can be readily calculated (see e.g. equation 30 in Hogg 1999). We then convert the cell size from comoving to physical length by using the corresponding scale factor $a$, and thus calculate the time it will take the photon to cross the cell, $dt$. The age of the universe after this cell crossing is now $t = t_1 + dt$. Assuming an Einstein-de Sitter universe, which is a good approximation at high redshift, the photon is now at redshift $z = (1 + z_i)(t/t_i)^{\frac{2}{3}} - 1$. We can now calculate the electron pressure of the cell at $z$ by interpolating the data at the snapshots on either side of this redshift, say $z_1$ and $z_2$ corresponding to snapshots 1 and 2, respectively.

The electron pressure is interpolated using a Sigmoid function

$$g = \frac{1}{1 + e^{-\beta z_p}}, \quad (9)$$

with $\beta = 2$ and $z_p = -10 + 20 \left( \frac{z-z_2}{z_1-z} \right)$, as follows

$$p_{ei} = (1 - g)p_{e1} + g p_{e2}, \quad (10)$$

where $p_e$ is the interpolated electron pressure at $z$, and $p_{e1}$ and $p_{e2}$ are the electron pressure of the cell at $z_1$ and $z_2$ (snapshots 1 and 2), respectively. We then repeat the process, calculating the cell size, the time taken for the photon to cross it, and so on at the updated redshift $z$. This allows us to obtain a lightcone of the electron pressure: a 3D grid of the field where each slice in the direction of light propagation reflects the state of the universe at that moment in time.

We then perform an integration along the light travel path in order to obtain the Comptonization parameter $y$. For each box light crossing, we numerically integrate the pressure over all lines of sight along each axis using the composite Simpson's rule. For the integral we use the previously calculated time intervals between each cell in the lightcone ($dt$). In order to avoid artificial amplification of the y-parameter due to structural repetition during the box replication after each light crossing, the volume was randomly shifted and rotated during the integration stage. This procedure yields a 2D grid of the y-parameter.

We create lightcones for the range of redshift for each simulation extending from the redshift of overlap, $z \sim 6$, to 12, corresponding to the duration of the bulk of EoR. An exception to this are the 10 cMpc boxes, which were not run all the way to redshift 6, and thus we constructed lightcones between $z \sim 12$ and their last available snapshots, ($z \sim 6.38$ for $10-\delta200-f0.2-\epsilon3.5$ and $z \sim 6.74$ for $10-\delta200-f0.25-\epsilon3.5-\eta15$). The number of light crossings required to construct the lightcones vary from $\sim 10.8$ (100 cMpc boxes), to 17.2 light crossings (CoDa II), and up to $\sim 88.4$ light crossings for $10-\delta200-f0.25-\epsilon3.5-\eta15$ and $\sim 97.8$ for $10-\delta200-f0.2-\epsilon3.5$.

We present the results in the form of maps of the y-parameter, probability density functions of $y$, and angular power spectra of the signal. The maps are simply the images of the 2D y-parameter grids obtained from the integration. We consider the full resolution results as well as with smoothed maps with Gaussian beams with full width at half-maximum (FWHM) of 1.2 arcmin and 1.7 arcmin. These beams correspond to the resolution of the South Pole Telescope (SPT)<fnref n="5" /> (Carlstrom et al. 2011) at 150 GHz and 95 GHz, respectively. We do not use the SPT 220 GHz beam as the tSZ effect disappears near this frequency. We smooth the data by converting the angular resolutions of the beams to comoving Mpc at the redshift of the lightcone at the end of each light crossing. We then calculated the number of cells spanning this size and use this value as our FWHM for the Gaussian beam. Since angular sizes vary negligibly at the redshifts considered here, we ignore this effect when constructing the lightcones.

We also consider the contribution to the y-parameter for each redshift of the lightcone in order to determine the period of dominant contribution. For these PDFs, we calculate the distribution of the y-values for the start of the lightcone to each redshift, i.e. from $z \sim 12$ to 11, $z \sim 12$ to 10, etc.

For comparison with previous results and data we compute the angular power spectra of the maps and show them along with the Shaw et al. (2010) power spectrum template scaled to fit the cosmological parameters as done by George et al. (2015), and the

<footnote n="5">⁵https://pole.uchicago.edu/</footnote>
SPT observations of the total CMB power spectrum at the 95 GHz and 150 GHz bandpowers.

## 4.2 Separate contribution lightcones

We also investigate the separate contributions to the y-parameter and its fluctuations from each constituent field-gas density, temperature, and ionization fraction. We do this as follows, using the volume-weighted gas density and mass-weighted ionization fraction.

The fluctuations in the temperature field were removed by using a fixed global value for the temperature, namely $T_e = 30\,000$ K. We then constructed an electron density lightcone via the same steps as for the electron pressure lightcones in Section 4.1 for the redshift range $z \sim 6$ to 12, and using the expression

$$ \rho_e = \frac{x_{HII}}{1 + x_{HII}} \rho_{gas} \quad (11) $$

for the electron density. After the lightcone was constructed, we numerically integrated the electron density over all lines of sight along each axis and along the direction of light propagation, as before. The lightcone integral then becomes

$$ y = \frac{\sigma_T k_B}{m_e m_p c} T_e \int \frac{x_{HII}}{1 + x_{HII}} \rho_{gas} dt, \quad (12) $$

where $m_p$ is the proton mass and comes from converting the electron mass density to number density ($n_e = n_{HII} = \frac{\rho_{HII}}{m_p}$). We also randomly shifted and rotated the box to avoid artificially boosting the signal. Again, this yields a 2D grid of y-parameter values for our volume.

Next, we removed the patchiness of the EoR in addition to the temperature fluctuations, by using a globally averaged value of the ionization fraction for each light crossing. We made a lightcone of gas density using only the volume-weighted gas density files for CoDa II for the interpolation process, for the redshift range $z \sim 6$ to 12. In this case, equation (12) becomes

$$ y = \frac{\sigma_T k_B}{m_e m_p c} T_e \langle x_{HII} \rangle \int \rho_{gas} dt, \quad (13) $$

where $\langle x_{HII} \rangle$ is the globally averaged ionization fraction. We used the value of $\langle x_{HII} \rangle$ corresponding to the redshift of the box at the end of each light crossing. Once again, we randomly shifted and rotated the box during integration, obtaining a 2D map of the y-parameter.

### 4.2.1 Phase diagrams

In order to further probe the underlying quantities contributing to the y-parameter, we also construct and compare the phase diagrams of the gas temperature and density. Since the electron pressure is dependent on these quantities, this allows for better understanding of the trends and dependence on different parameters. We present these phase diagrams and an interpretation of them with respect to the y-parameter values in Section 5.3.

## 4.3 Additional corrections and tests: y-parameter for uniform IGM

If we assume instantaneous reionization, a constant post-reionization IGM temperature, and uniform density at the mean value for the universe, then the y-parameter mean value can be calculated analytically. While clearly unrealistic, this yields an useful reference value for our more detailed results later.

We must first rewrite the integral with respect to redshift, since we are concerned with the redshift evolution of the IGM density, The Friedmann equation:

$$ H^2 = H_0^2 [\Omega_{r,0} (1+z)^4 + \Omega_{m,0} (1+z)^3 + \Omega_{k,0} (1+z)^2 + \Omega_{\Lambda,0}] . \quad (14) $$

for a flat universe, $\Omega_{k,0} = 0$, and at high redshift, $\Omega_{r,0} \sim 0$ and $\Omega_{m,0} \gg \Omega_{\Lambda,0}$ reduces to

$$ \left(\frac{\dot{a}}{a}\right)^2 = H_0^2 [\Omega_{m,0} (1+z)^3]. \quad (15) $$

Then, equation (6) yields

$$ y = - \frac{\sigma_T}{m_e c^2} \frac{1}{H_0 \Omega_{m,0}^{\frac{1}{2}}} \int p_e (1+z)^{-\frac{5}{2}} dz. \quad (16) $$

Assuming the IGM is an ideal gas, $p_e = k_B n_e T_e$, where $n_e$ and $T_e$ are the electron number density and temperature, respectively. Also, for pure-hydrogen gas we have

$$ n_e = n_{HII} = \frac{\rho_{HII}}{m_{HII}} \quad (17) $$

and

$$ \rho_{HII} = \rho_{gas} \frac{x_{HII}}{1 + x_{HII}}. \quad (18) $$

Since $m_{HII}$ is simply the proton mass, we shall use the notation $m_p$ henceforth. For instantaneous reionisation, $x_{HII}$ can only be either 0 or 1 if neutral or ionized, respectively. Hence, $\frac{x_{HII}}{1+x_{HII}} = \frac{1}{2}$ for a reionized IGM. Thus the electron number density for uniform IGM is

$$ n_e = \frac{1}{2} \frac{\rho_{gas}}{m_p} = \frac{1}{2} \frac{\rho_0 \Omega_{b,0} (1+z)^3}{m_p}, \quad (19) $$

where $\rho_0 = 10^{-30} \text{ g cm}^{-3}$ is the mean density of the universe at $z=0$. Converting the integrand from electron pressure to density for an ideal gas, we get

$$ y = - \frac{1}{2} \frac{\sigma_T k_B}{m_e m_p c} \frac{\rho_0 \Omega_{b,0}}{H_0 \Omega_{m,0}^{\frac{1}{2}}} T_e \int (1+z)^{\frac{1}{2}} dz. \quad (20) $$

Evaluating this integral for the redshift range of the lightcones ($z=12$ to 6) and the cosmological parameters for CoDa II, and assuming an average IGM temperature of 30 000 K, we obtain a Comptonization parameter $y \approx 4.22 \times 10^{-8}$. For the cosmology of the other simulations, we find a very similar mean value of $y \approx 4.18 \times 10^{-8}$.

## 4.4 Helium reionization

So far, we have assumed that the baryonic matter is comprised of hydrogen only, ignoring helium and metals. While the presence of metals is insignificant to this study, the presence of helium will contribute to the tSZ signal. Neutral helium (He) requires photons with energy of at least 24.6 eV to be singly ionized to He II, with He II recombining at roughly the same rate as H II. However, He II requires photons of at least 54.4 eV energy to fully ionize He III, with a recombination rate over 5 times that of hydrogen. Thus, when accounting for the presence of helium, the sources that ionize hydrogen are often considered to singly ionize helium. On the other hand, the second ionization of helium occurs at later redshifts, with the reionization of helium ending at $z \sim 3$ (e.g. Barkana & Loeb 2001). Therefore, during the redshift range considered in this study, we assume that each helium atom contributes one electron to the
IGM whenever hydrogen is ionized. Equation (17) thus becomes

$$n_e = n_{\text{HII}} + n_{\text{HeII}}, \quad (21)$$

where $n_{\text{HeII}}$ is the number density of singly ionized helium. When accounting for doubly ionized helium, the term $2n_{\text{HeII}}$ is added.

Assuming the primordial abundances of hydrogen and helium by mass ($X = 0.76$ and $Y = 0.24$, respectively), and ignoring any isotopes, we get

$$n_e = x_{\text{HII}}X \frac{\rho_{\text{gas}}}{m_H} + x_{\text{HeII}}Y \frac{\rho_{\text{gas}}}{4m_H} = \left(x_{\text{HII}}X + x_{\text{HeII}}\frac{Y}{4}\right) \frac{\rho_{\text{gas}}}{m_H}, \quad (22)$$

where $x_{\text{HeII}} = n_{\text{HeII}}/n_{\text{He}}$ is the ionization fraction of singly ionized helium. For instantaneous reionization, we have $x_{\text{HII}}$ and $x_{\text{HeII}}$ jump from 0 to 1 at the transition, and we have

$$n_e = 0.82 \frac{\rho_{\text{gas}}}{m_H}. \quad (23)$$

This results in an 18 per cent decrease in the tSZ signal, compared to the hydrogen-only numbers. The estimated mean Comptonization parameters for CoDa II and the auxiliary simulations therefore become $y \approx 3.46 \times 10^{-8}$ and $y \approx 3.43 \times 10^{-8}$, respectively. In addition, RAMSES-CUDATON only tracks the ionization and cooling processes for atomic hydrogen. In doing so, it assumes that hydrogen makes up 76 per cent of the baryonic matter, thus preserving its primordial abundance. Nevertheless, the temperature field generated by the simulation is given in units of $K\mu^{-1}$, where $\mu$ is the mean molecular weight of the gas. Since $\mu$ depends on the ionization fractions,

$$\frac{1}{\mu} = (1 + x_{\text{HII}})X + (1 + x_{\text{HeII}} + 2x_{\text{HeII}})\frac{Y}{4}, \quad (24)$$

the temperature field (and, hence, the pressure field) can be updated, in post-processing, to roughly account for helium reionization, by interpolating $\mu$ values for the required redshift range.

## 4.5 Quadratic (second-order) Doppler contribution

We estimate the contribution to the y-parameter due to the second-order Doppler distortions using a similar methodology to the one used for calculating the mean Compton-y from tSZ effect. We start by isolating the $O(v^2)$ term in equation (8) and recast the integral to be dependent on $z$

$$y = \frac{\sigma_T}{3cH_0} \int_{z_0}^{z_1} \frac{n_e \langle v^2 \rangle}{(1+z)E(z)} dz, \quad (25)$$

where $n_e$ and $\langle v^2 \rangle$ are extracted from the simulations, as follows. For each cell, we first compute

$$n_e = \frac{\rho_{\text{HII}}}{m_{\text{HII}}} = \frac{X \rho_x x_{\text{HII}}}{m_p} \quad (26)$$

and

$$v^2 = v_x^2 + v_y^2 + v_z^2, \quad (27)$$

where we have dropped the angled brackets for clarity, before taking the average over a snapshot to obtain $\langle n_e v^2 \rangle$ for a given $z$. With these values of $\langle n_e v^2 \rangle$ as a function of $z$ in hand, we numerically integrate equation (25) between $z = 6$ and $z = 12$, linearly interpolating between the available snapshots. We confirm that the result is converged with respect to the interpolation sampling.

The power of the velocity field peaks at fairly large scales ~ 100 Mpc (for more detailed discussion see e.g. Iliev et al. 2007), thus a significant fraction of that power is missing in our relatively small simulation volumes. However, at large scales the velocity field

<table><caption>**Table 3.** The mean Compton y-parameter for each simulation.</caption><thead><tr><th>Simulation</th><th>⟨y⟩ × 10<sup>8</sup></th></tr></thead><tbody><tr><td>CoDa II</td><td>3.67</td></tr><tr><td>100-δ20-f0.25-ε3.5</td><td>1.22</td></tr><tr><td>100-δ20-f0.3-ε5</td><td>1.96</td></tr><tr><td>100-δ20-f0.25-ε6</td><td>1.91</td></tr><tr><td>100-δ20-f0.3-ε8</td><td>2.70</td></tr><tr><td>50-δ30-f0.25-ε2.5</td><td>2.47</td></tr><tr><td>10-δ200-f0.2-ε3.5</td><td>2.62</td></tr><tr><td>10-δ200-f0.25-ε3.5-η15</td><td>2.46</td></tr><tr><td>25-δ30-f0.1-ε4</td><td>2.86</td></tr><tr><td>25-δ30-f0.14-ε4</td><td>2.18</td></tr><tr><td>25-δ30-f0.2-ε4</td><td>3.63</td></tr><tr><td>25-δ30-f0.2-ε2</td><td>1.97</td></tr><tr><td>25-δ30-f0.2-ε2-η15</td><td>2.03</td></tr><tr><td>25-δ30-f0.25-ε2.5</td><td>2.96</td></tr><tr><td>25-δ50-f0.3-ε3</td><td>3.09</td></tr><tr><td>25-δ50-f0.5-ε1</td><td>1.82</td></tr></tbody></table>

is linear and thus the effect of this missing power can be readily estimated and included in our calculations, as follows. Starting from the power spectrum of density fluctuations, defined as

$$\langle \delta(\mathbf{k}, z) \delta(\mathbf{k}', z) \rangle = (2\pi)^3 \delta_D(\mathbf{k} + \mathbf{k}') P(\mathbf{k}, z), \quad (28)$$

where $\delta(\mathbf{k}, z)$ is the Fourier transform of the real-space density fluctuations and $\delta_D$ is the Dirac delta function, we use the continuity relation

$$v(\mathbf{k}, z) \approx \frac{iH(z)f[\Omega_m(z)]k}{(1+z)k^2} \delta(\mathbf{k}, z), \quad (29)$$

where $f[\Omega_m(z)] \approx \Omega_m(z)^{0.6}$ (Lahav et al. 1991), to write the total mean-square fluctuations for the velocity in linear theory as

$$\langle v^2 \rangle_{\text{tot}} = \frac{H(z)^2 f[\Omega_m(z)]^2}{2\pi^2(1+z)^2} \int_0^\infty P(k)dk. \quad (30)$$

Finally, we estimate the mean-square velocity missing from a simulation box through

$$\langle v^2 \rangle_{\text{missing}} = \langle v^2 \rangle_{\text{tot}} - \langle v^2 \rangle_{\text{box}}, \quad (31)$$

where $\langle v^2 \rangle_{\text{box}}$ is the mean-square velocity calculated directly from the simulation box. This missing power can then be added as a correction to the simulation results from sub-box scales.

# 5 RESULTS

## 5.1 Electron pressure lightcones

The main results of our study is the tSZ effect calculated from the electron pressure lightcones calculated as discussed in Section 4.1. The mean y-parameter values, $\langle y \rangle$, found from these lightcones for the 16 simulations are listed in Table 3. All values of $\langle y \rangle$ we find are of order a few × $10^{-8}$, well below the COBE-FIRAS limit of $\langle y \rangle < 1.5 \times 10^{-5}$ (Fixsen et al. 1996), and 2 orders of magnitude lower than the total mean y-parameter estimated by Refregier et al. (2000), Zhang et al. (2004), and Hill et al. (2015). This is expected, since galaxy clusters provide the dominant contribution to this quantity, at $\langle y \rangle_{\text{ICM}} = 1.58 \times 10^{-6}$ (Hill et al. 2015). Compared with the reionization contribution estimated by Hill et al. (2015), namely $\langle y \rangle_{\text{EoR}} = 9.8 \times 10^{-8}$), our values are slightly lower, but of the same order of magnitude. This discrepancy is most likely due to differences between our simulations and their reionization model. Nevertheless, these values for the mean Comptonization parameter for the EoR are comparable, albeit approximate models
<figure><img src="image_6.png" /><figcaption>**Figure 3.** PDF distributions of the Comptonization y-parameter for our simulations, some of which are shown individually (as labelled) and for clarity the rest are grouped in two bands – turquoise one including 100-δ20- *f*0.3-ε5, 100-δ20- *f*0.25-ε6, 25-δ30-*f*0.1-ε4, 25-δ30-*f*0.2-ε2, 25-δ30-*f*0.2-ε2-η15, and 25-δ50-*f*0.5-ε1, and coral one for 25-δ30-*f*0.14-ε4, 25-δ30-*f*0.25-ε2.5, and 25-δ50-*f*0.3-ε3.</figcaption></figure>

might sometimes overestimate the mean EoR signal. As we will see below, the value of the y-parameter closely correlates with the end of reionization redshift, with early overlap yielding higher y-values.

### 5.1.1 The tSZ signal and the timing of the EoR

For further insight and understanding of these values, we consider the full probability density functions (see Fig. 3). For clarity, we only show a few of the simulations individually (as labelled), and we group the rest into two bands – first includes 100-δ20-*f*0.3-ε5, 100-δ20-*f*0.25-ε6, 25-δ30-*f*0.1-ε4, 25-δ30-*f*0.2-ε2, 25-δ30-*f*0.2-ε2-η15, and 25-δ50-*f*0.5-ε1 (turquoise) and second – 25-δ30-*f*0.14-ε4, 25-δ30-*f*0.25-ε2.5, and 25-δ50-*f*0.3-ε3 (coral). Simulations in the former group have lower y-parameter values (y ≲ 2 × 10⁻⁸) and generally late end of reionization. The only simulation with an even smaller tSZ signal is 100-δ20-*f*0.25-ε3.5, which is also the case with the latest end of EoR. The second band (coral) groups cases with higher y-values and early-finishing reionization.

At the high-end of y-values are CoDa II and 25-δ30-*f*0.2-ε4, which is the simulation with the earliest overlap and highest y-parameter values among the auxiliary simulations. Their PDFs overlap, save for the long high-y tail in the case of CoDa II, which is a result of its higher resolution, and different supernova feedback model. This allows resolving small, high-temperature regions around supernova explosions, resulting in small areas of very high pressure, and therefore bright in tSZ. However, we note that despite their similar PDFs, the reionization histories of those two simulations differ significantly, with 25-δ30-*f*0.2-ε4 reionizing considerably earlier (z<sub>reion</sub> ~ 6.6) than CoDa II (z<sub>reion</sub> ~ 6.1).

We also note that earlier and faster reionization correlates with narrower y-parameter distributions, since such scenarios mean most of the gas being ionized at similar redshifts, and the gas has more time to equilibrate after overlap, resulting in less fluctuations in the ionized hydrogen fraction and more narrow range in y-values. These results indicate a link between the timing of the EoR and the strength of the tSZ effect stemming from that cosmic period. When reionization ends earlier, there is a longer period of fully ionized gas, thus higher y-parameter overall.

In order to better visualize the relationship between the tSZ effect and the timing of the EoR, in Fig. 4 we show the mean y-parameter

<figure><img src="image_7.png" /><figcaption>**Figure 4.** The mean y-parameter for each simulation (as indicated in the legend) versus z<sub>(xH<sub>I</sub>)</sub>~0.001, the redshift at which x<sub>H<sub>I</sub></sub> ~ 0.001.</figcaption></figure>

versus the redshift at which the average neutral hydrogen fraction reaches 10⁻³, referred to as z<sub>(xH<sub>I</sub>)</sub>~0.001. We picked this late point of the reionization history as it is relatively more accessible by observations compared to earlier stages. We see that there is indeed a clear positive correlation between 〈y〉 and z<sub>(xH<sub>I</sub>)</sub>~0.001, albeit with significant scatter. Curiously, the most significant outliers, namely CoDa II and 10-δ200-*f*0.25-ε3.5-η15, are among the best-resolved simulations in our set. However, we do not believe this is a resolution effect, since these two outliers are at opposite extremes, while our other high-resolution simulation, 10-δ200-*f*0.2-ε3.5, is exactly on the trend, as are the rest of our simulations, despite their varying resolutions. For CoDa II, this is most likely due to the different simulation settings used in that run, while the lightcone for 10-δ200-*f*0.25-ε3.5-η15 has the shortest redshift range, resulting in lower mean y-parameter.

In addition to their different EoR timings, these boxes have varying SFR histories, which also influences the observed scatter. For two simulations with similar y-parameters but different EoR timings, the SFR of the simulation with a later EoR is higher than that of the one whose EoR ends earlier. This suggests that, although a simulation may reionize later, if the corresponding SFR is relatively high, the tSZ effect is boosted. A higher SFR leads to more stars which then go supernova and raise the local gas temperature and pressure. The electron temperature therefore increases, boosting the inverse Compton scattering.

Finally, we consider the differential contributions of different redshift bands to the total, integrated Compton y-parameter. The corresponding PDFs are shown in Fig. 5. Over time, the cumulative PDFs show a clear evolution. It starts as a very wide, double-peaked distribution, with a peak at low values for not yet reionized lines of sight and a broad peak corresponding to the first H II regions. Over time, the signal naturally accumulates, moving towards larger y-values and the low-y peak disappears as more and more ionized patches appear. Less obviously, the distribution grows much narrower over time, sharply peaked around the mean value, and an ever larger high y-values tail appears, due to the incleasing and wider-spread contribution of superheated gas due to supernovae feedback at later times.

### 5.1.2 Impact on galaxy cluster measurements

Another objective of our study is to probe the effect of the EoR tSZ signal on galaxy cluster measurements. Although the mean tSZ signal from the EoR is sub-dominant to that of clusters, this might not be the true at all scales. Furthermore, even if sub-dominant, it
<figure><img src="image_8.png" /><figcaption>**Figure 5.** PDFs of the contributions of different redshift bands to the Comptonization y-parameter for the Cosmic Dawn II simulation.</figcaption></figure>

is useful to quantify how much error the former might contribute toward the latter. We start by showing sample Compton y-parameter maps from several of our simulations in Fig. 6 and the angular power spectra for all simulations in Fig. 7. In Fig. 6, we present the full-resolution map (left panels) as well as maps smoothed with Gaussian beams of FWHM of 1.2 arcmin (middle panel) and 1.7 arcmin (right panel), corresponding to those of the 150 GHz and 95 GHz SPT channels, respectively. The fluctuations in the full-resolution maps are strongest in the sub-arcmin scale, with the largest y-parameters reaching $y \sim 10^{-6}$. When the maps are smoothed, these fluctuations are largest at $\sim 1$ arcmin, with values of order $y \sim \text{few} \times 10^{-8}$. On the other hand, the y-parameter of galaxy clusters is typically $y < 10^{-5}$, and their angular size ranges from from tens of arcsec to tens of arcmin. From the power spectra, we see that, at cluster scales ($l < 10^4$), the EoR accounts for <1 per cent of the cluster signal, but could be significantly larger fraction at smaller scales. These observations are more quantitatively confirmed by the angular power spectra from our full-resolution maps

The angular power spectra of the full-resolution maps shown in Fig. 7 along with the Bolliet et al. (2018) tSZ template (1-halo + 2-halo terms; solid black line) and the total CMB power spectrum measured by the SPT (George et al. 2015). Our predicted sky power spectra are flatter in shape than the Bolliet et al. (2018) template. On smaller scales ($l > 10^4$), the tSZ signals from CoDa II and the 10 Mpc boxes become stronger, possibly surpassing the post-EoR signal. The other simulations have lower resolutions and do not allow for the signal at these small scales to be calculated. For the current instruments like SPT with an angular resolution of $\sim 1$ arcmin, the areas of highest electron pressure in the EoR produce a tSZ signal which is about 3 orders of magnitude smaller than that of clusters. However, a more precise instrument, with arcsec resolution and higher sensitivity, would be able to detect the small regions where $y \sim 10^{-6}$ in our full-resolution maps, and extend the power spectra to $l \sim 10^5$, where the EoR contribution becomes stronger, and possibly dominant. Although the EoR is a sub-dominant contributor to the total tSZ effect on larger scales, it nevertheless should be considered for precision cosmology.

## 5.2 Density lightcones

The tSZ signal from the reionization epoch is formed by the complex interactions of a number of factors and quantities, making it difficult to fully understand. In order to disentangle and evaluate these different contributions to the signal, we next consider them MNRAS 540, 1432–1448 (2025)

separately, as discussed in Section 4.2. Specifically, we compare the full tSZ results to series of simplified cases where we: (1) fix the electron temperature to $T_e = 30,000$ K throughout the volume; (2) fix the electron temperature as well as the ionized fraction to its globally averaged value at that redshift; and (3) the analytic case in Section 4.3, of instantaneous reionization and constant values of all fields. All calculations are done based on the CoDa II data and results are shown in Fig. 8.

The peak of each distribution lies near to the analytic result, with the two patchy reionization cases peaking slightly below it, and the uniform reionization case slightly above it. Furthermore, the patchy reionization probability distribution functions (PDFs) are wider, with significant fraction of the lines of sight having low y-values, below $y = 4.22 \times 10^{-8}$. In contrast, the uniform reionization case yields much narrower distribution, with almost no lines of sight below the analytical result. This is due to the presence of cold, neutral cells in the patchy reionization cases, contributing little to the integrated Compton y parameter.

Conversely, the high-y tail we observed in the full CoDa-II results is primarily due to the highly heated gas from SNe explosions and structure formation shocks, resulting in local temperatures up to millions of K. When the electron temperature fluctuations are removed (blue) this high-y tail completely disappears, leaving an almost Gaussian distribution. When also the reionization patchiness is removed a small tail re-appears at $log_{10}(y) \sim -7$ to $-6.5$ – this is due to the high-density peaks not being fully ionized in this scenario. Patchy reionization proceeds in an inside-out fashion (Iliev et al. 2006), with the denser regions ionized earlier, on average, compared to the low-density ones.

In summary, compared to the base density-fluctuations-only scenario, the ionization fraction fluctuations are responsible for the broad peak of the PDF distribution, while eliminating the modest high-y tail due to the density peaks being ionized first, while the temperature fluctuations yield the high-y tail in the distribution, with some values reaching as high as $log_{10}(y) \sim -6$ to $-5.5$.

## 5.3 Phase diagrams

In this section, we present phase diagrams of the gas temperature versus density (equation of state) for a selected sub-set of our simulations (Figs 9 and 10). These plots reflect the local state of the gas and ultimately determine the range and variety of y-parameters we observed. The distributions are considered at $z \sim 6$, i.e. around, or soon after, overlap.

There are a number of features common for all cases, along with some instructive differences. The majority of cells lie along a diagonal line from $log_{10} T \sim 3.5$ to $4.5$, corresponding to photoionization equilibrium in the underdense regions ($log_{10} \rho < -26.5 \text{ cm}^{-3}$), where gas atomic line cooling is not very strong. Volumes with higher than average density have shorter cooling times, but are also continuously heated by stellar radiation, keeping the temperature to a roughly constant value, independent of the local density, but varying between simulations. We note that the observed bifurcation of the distribution at high densities in some of the cases is an (unphysical) resolution effect, which largely disappears as the resolution increases – comparing $100-\delta_{20}-f0.25-\epsilon_{6}$ versus $50-\delta_{30}-f0.25-\epsilon_{2.5}$ versus the high-resolution runs CoDa II and $10-\delta_{200}-f0.2-\epsilon_{3.5}$.

A relatively small fraction of the cells have very high temperatures ($T > 10^5$ K), which cannot be reached through photoionization and are instead due to local heating from SNe explosions and strong structure formation shocks. Although rare, these regions are important since they are the origin of the very high y-parameter values we find along certain lines of sight. Such regions are present in all cases,
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Maps of the Comptonization y-parameter for CoDa II, 100 $\epsilon$6, and 100 $\epsilon$8 simulations. Left: full resolution maps sharing the same colour bar scale for y. Middle: maps smoothed with a Gaussian beam of FWHM 1.2 arcmin FWHM, corresponding to the resolution of the 150 GHz channel of the SPT. Right: maps smoothed with a Gaussian beam of FWHM 1.7 arcmin, corresponding to the 95 GHz channel of the SPT. Smoothed maps share the same colour bar scale.</figcaption></figure>

however the detailed distribution depends on the volume, resolution and parameters of the simulation. Some cells in the higher-resolution simulations reach significantly higher temperatures, up to $T > 10^6$ K or even $10^7$ K, which is not the case for the lower-resolution cases, where the energy input is averaged over larger, coarser cells. For the same reason the later cases also lack high-density, high-temperature cells, which are the ones that have the highest gas pressure and thus yield largest y-parameter values. This observation explains the high-y tail in the y-parameter PDF distribution for CoDa-II. The small, high-resolution volumes, which are better resolved than CoDa II, do not reach quite such high y-values due to their much smaller volumes, which limits the statistical sampling of different environments.

Turning our attention to the further cases in Fig. 10, these provide additional insight on how the variations in simulation input parameters affect these phase distributions. Here, we focus on comparing simulations at fixed resolution, as opposed to the effects of varying the numerical resolution, which we have already discussed above.

The first three cases ($25-\delta30-f0.14-\epsilon4$, $25-\delta30-f0.2-\epsilon4$, $25-\delta50-f0.3-\epsilon3$) correspond to early reionization, while the last one ($25-\delta50-f0.5-\epsilon1$) is a late-reionization one. The main variation among

the former cases is in the thickness of the photoionization equilibrium 'line' discussed above. That feature is considerably thicker for $25-\delta30-f0.14-\epsilon4$, the reason being that it completed reionization later ($z \sim 6.2$), than $25-\delta30-f0.2-\epsilon4$ and $25-\delta50-f0.3-\epsilon3$ ($z \sim 6.6$) and thus in the latter cases by $z \sim 6$ many more cells have had time to cool down to lower equilibrium temperatures than those reached immediately post-ionization. Other, minor differences are observed in the maximum temperatures reached ($T > 10^6$ K) in the first two cases versus $T < 10^6$ K for $25-\delta50-f0.3-\epsilon3$, related to the lower star formation efficiency in the latter case. Finally, in the late reionization case, $25-\delta50-f0.5-\epsilon1$, there is a (weak) tail of low-temperature, still neutral cells due to reionization just completing at that time. We also observe similar trends as above, with a thicker equilibrium 'line' due to incomplete equilibration and relatively lower maximum temperatures due to the even lower star formation efficiency in this case.

## 5.4 Second-order Doppler contributions

In Fig. 11, we show the additional contribution to the y-parameter from the quadratic Doppler distortions, both as computed directly
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The angular temperature power spectra of the tSZ signals yielded by our simulations, as listed. Shaded region combines the power spectra from the simulations which are not explicitly shown in the legend. Power spectra calculated for the two SPT band-powers in which the tSZ effect is visible: 150 GHz (**left**) and 95 GHz (**right**). The black solid line is the 1-halo and 2-halo combined contribution of the theoretical y-power spectrum by Bolliet et al. (2018), while the points with error bars are the SPT data for the total CMB power spectrum at the respective frequencies (George et al. 2015).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Distributions of the y-parameter for full CoDa II results (widest distribution, green), uniform global electron temperature (middle distribution, blue), both uniform global electron temperature and uniform ionized fraction (narrowest distribution, orange). Finally, the vertical black dashed line shows the analytical result for instantaneous reionization (Section 4.3, $\log_{10} y = -7.38$, $y = 4.22 \times 10^{-8}$).</figcaption></figure>

from our simulations (unfilled markers) and when accounting for the additional velocity power missing in our simulation volumes (filled markers). In both cases there is a clear positive correlation between the redshift of the end of reionization $z_{xH1} \sim 0.001$ and y from quadratic Doppler. There is a notable scatter around similar $z_{xH1} \sim 0.001$ between the different simulation volumes. However, after accounting for the missing velocity power, this scatter is reduced significantly, indicating that it is largely due to the different box sizes used.

The relationship between this additional Compton-y contribution and the redshift of the end of reionization can be easily understood through equation (25), where the only terms that vary as a function of

$z$ between simulations are $\langle v^2 \rangle$ and $n_e$. After accounting for missing power, the behaviour of $\langle v^2 \rangle$ will be largely similar between all simulations, and so the key determining factor is the evolution of $n_e$ — simulations that reionize earlier will have larger values of $n_e$ throughout cosmic time and thus will have a larger y. The redshift at which reionization ends does not tell the whole story, however. Simulations can complete reionization at the same time but have a different evolution of $n_e$, leading to the remaining scatter we see in the y-parameter for simulations with similar $z_{xH1} \sim 0.001$. For example from Fig. 1 we see that the 100-$\delta$20-$f$0.25-$\epsilon$6 and 50-$\delta$30-$f$0.25-$\epsilon$2.5 simulations have very similar reionization histories, and end up with very similar y-parameters (after accounting for missing velocity power). In contrast to these two simulations, reionization begins more gradually in the 25-$\delta$50-$f$0.5-$\epsilon$1 simulation before ending at a similar time to both 100-$\delta$20-$f$0.25-$\epsilon$6 and 50-$\delta$30-$f$0.25-$\epsilon$2.5. Consequently, the y-parameter is smaller in 25-$\delta$50-$f$0.5-$\epsilon$1 than in 100-$\delta$20-$f$0.25-$\epsilon$6 and 50-$\delta$30-$f$0.25-$\epsilon$2.5.

In summary, we find that, after accounting for missing velocity power, the Compton y-parameter from the quadratic Doppler distortions is of order 10 per cent of the corresponding EoR tSZ signal, with the additional contribution ranging from $1.8 < y/10^{-9} < 3.8$, where the exact value depends on the detailed reionization history. This means the second-order contribution is a sub-dominant, but non-trivial contribution to the EoR tSZ signal.

# 6 SUMMARY AND CONCLUSIONS

The tSZ effect traces the integrated gas pressure history of the universe. The full tSZ signal is generally dominated by the contribution of galaxy clusters at lower redshifts, however we have shown here that it also receives a notable addition from the dense, ionized IGM during the EoR. Consequently, we argued that it is essential to investigate the extent of this effect at this stage in cosmic history.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Phase diagrams showing the gas density versus temperature at the end of the reionization ($z \sim 6$), for selected subset of our simulations: (top left) CoDa II, (top right) 100-$\delta$20-$f$0.25-$\epsilon$6, (bottom left) 50-$\delta$30-$f$0.25-$\epsilon$2.5, and (bottom right) 10-$\delta$200-$f$0.2-$\epsilon$3.5.</figcaption></figure>

A proper, in-depth understanding of the EoR contribution enables us to better model the total tSZ signal and to interpret it correctly. Since the tSZ power spectrum is sensitive to the underlying cosmological model, advancements in our measurements and models will allow for tighter constraints to be placed on cosmological parameters, e.g. $\sigma_8$. In this work, we utilized the data from 16 RAMSES-CUDATON full radiative hydrodynamics cosmological simulations to compute the tSZ signal arising from the EoR. All simulations were well calibrated to yield plausible reionization histories which agree with the available observational constraints. We constructed lightcones of the electron pressure in the redshift range of reionization (roughly $z \sim 6-12$) by interpolating ionized fraction and gas pressure fields between simulation snapshots. From these lightcones, we calculated the Compton y-parameter by numerically integrating along lines of sight in the direction of light propagation. We obtained both high-resolution and SPT beam-smoothed maps of the y-parameter in multiple frequency channels, which we then used to construct PDFs and angular power spectra. We also estimated the contribution to the y-parameter due to quadratic Doppler distortions, correcting for missing velocity power due to small box sizes.

For further understanding and additional tests of the data we evaluated the contributions to the Compton y-parameter from different stages of EoR. We also probed the extent to which the temperature and density fluctuations, and patchiness of the EoR were each

responsible for the values found, and considered in detail the gas equation of state phase diagrams.

All simulations yielded mean Compton y-parameter values in the order of $\langle y \rangle \sim \text{few} \times 10^{-8}$, in rough agreement with previous estimate made by Hill et al. (2015). The magnitude of the tSZ signal originating from the EoR shows clear dependence, roughly proportional, to the redshift of overlap/end of reionization, $z_{0.001}$ (Fig. 4). Simulations which reionize early yield higher y-parameters, by up to factor of 2–3, than those where EoR ended later. For example when the EoR ended at $z_{ov} \sim 6.6$, y-parameter is up to $y \sim 3.5 \times 10^{-8}$, while for $z_{ov} \sim 5.6$, the value of y can be as low as $y \sim 1.22 \times 10^{-8}$. We also note that other parameters which influence the detailed reionization history also influence the tSZ results, increasing the scatter in the mean y-parameter values at a given $z_{ov}$. For example higher star formation efficiency yields a somewhat stronger tSZ effect, likely due to increased supernovae rate. We find that the contribution of quadratic Doppler distortions to the y-parameter is sub-dominant to the electron pressure contribution, at the order of $\sim \text{few} \times 10^9$, or an additional $\sim 10$ per cent contribution to the EoR tSZ effect. The variation in this quadratic Doppler y-parameter is driven by reionization history, where scenarios that reionize earlier have a larger y.

For the purpose of consistency, we attempted to use the same redshift range for the lightcones. However, we note that in doing so, we did not capture the full extent of the EoR for some of our
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Phase diagrams showing the gas density versus temperature at the end of the reionization ($z \sim 6$), for selected subset of our simulations with same resolution, but varying parameters: (top left) 25-$\delta$30-$f$0.14-$\epsilon$4, (top right) 25-$\delta$30-$f$0.2-$\epsilon$4, (bottom left) 25-$\delta$50-$f$0.3-$\epsilon$3, and (bottom right) 25-$\delta$50-$f$0.5-$\epsilon$1.</figcaption></figure>

simulations. As we saw from separating the redshift contributions for the CoDa II simulation, the majority of the tSZ signal arises around the late stages of the EoR. Another limiting factor in this study is the shorter lightcones of the 10 Mpc boxes since we were not able to run one of these simulations to the complete end of the EoR, thus somewhat underestimating the signal in that case.

The tSZ signal from EoR differs significantly from the post-reionization one in both spatial scales and magnitude, as evidenced by the maps and power spectra of the signal that we produced. The mean signal is of order up to few per cent of the cluster one, but its power spectrum shape is quite different, with peak power at smaller scales ($\ell \sim 10^5$ as opposed to $\ell \sim \text{few} \times 10^3$ for post-EoR signal. Overall, the EoR tSZ is a modest contaminant at larger scales, but potentially a dominant one at smaller ones. The Compton y-parameter maps smoothed with Gaussian beams of FWHM 1.2 arcmin and 1.7 arcmin (corresponding to South Pole Telescope (SPT), bands) peak at $y \sim \text{few} \times 10^{-8}$, roughly 3 orders of magnitude lower than the typical cluster values. At these scales, the EoR contribution to the angular power spectrum of the tSZ signal is below 1 per cent that of clusters. However, for the large-volume, high-resolution CoDa II simulation the power spectrum increases by about 1 order of magnitude at smaller scales, while the post-EoR signal decreases significantly, suggesting a more significant and potentially dominant EoR contribution there, although we note that the post-reionization tSZ is not yet fully understood at those scales. The EoR tSZ signal

peaks at such small scales in significant part due to strong local heating by supernovae. Based on the phase diagrams and our results from the density lightcones, we see that the high cell temperatures in dense cells reached by CoDa II are responsible for its locally higher y-parameter values. This tail in its PDF is a result of the spatially resolved cells heated by supernova explosions. Removing the variations in the temperature and ionization fraction fields raises the mean Compton parameter. This is because they widen the range of y-parameter values in the signal. The temperature fluctuations are responsible for the distributions being positively skewed, due to supernova explosions driving up the temperatures of small areas. In conclusion, while we find that the EoR contribution is generally sub-dominant to that of galaxy clusters, it is still essential to obtain an estimate of it as more sensitive technologies are developed in the advancing era of precision cosmology.

## ACKNOWLEDGEMENTS

ITI was supported by the Science and Technology Facilities Council (grant numbers ST/T000473/1, ST/F002858/1, and ST/I000976/1) and the Southeast Physics Network (SEPNet). LC was supported by the Science and Technology Facilities Council (grant number ST/X000982/1). KA was supported by NRF-2021R1A2C1095136 and RS-2022-00197685. The CoDa II simulation was performed at Oak Ridge National Laboratory/Oak Ridge Leadership Com-
<figure><img src="image_14.png" /><figcaption>**Figure 11.** As in Fig. 4, but for the contribution to the y-parameter from second-order Doppler distortions. Unfilled markers correspond to the y-parameter computed using only the velocity of the box, while filled markers include an extra contribution to the velocity to account for missing large-scale power due to small box sizes. Note that the 100-$\delta$20-$f$0.25-$\epsilon$6 and 50-$\delta$30-$f$0.25-$\epsilon$2.5 cases fall directly on top of one another when including the missing power.</figcaption></figure>

puting Facility on the Titan supercomputer, under DOE INCITE award AST031. PRS acknowledges the support of NASA grant No. 80NSSC22K1756. Some of the analysis was done on the Apollo2 cluster at The University of Sussex. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC). We also acknowledge support from Partnership for Advanced Computing in Europe (PRACE) though the DECI-14 funded project subgridEoR. Authors thank Antony Lewis for useful discussions and information on future CMB experiments.

## DATA AVAILABILITY

The data underlying this article are partly a product of The Cosmic Dawn ('CoDa') Project. The derived data generated in this research will be shared upon reasonable request to the corresponding author.

## REFERENCES

* Arnaud M., Pointecouteau E., Pratt G. W., 2005, A&A, 441, 893

* Aubert D., Teysier R., 2008, MNRAS, 387, 295

* Barbosa D., Bartlett J. G., Blanchard A., Oukbir J., 1996, A&A, 314, 13

* Barkana R., Loeb A., 2001, Phys. Rep., 349, 125

* Battaglia N., Bond J. R., Pfrommer C., Sievers J. L., Sijacki D., 2010, ApJ, 725, 91

* Battaglia N., Trac H., Cen R., Loeb A., 2013, ApJ, 776, 81

* Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023

* Birkinshaw M., 1999, Phys. Rep., 310, 97

* Bolliet B., Comis B., Komatsu E., Macías-Pérez J. F., 2018, MNRAS, 477, 4957

* Bouwens R. J. et al., 2015, ApJ, 803, 34

* Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011, MNRAS, 412, 2543

* Carlstrom J. E., Holder G. P., Reese E. D., 2002, ARA&A, 40, 643

* Carlstrom J. E. et al., 2011, PASP, 123, 568

* Challinor A., Lasenby A., 1998, ApJ, 499, 1

* Chluba J., Switzer E., Nelson K., Nagai D., 2013, MNRAS, 430, 3054

* Crawford T. M. et al., 2014, ApJ, 784, 143

* D'Aloisio A., McQuinn M., Davies F. B., Furlanetto S. R., 2018, MNRAS, 473, 560

* da Silva A. C., Barbosa D., Liddle A. R., Thomas P. A., 2001, MNRAS, 326, 155

* Davies F. B. et al., 2018, ApJ, 864, 142

* Eisenstein D. J., Hu W., 1999, ApJ, 511, 5

* Fan X., Carilli C. L., Keating B., 2006, ARA&A, 44, 415

* Faucher-Giguère C.-A., Lidz A., Hernquist L., Zaldarriaga M., 2008, ApJ, 688, 85

* Fixsen D. J., Cheng E. S., Gales J. M., Mather J. C., Shafer R. A., Wright E. L., 1996, ApJ, 473, 576

* George E. M. et al., 2015, ApJ, 799, 177

* Greig B., Mesinger A., 2017, MNRAS, 465, 4838

* Hill J. C. et al., 2014, preprint (arXiv:1411.8004)

* Hill J. C., Battaglia N., Chluba J., Ferraro S., Schaan E., Spergel D. N., 2015, Phys. Rev. Lett., 115, 261301

* Hoag A. et al., 2019, ApJ, 878, 12

* Hogg D. W., 1999, preprint (arXiv Astrophysics e-prints, astro-ph/9905116)

* Horowitz B., Seljak U., 2017, MNRAS, 469, 394

* Hu W., Scott D., Silk J., 1994, Phys. Rev. D, 49, 648

* Iliev I. T., Mellema G., Pen U.-L., Merz H., Shapiro P. R., Alvarez M. A., 2006, MNRAS, 369, 1625

* Iliev I. T., Pen U.-L., Bond J. R., Mellema G., Shapiro P. R., 2007, ApJ, 660, 933

* Iliev I. T., Mellema G., Pen U.-L., Bond J. R., Shapiro P. R., 2008, MNRAS, 384, 863

* Jelić V. et al., 2010, MNRAS, 402, 2279

* Kitayama T., 2014, Prog. Theor. Exp. Phys., 2014, 06B111

* Komatsu E., Seljak U., 2002, MNRAS, 336, 1256

* Lahav O., Lilje P. B., Primack J. R., Rees M. J., 1991, MNRAS, 251, 128

* Mason C. A., Treu T., Dijkstra M., Mesinger A., Trenti M., Pentericci L., de Barros S., Vanzella E., 2018, ApJ, 856, 2

* Mesinger A., McQuinn M., Spergel D. N., 2012, MNRAS, 422, 1403

* Mroczkowski T. et al., 2019, Space Sci. Rev., 215, 17

* Nozawa S., Itoh N., Suda Y., Ohhata Y., 2006, Nuovo Cimento B Serie, 121, 487

* Ocvirk P. et al., 2020, MNRAS, 496, 4087

* Ouchi M. et al., 2010, ApJ, 723, 869

* Park H., Shapiro P. R., Komatsu E., Iliev I. T., Ahn K., Mellema G., 2013, ApJ, 769, 93

* Planck Collaboration VI, 2018, A&A, 641, A6

* Planck Collaboration XVI, 2014, A&A, 571, A16

* Planck Collaboration XXII, 2015, A&A, 594, A22

* Press W. H., Schechter P., 1974, ApJ, 187, 425

* Refregier A., Komatsu E., Spergel D. N., Pen U.-L., 2000, Phys. Rev. D, 61, 123001

* Rephaeli Y., 1995, ARA&A, 33, 541

* Rybicki G. B., Lightman A. P., 1979, Radiative Processes in Astrophysics. Wiley, New York

* Sazonov S. Y., Sunyaev R. A., 2000, ApJ, 543, 28

* Sehgal N., Bode P., Das S., Hernandez-Monteagudo C., Huffenberger K., Lin Y.-T., Ostriker J. P., Trac H., 2010, ApJ, 709, 920

* Seljak U., Burwell J., Pen U.-L., 2001, Phys. Rev. D, 63, 063001
Shaw L. D., Nagai D., Bhattacharya S., Lau E. T., 2010, ApJ, 725, 1452

Sievers J. L. et al., 2013, J. Cosmol. Astropart. Phys., 2013, 060

Sunyaev R. A., Zeldovich Y. B., 1972, Comments on Astrophysics and Space Physics, 4, 173

Sunyaev R. A., Zeldovich Y. B., 1980a, ARA&A, 18, 537

Sunyaev R. A., Zeldovich Y. B., 1980b, MNRAS, 190, 413

Tashiro H., Aghanim N., Langer M., Douspis M., Zaroubi S., Jelic V., 2010, MNRAS, 402, 2617

Teyssier R., 2002, A&A, 385, 337

Toro E. F., 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd edn. Springer-Verlag Berlin

Trac H., Bode P., Ostriker J. P., 2011, ApJ, 727, 94

Wyithe J. S. B., Bolton J. S., 2011, MNRAS, 412, 1926

Zeldovich Y. B., Illarionov A. F., Syunyaev R. A., 1972, Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki, 62, 1217

Zhang P., Pen U.-L., Trac H., 2004, MNRAS, 355, 451

This paper has been typeset from a LaTeX file prepared by the author.