Do Cosmic Rays Heat the Early Intergalactic Medium?

TL;DR

Cosmic rays likely played a significant role in heating the early intergalactic medium before reionization, with implications for understanding early universe magnetic fields and supernova activity.

Contribution

This study models the impact of low-energy cosmic rays on the early IGM, highlighting their heating effect and the importance of CR diffusion and self-confinement.

Findings

CRs heat the IGM by 10-200 K by redshift 10.

CRs contribute negligibly to ionization of the IGM.

The heating pattern depends on CR diffusion and magnetic field structure.

Abstract

Cosmic rays (CRs) govern the energetics of present-day galaxies and might have also played a pivotal role during the Epoch of Reionization. In particular, energy deposition by low-energy ($E \lesssim 10$ MeV) CRs accelerated by the first supernovae, might have heated and ionized the neutral intergalactic medium (IGM) well before ($z \approx 20$) it was reionized, significantly adding to the similar effect by X-rays or dark matter annihilations. Using a simple, but physically motivated reionization model, and a thorough implementation of CR energy losses, we show that CRs contribute negligibly to IGM ionization, but heat it substantially, raising its temperature by $\Delta T=10-200$ K by $z=10$, depending on the CR injection spectrum. Whether this IGM pre-heating is uniform or clustered around the first galaxies depends on CR diffusion, in turn governed by the efficiency of self-confinement due to plasma streaming instabilities that we discuss in detail. This aspect is crucial to interpret future HI 21 cm observations which can be used to gain unique information on the strength and structure of early intergalactic magnetic fields, and the efficiency of CR acceleration by the first supernovae.

Equations

$$t_{\mathrm{ff}}=\sqrt{\frac{3\pi }{32G_N{\rho }_m}},$$

$$M_{\mathrm{Ly}\alpha }(z)\sim 10^8{M}_{\odot }{\left(\frac{10}{1+z}\right)}^{3/2}.$$

$${\bar{V}}_c=24\text{km/s}{\left(\frac{M_{\mathrm{rf}}}{10^8{\text{M}}_{\odot }}\right)}^{1/3}{\left(\frac{1+z}{10}\right)}^{1/2}.$$

$${\dot{\rho }}_{\ast }(z)=f_{\ast }\frac{{\Omega }_b}{{\Omega }_m}{\int }_{M_{\min }(z)}^{\infty }d{M}_h\frac{M_h}{t_{\mathrm{ff}}(M_h)}\frac{dN}{d{M}_{h}dV},$$

$$\frac{d{x}_{\mathrm{HII}}}{dz}=\frac{dt}{dz}\left[x_{\mathrm{HI}}{\Gamma }_{\mathrm{HI}}-R\right].$$

$${\Gamma }_{\mathrm{HI}}(z)={\int }_{{\nu }_0}^{\infty }d\nu {\lambda }_{\mathrm{HI}}(z,\nu ){\sigma }_{\mathrm{HI}}(\nu ){\dot{n}}_{\gamma }(z,\nu ),$$

$${\lambda }_{\mathrm{HI}}(z,\nu )={\lambda }_{{\nu }_0}(z){\left(\frac{\nu }{{\nu }_0}\right)}^{3/2},$$

$${\lambda }_{{\nu }_0}(z)\approx 39{\left(\frac{1+z}{4}\right)}^{-5}\mathrm{Mpc}.$$

$${\dot{n}}_{\gamma }(z,\nu )=f_{\mathrm{esc}}\frac{d{N}_{\gamma }}{dMd\nu }{\dot{\rho }}_{\ast }(z)(1+z{)}^3;$$

$${\Gamma }_{\mathrm{HI}}(z)$$

$$R(z)$$

$$\frac{d{T}_k^{\mathrm{i}}}{dz}={\frac{d{T}_k^{\mathrm{i}}}{dz}|}_{\mathrm{ex}}+{\frac{d{T}_k^{\mathrm{i}}}{dz}|}_{\mathrm{ion}}+{\frac{d{T}_k^{\mathrm{i}}}{dz}|}_{\mathrm{ph}},$$

$${\frac{d{T}_k^{\mathrm{i}}}{dz}|}_{\mathrm{ex}}=\frac{2T_k^{\mathrm{i}}}{1+z},$$

$${\frac{d{T}_k^{\mathrm{i}}}{dz}|}_{\mathrm{ion}}=-\frac{T_k^{\mathrm{i}}}{1+x_e}\frac{d{x}_e}{dz},$$

$${\frac{d{T}_k^{\mathrm{i}}}{dz}|}_{\mathrm{ph}}=\frac{2}{3k_B(1+x_e)}\frac{dt}{dz}{\mathcal{H}}_{\ast },$$

$${\mathcal{H}}_{\ast }(z)=x_{\mathrm{HI}}(z){\int }_{{\nu }_0}^{\infty }d\nu {\lambda }_{\mathrm{HI}}(z,\nu ){\sigma }_{\mathrm{HI}}(\nu ){\dot{n}}_{\gamma }(z,\nu )h\left(\nu -{\nu }_0\right),$$

$$\begin{array}{rl}{\mathcal{H}}_{\ast }(z)& =\left(\frac{\beta -1}{{\beta }^2-1/4}\right)x_{\mathrm{HI}}(z){\lambda }_{{\nu }_0}(z){\sigma }_0{f}_{\mathrm{esc}}h{\nu }_0\frac{d{N}_{\gamma }}{dM}{\dot{\rho }}_{\ast }(z)(1+z{)}^3\\ & \approx 3\times 10^{-19}{\mathrm{erg}\mathrm{s}}^{-1}{x}_{\mathrm{HI}}(z)\left(\frac{f_{\mathrm{esc}}}{10^{-2}}\right)\\ & \times \left(\frac{{\dot{\rho }}_{\ast }(z)}{M_{\odot }{\mathrm{Mpc}}^{-3}{\mathrm{yr}}^{-1}}\right)(1+z{)}^{-2}.\end{array}$$

$${\tau }_e(z)={\int }_0^z{n}_e(z^{\prime }){\sigma }_T\frac{dl}{d{z}^{\prime }}d{z}^{\prime },$$

$${\dot{E}}_{\mathrm{p}}(z)=ϵE_{\mathrm{SN}}\mathrm{SNR}(z)(1+z{)}^3,$$

$$f_{\mathrm{SN}}=\frac{{\int }_8^{50}\varphi (m)dm}{{\int }_{0.1}^{100}m\varphi (m)dm}\sim 10^{-2}{M}_{\odot }^{-1},$$

$$\varphi (m)\propto m^{-1+x}\exp \left(-\frac{m_c}{m}\right),$$

$$\begin{gathered} {\dot{E}}_{\mathrm{p}}(z)\sim 10^{-33}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}{\mathrm{s}}^{-1}\left(\frac{ϵ}{0.1}\right)\left(\frac{E_{\mathrm{SN}}}{10^{51}\mathrm{erg}}\right) \\ \left(\frac{f_{\mathrm{SN}}}{10^{-2}{M}_{\odot }^{-1}}\right)\left(\frac{{\dot{\rho }}_{\ast }(z)}{M_{\odot }{\mathrm{yr}}^{-1}{\mathrm{Mpc}}^{-3}}\right)(1+z{)}^3. \end{gathered}$$

$$q_p(E,z)=\frac{C(z)}{\beta (E)}{\left(\frac{E^2+2E{m}_p{c}^2}{E_0^2+2E_0{m}_p{c}^2}\right)}^{-\frac{\alpha }{2}},$$

$${\dot{E}}_{\mathrm{p}}(z)={\int }_{E_{\min }}^{E_{\max }}E{q}_p(E,z)dE.$$

$${-\frac{dE}{dt}|}_I=\frac{4\pi e^4}{m_e\beta c}\sum _{Z}Z{x}_{\mathrm{HI}}{n}_Z\left[\ln \left(\frac{2m_e{c}^2}{I_Z}{P}^2\right)-{\beta }^2\right],$$

$${-\frac{dE}{dt}|}_C=\frac{4\pi e^4{n}_e}{m_e\beta c}\left[\ln \left(\frac{2m_e{c}^2\beta }{\hslash {\omega }_{pl}}P\right)-\frac{{\beta }^2}{2}\right],$$

$${-\frac{dE}{dt}|}_a=\frac{E(E+2m_p{c}^2)}{E+m_p{c}^2}\frac{1}{1+z}\frac{dz}{dt}.$$

$$t_i=\frac{E}{d{E}_i/dt}\le t_H,$$

$$t_H(z)={\int }_{\infty }^{z}d{z}^{\prime }\frac{dt}{d{z}^{\prime }}\simeq \frac{2(1+z{)}^{-3/2}}{3H_0{\Omega }_m^{1/2}}\simeq 0.2{\left(\frac{1+z}{21}\right)}^{-3/2}\mathrm{Gyr}.$$

$$\frac{\partial }{\partial t}{N}_p+\frac{\partial }{\partial E}\left(b{N}_p\right)+\frac{N_p}{t_D}=Q_p(E,z),$$

Full text

††pagerange: Do cosmic rays heat the early intergalactic medium?–References††pubyear: 2017

1 Introduction

Sometimes during the first billion cosmic years the first stars formed. Ionizing photons emitted by these sources induced a major phase transition from the otherwise cold and neutral state in which the intergalactic medium (IGM) was left after recombination ($z\lesssim 500$) to the warm/ionized state we measure today. The transition phase is known as the Epoch of Reionization (EoR). Understanding the detailed EoR physics is one of the primary goals of present-day cosmology (Barkana & Loeb, 2001; Zaroubi, 2013; Ferrara & Pandolfi, 2014; Mesinger, 2016).

To date, we only have a broad-brush picture of such important event. Nevertheless sizable steps forward have been made in the last decade. High redshift ($z\lesssim 6$) quasar absorption line experiments, probing the physical state of the neutral IGM, indicate that reionization was complete by $z\sim 5.7$ (Fan et al., 2006). The CMB polarization measurements made by PLANCK (Planck Collaboration et al., 2016) have set tight bounds to the value of the free electron scattering optical depth, $\tau_{\rm e}=0.066\pm 0.016$. Translating this value into a reionization redshift, $z_{r}$, yields uncertain results as the conversion depends from the (yet unknown) reionization history. In the popular instantaneous reionization model it is $7.8<z_{r}<8.8$, with an upper limit on the EoR duration of $\Delta z<2.8$. Thus, reionization might have been a remarkably fast ($\leq 400$ Myr) process.

Great advances in the understanding of the EoR history and physics are expected from a number of upcoming observations of the redshifted HI 21-cm line signal from these epochs (Ciardi & Madau, 2003; Furlanetto et al., 2006; Morales & Wyithe, 2010; Pritchard & Loeb, 2012).

Several experiments are attempting to measure the 21-cm signal from the EoR using low-frequency radio interferometers. These include GMRT111http://gmrt.ncra.tifr.res.in, LOFAR222http://www.lofar.org, MWA333http://www.mwatelescope.org, PAPER444http://eor.berkeley.edu, HERA555http://reionization.org and, in the future, SKA666http://www.skatelescope.org. Measuring such signal will provide direct information on the pre-reionization ages. The 21-cm signal will inform us of the HI spin temperature as a function of redshift. As such it will probe the efficiency of gas heating and ionization mechanisms that might have played a role from recombination to the time at which the first light sources appeared (the so-called ”Dark Ages”).

The first stars and galaxies are likely to be the main sources of the ionizing UV  (Ciardi & Ferrara, 2005). Both theoretical arguments and numerical studies suggest that the first generation stars (known as Pop III stars) were more massive than present day ones and metal-free. These short-lived stars polluted the surrounding gas with metals, inducing a rapid transition to a cosmological star formation rate (SFR) dominated by the present-day, Pop II/I, stars (Schneider et al., 2002, 2006).

However, sources of higher energy X-ray photons might also have been present. X-rays have a far longer mean free path than lower energy UV photons. Therefore, these photons are able to travel significantly larger distances in the neutral IGM and release a significant amount of energy, eventually increasing the temperature of the intergalactic gas. Potential X-ray sources are quasars, supernovae and X-ray binaries. However, very little is known about their abundances, evolution and spectra, especially at these high redshifts (McQuinn, 2012).

An additional contribution to IGM heating has been proposed by Chang et al. (2012); Chang et al. (2016). The deposition of kinetic energy into the IGM via plasma instabilities (Schlickeiser et al., 2012) triggered by TeV photons from blazars may yield a higher heating rate than photoheating for $z\lesssim 6$. Such an effect is more relevant in the intergalactic voids, i.e. in less dense regions, of the ionized medium and allows to better reproduce Lyman-$\alpha$ forest observables over the redshift range $2<z<3$ (Puchwein et al., 2012). On the other hand, recent numerical simulations have been suggesting that this mechanism could be fairly ineffective in the IGM conditions (Sironi & Giannios, 2014), although we note that there is still no general consensus on this.

Similarly to X-ray photons, cosmic rays (CRs) – accelerated from the shocks created by the exploding supernovae and promptly released in the IGM – could also efficiently deposit their energy in the form of gas heating (Nath & Biermann, 1993; Sazonov & Sunyaev, 2015). Curiously, such heating source has received relatively little attention in the literature so far, in spite of the fact that it is known since decades that in local galaxies low-energy CRs play a key overall role in the energetics of these systems and, in particular, in regulating the ionization and thermal state of the interstellar gas.

At higher redshifts, CRs should be able to escape their host galaxies by advection or by diffusion before they lose a significant fraction of their energy within the halo (Rollinde et al., 2008). As they propagate in the IGM, CRs interact with the surrounding environments mainly via H/He photoionization and Coulomb collisions with free electrons, and by doing so they deposit thermal energy in the gas. Both mechanisms imply that CRs could contribute to the thermal history of the IGM (Samui et al., 2005).

Earlier studies of the impact of CRs on the high-redshift IGM mainly focused on the cosmological CRs originated by Pop III stars (e.g., Rollinde et al., 2006). Recently, the authors of Sazonov & Sunyaev (2015) found that low-energy particles ($E\lesssim 30$ MeV) are capable of increasing the IGM temperature to 10-100 K before standard heating sources, such as galaxies and quasars, appear.

In the present work, we reanalyze the role played by CRs from high-$z$ Pop II stars. In fact, it has been shown that in realistic models of galaxy formation, chemical feedback suppresses metal-free star formation in the self-enriched progenitors, and although Pop III star formation can in principle persist down to $z\sim 3-4$, Pop II stars dominate the SFR at any redshift (Tornatore et al., 2007; Pallottini et al., 2014). Moreover, Pop II stars are observed in the local Universe and their properties can be much more robustly constrained.

In §2, we build a simple reionization model and link it to the cosmic star formation rate. In §3 we introduce the treatment of CR production, energy losses and propagation. With these ingredients, in §4 we compute the CR contribution to IGM heating, and explore the implications of the spatial dependence of the temperature increment. Results and assumptions are discussed in §5. We use the cosmological parameters: $h=0.678$, $\Omega_{m}=0.308$, $\Omega_{\Lambda}=0.692$, $\Omega_{b}h^{2}=0.0223$, $n_{s}=0.968$ and $\tau_{e}=0.066$ (Planck Collaboration et al., 2016).

2 Reionization by early galaxies

In this Section, we assume that galaxies are the primary reionization sources. This scenario is supported by observations of high-redshift galaxies, which would be able to reionize the Universe by $z=6$, provided that a substantial fraction of their ionizing emission escapes into the IGM (Robertson et al., 2010).

Our model implements the basic physical processes required to properly model IGM evolution in the presence of the ionizing radiation from galaxies. Such formalism allows us to track the IGM thermal and ionization history, and to reproduce the available observational EoR data. A more detailed treatment can be obtained by means of numerical simulations, e.g. in Iliev et al. (2007); Ciardi et al. (2012); Gnedin (2016).

2.1 Star Formation Rate

We assume that the star formation rate per unit of stellar mass and comoving volume (SFR) inside a DM halo is proportional to its mass, $\dot{\rho}_{*}(M_{h})\propto M_{h}$, and it occurs on a free-fall time-scale,

[TABLE]

where $G_{N}$ is the gravitational constant and $\rho_{m}$ is the average mass density inside the virial radius of the halo (Bryan & Norman, 1998). We further assume that stars only form in Ly$\alpha$ cooling halos (Barkana & Loeb, 2000), i.e. those above a mass

[TABLE]

Radiative feedback is expected to quench SF in halos with circular velocity at the virial radius, $V_{c}$ smaller than a critical value $\bar{V}_{c}$. The corresponding minimum mass, $M_{\rm rf}(z)$, is obtained according to the relationship (Madau et al., 2001)

[TABLE]

We can then write the cosmological SFR per unit comoving volume at redshift $z$ as

[TABLE]

where $f_{*}$ is the star-forming efficiency, $M_{\rm min}$ is the maximum between $M_{\rm Ly\alpha}$ and $M_{\rm rf}$. $\Omega_{b}$ and $\Omega_{m}$ are the density parameters of baryonic and total matter, respectively, and $dN/(dM_{h}dV)$ gives the number density of haloes within a mass range ($M_{h},M_{h}+dM_{h}$). To estimate the DM halo mass function we use the Press-Schechter formalism augmented by the Sheth-Tormen correction for ellipsoidal collapse (Sheth et al., 2001).

As we show in Fig. 1, a value of $f_{*}=0.02$ and $\bar{V}_{c}=100$ km/s allows us to reproduce the SFR measurements reported by Madau & Dickinson (2014).

2.2 Reionization history

We denote the neutral hydrogen (HI) fraction by $x_{\rm H\small{I}}=n_{\rm H\small{I}}/n_{\rm H}$ and the ionized fraction by $x_{\rm H\small{II}}=n_{\rm H\small{II}}/n_{\rm H}$ with $x_{\rm H\small{I}}+x_{\rm H\small{II}}=1$ where $n_{x}$ is the number density of species $x$.

It is further assumed that the ionization fraction of singly ionized helium and hydrogen are equal, $x_{\rm He\small{II}}=x_{\rm H\small{II}}$.

We can then write the evolution equation for the ionization fraction in terms of the photo-ionization ($\Gamma_{\rm HI}$) and the recombination ($R$) rates as

[TABLE]

Since stellar radiation contributes to ionize the IGM with photons of energy higher than the ionization threshold $I_{H}=h\nu_{0}=13.6$ eV, one can write the ionization rate as (Choudhury & Ferrara, 2005)

[TABLE]

where the ionization cross-section is $\sigma_{\rm HI}(\nu)=\sigma_{0}(\nu/\nu_{0})^{-3}$, with $\sigma_{0}=6.3\times 10^{-18}$ cm2.

The mean free path of hydrogen ionizing photons depends on the distribution of absorbing Lyman limit systems, which can be computed through the distribution of the column density $N_{H\small{I}}$ (Petitjean et al., 1993). Assuming the distribution of absorbers to be given by $f(N_{H\small{I}})dN_{H\small{I}}\propto N_{H\small{I}}^{-3/2}dN_{H\small{I}}$ (Faucher-Giguère et al., 2008), the mean free path increases with frequency as in Miralda-Escudé (2003),

[TABLE]

with $\lambda_{\nu_{0}}(z)$ depending on the size and topology of the ionized regions.

In Schirber & Bullock (2003), $\lambda_{\nu_{0}}(z)$ is derived from the column density distribution of the Lyman-limit systems and is found to be a rapidly evolving function of $z$. We then assume

[TABLE]

Finally in Eq. (6), the proper specific density rate of ionizing photons production is given by

[TABLE]

$f_{\rm esc}$ is the escape fraction of ionizing photons from galaxies, and $dN_{\gamma}/(dMd\nu)$ the specific number of photons produced per unit mass of Pop II stars formed. Here we assume that the UV stellar spectrum is a power law $\propto\nu^{-\beta}$ and that, integrating over frequency, $dN_{\gamma}/dM=8.05\times 10^{60}M_{\odot}^{-1}$ (Choudhury & Ferrara, 2005).

Putting all together, Eq. (6) can then be simplified to

[TABLE]

where $\beta=5$ is a typical value for Pop II stars.

The recombination rate can be expressed as

[TABLE]

where $\alpha_{A}$ is the case-A recombination coefficient, $n_{e}(z)=x_{\rm HII}(z)(1+\chi_{\rm He})n_{\rm H}(z)$ is the total electron number density and $\chi_{\rm He}$ the cosmic helium fraction (in density). For the clumping factor $C\equiv\langle n_{\rm H\small{II}}^{2}\rangle/\langle n_{\rm H\small{II}}\rangle^{2}\approx 2$ we take the fiducial average given by Finlator et al. (2011).

The evolution of the gas temperature in the ionized regions ($T_{k}^{\rm i}$) is given by the combination of all the cooling and heating processes (see, e.g., Nath & Silk, 2001)

[TABLE]

where the cooling due to the expansion of the universe can be written as

[TABLE]

and the heating due to the change in the internal energy, corresponding to the change in the total number of gas particles due to ionizations, where He ionizations are assumed to be negligible, is

[TABLE]

with $x_{e}=n_{e}/n_{H}$. The last term is the heat gain by the gas particles from the surrounding radiation field

[TABLE]

where $k_{B}$ is the Boltzmann constant and $\mathcal{H}_{*}$ is the photoheating rate per baryon.

Analogously to Eq. (6), the latter can be written as

[TABLE]

which we can integrate over frequency, obtaining

[TABLE]

Finally, we choose the value for $f_{\rm esc}$ by assuming that the total optical depth, defined as

[TABLE]

does not exceed the $3\sigma$ value measured by PLANCK, $\tau_{e}=0.066\pm 0.016$ (Planck Collaboration et al., 2016). In Fig. 2 we show the optical depth corresponding to $f_{\rm esc}=10^{-2}$ and compared to the observed value for $\tau_{e}$. We also show that our model predicts a fully reionized IGM by $z\sim 6$ as inferred from the Gunn-Peterson trough detections (Fan et al., 2006).

Finally, we note that a very similar value for $f_{\rm esc}$ has been found by means of a more sophisticated approach and tested against more observables by the authors of Choudhury & Ferrara (2005).

3 Cosmic rays in the IGM

CRs accelerated in early galaxies can act in principle as an additional source of non-thermal energy for the IGM. In the Milky Way, most of the CR energy is in protons. They diffuse or advect out of the Galaxy on timescales of about 30 Myr that can be directly inferred from secondary-over-primary ratios (Blasi, 2013), with only a few percent of the energy lost in pion production and ionization (Strong et al., 2007). Thus, a large fraction of the power injected in CRs in our Galaxy ends up in the surrounding IGM.

In our model CRs are accelerated by star-forming galaxies with an universal energy spectrum and their energy released far beyond the circumgalactic gas, i.e. in the IGM. This conclusion is motivated by the fact that earlier structures are expected to be less confining than the present galaxies, since they were smaller and had a weaker magnetic field. In fact, Lacki (2015) and Rollinde et al. (2006) argued that primary CRs escape from parent galaxies on a timescale short enough so that they do not suffer any energy loss.

To follow the propagation of an homogeneous CR population in an expanding universe for a continuous source of CRs, we generalize the classical work of Montmerle (1977), including all the relevant energy loss processes.

3.1 Production in the early galaxies

Star formation pumps energy into CR protons at a rate

[TABLE]

where $E_{\rm SN}\sim 10^{51}$ erg is the average explosion energy for a Type II supernova (SN) not going in neutrinos, $\epsilon\sim 0.1$ is the fraction of the kinetic energy transferred to CRs by a single SN, and SNR is the comoving SN rate. In principle, one should account also for Helium nuclei. For simplicity, we assume here that $\alpha$-particles can be treated as four protons, and hence be absorbed in the proton spectrum efficiency.

In addition to the SFR, to derive SNR we need to know the number of SNe explosions per solar mass of forming stars, which is given by

[TABLE]

where $\phi(m)$ is the Initial Mass Function (IMF) of Population II/I stars, for which we assume the following form

[TABLE]

with $x=-1.35$, $m_{c}=0.35$ $M_{\odot}$; $m$ lies in the range $[0.1,100]$ $M_{\odot}$ (Larson, 1998). We neglect the contribution from Pop III stars. Their SF history is very debated and still highly uncertain. In fact, detailed studies, exploiting cosmological hydrodynamical simulations implementing chemical feedback effects, have shown that Pop II/I stars dominate the global SFR at any redshift (Tornatore et al., 2007; Maio et al., 2010). Moreover, the recent surveys hunting for PopIII stars in the Milky Way have found no metal free stars so far, implying that they are rare in the Milky Way even if they exist (Beers & Christlieb, 2005).

Combining the above formulae, Eq. (19) becomes

[TABLE]

Particle acceleration in SN explosions is believed to occur through diffusive shock acceleration, which leads to momentum power-law spectra of accelerated particles. With this in mind, we assume that the source function of volume averaged CR protons injected by SNe (defined as a rate per unit energy and volume) is

[TABLE]

where $E$ is the proton kinetic energy, $E_{0}=1$ GeV, $m_{p}$ is the proton mass, $\beta=v/c$ is the dimensionless velocity of the particle, $\alpha\geq 2$ is the slope of the differential spectrum of accelerated particles and $C(z)$ is a redshift-dependent normalization obtained by imposing that the total kinetic energy rate equals $\dot{E}_{\rm p}(z)$, i.e.

[TABLE]

In Eq. (24) we fix $E_{\rm max}=10^{6}$ GeV (Lagage & Cesarsky, 1983) and we verify a posteriori that our conclusions are not strongly dependent on our choice of $E_{\rm min}=10$ keV. From Eq. (23) one can easily realize that $\alpha$ is the key parameter determining the fraction of the total kinetic energy released that goes into protons with $E\ll 1$ GeV. We keep $\alpha$ as the only free parameter of the model.

In Fig. 3 we plot the source function as a function of the kinetic energy for three different values of $\alpha$. We found that protons with kinetic energies below $10$ MeV represent $0.5$%, $2.7$%, $13$% of the total kinetic energy released for $\alpha=2,\,2.2,\,2.5$, respectively.

3.2 Energy losses in the IGM

CRs can be an efficient heat source especially for a low density gas. When the CR proton ionizes an atom, it transfers a certain fraction of its kinetic energy to the electron, which is either used for further atomic excitation and ionization, or distributed via elastic collisions to other species of the medium. The latter process increases the kinetic temperature of the gas.

Ionization losses can be taken into account using the Bethe-Bloch equation, that for $\gamma\ll m_{p}/m_{e}$ can be approximated as (Enßlin et al., 2007)

[TABLE]

where $m_{e}$ is the electron mass, $n_{Z}$ is the number density of the elements with atomic number $Z$, $I_{Z}$ is the ionization potential ($I_{\rm H}=13.6$ eV and $I_{\rm He}=24.6$ eV), and $P=p/(m_{p}c^{2})=\sqrt{\gamma^{2}-1}$ is the dimensionless particle momentum.

Losses due to Coulomb interactions, which describe the fact that the energy lost by protons in Coulomb interactions is directly transferred to momentum of the plasma electrons (and hence heating), can be expressed as (Gould, 1972)

[TABLE]

where $\omega_{pl}=(4\pi e^{2}n_{e}/m_{e})^{1/2}$ is the plasma frequency. The number density of free electrons $n_{e}$ is computed from our reionization model described in §2.2 assuming only stellar radiation ionizations. Indeed, we verify a posteriori that CR ionizations are a subdominant contribution (see §4.1).

Inverse Compton scattering with respect to CMB photons can be safely neglected, since its timescale is much longer compared to collisional processes (Evoli et al., 2012).

Finally, the adiabatic energy losses caused by Hubble expansion can be taken into account as (Montmerle, 1977)

[TABLE]

For CR energy deposition to be effective, its timescale must be shorter than the Hubble time,

[TABLE]

where

[TABLE]

In Fig. 4 we plot $t_{i}/t_{H}$ as a function of redshift for the loss mechanisms discussed above. Energy losses are efficient for kinetic energies $\lesssim 10$ MeV. Ionization losses dominate over Coulomb losses at earlier epochs when the IGM was mainly neutral.

3.3 Propagation in the IGM

The evolution equation of the CR proton number density (averaged over the volume), $n_{p}(E,z)$, can be written as follows (Montmerle, 1977; Rollinde et al., 2006; Evoli et al., 2008)

[TABLE]

where the number density of protons and the source term are now normalized to $n_{\rm H}(z)$, being $N_{p}=n_{p}/n_{H}$ and $Q_{p}=q_{p}/n_{H}$. We also assume that the total energy loss rate $b\equiv dE/dt$ is given by the sum of the loss processes described by Eqs. (25), (26) and (27).

Proton-proton interactions cause CR energy loss in a timescale (Gabici et al., 2007)

[TABLE]

where the inelasticity coefficient is $\kappa\approx 0.45$ and $\sigma_{\rm pp}\approx 35$ mb is the total inelastic cross-section for proton-proton interaction. This results in $t_{\rm pp}\approx 10^{9}(1+z)^{3}$ Gyr, making this process negligible for our purposes.

Eq. (30) is solved numerically using the Crank-Nicolson implicit method described in Press et al. (1992). The results are presented in Fig. 5, showing the redshift evolution of the proton spectrum. The effect of energy losses (solid vs. dashed lines) is evident mainly at low-energies ($E\lesssim 10$ MeV).

3.4 IGM ionization and heating

We come now to the central question: can CR energy losses sensibly affect the IGM ionization state and/or temperature? The primary ionization rate for H is

[TABLE]

where $W_{\rm H}\simeq 36.3$ eV is the mean energy expended by a CR proton to create an ion pair (Jasche et al., 2007).

Following Spitzer & Scott (1969), we account for all secondary and higher generation ionizations by multiplying the primary ionization rate in Eq. (32) by a factor $\xi(x_{e})$. For this we use a linear interpolation between the extreme case $\xi(1)=1$ and $\xi(0)=5/3$, which reads as

[TABLE]

In the case of Coulomb losses, we can assume that all the lost energy is entirely converted into background heat. The corresponding heating rate can be calculated by using Eq. (26)

[TABLE]

The contribution to heating by secondary electrons from ionization can be divided in three regimes according to their energy: for $E>I_{\rm H}$ ionization or excitation of H I can occur; for $3\,I_{\rm H}/4<E<I_{\rm H}$ the electron can suffer losses by Coulomb and excitation collisions; for $E<3\,I_{\rm H}/4$ the energy is transferred directly into heating.

An approximate general formula for the heating is given by Jasche et al. (2007), leading to a total heating rate by CR as

[TABLE]

It follows that, in a neutral medium, a heat input of $\Delta E=W_{\rm H}-5/3\,I_{\rm H}\sim 13.6$ eV for every ionization of hydrogen via CR protons is transferred to the IGM. We note, however, that the above expression is likely to overestimate the heating rate as electron energy losses via excitations are not accounted for.

4 Results

We are now ready to discuss the effects of CRs on the IGM ionization fraction and temperature from our model.

4.1 Impact of CRs on reionization

The ionization rate computed with Eq. (32) is shown in Fig. 6 and compared to the ionization rate by UV photons through Eq. (6). The CR ionization rate is several orders of magnitude smaller than the UV photoionization rate. This justifies the fact that CR were not included when we described our reionization model in § 2.2. This result can be better understood by comparing the stellar and CR emissivities in a more simplified scenario. We recall that the ionizing photon emissivity by galaxies is given by

[TABLE]

where the energy in photons is given by

[TABLE]

where $\dot{N}_{\gamma}\sim 5\times 10^{47}$ s*-1* is the rate of ionizing photons, and $t_{*}\sim 30$ Myr is the stellar lifetime. In doing so, we are assuming that the UV emission by galaxies is dominated by the same stars that go supernova, with mass $\sim 10$ M⊙.

On the other hand, the CR emissivity can be written as

[TABLE]

where $f_{d}\sim 10^{-3}$ is the fraction of energy deposited in the IGM and, as we discussed in § 3.2, corresponds to the fraction of energy in CR protons with $E\lesssim 10$ MeV. We are thus assuming that all the deposited energy is used to ionize IGM atoms.

The ratio between the corresponding fluxes is then

[TABLE]

where $\lambda_{i}$ designates the corresponding mean free path.

Finally, the ratio of the ionization rates, $J_{i}/\tau_{i}$, with $\tau_{i}$ the corresponding energy loss time, can then be roughly estimated as $10^{-3}$ by approximating $\tau_{i}\propto\lambda_{i}$.

4.2 IGM heating

While there is a consensus that UV stellar radiation is largely responsible for cosmic reionization, its impact on the global IGM temperature is limited to fully ionized regions. For quasi-neutral IGM regions, X-rays are clearly more relevant due to their larger mean free path. These energetic photons might come from an early population of relatively soft X-ray binaries  (Mesinger et al., 2013). Alternatively, if sourced by black holes, they could pre-heat the gas up to $10^{4}$ K  (Ricotti & Ostriker, 2004). However, uncertainties related to the nature and abundance of their sources at high redshift make predictions very uncertain (McQuinn, 2012; Pritchard & Furlanetto, 2007). For this reason, CRs might represent a competitive, alternative source of thermal input for the neutral IGM.

The IGM temperature increase produced by CR heating (Fig. 7) is

[TABLE]

where $\mathcal{H}_{\rm h}^{\rm CR}$ is given by Eq. (35). The IGM temperature can be raised up to $\sim 3\times 10^{3}$ K before reionization is complete and it exceeds the CMB one, $T_{\rm CMB}(z)=2.725\,(1+z)$ K, at $z\lesssim 9(12)$ for $\alpha=2(2.5)$. These results imply that the IGM is pre-heated well before being reionized. Depending on the efficiency of the CR diffusion mechanism, pre-heating might be confined in regions around star forming galaxies, or, if diffusion is very efficient, it might give rise to a more distributed, quasi-uniform warmer floor (see § 4.3).

We can directly compare our predictions with the results presented in Fig. 1 of Sazonov & Sunyaev (2015). These authors also study the IGM heating from low-energy CRs, finding a temperature increament between (1-103) K, mostly depending on the minimum halo mass allowed to form stars. The highest temperature was found for a supernova energy explosion $E_{\rm SN}=10^{53}$ erg (hence corresponding to a PopIII supernova) out of which 5% is pumped into low-energy CRs, and a minimum star-forming halo mass of $3\times 10^{5}\,M_{\odot}$. Our model predicts a similar $\Delta T$ but relying on fairly standard stellar populations (and supernova explosion energies) whose SFR has been calibrated with both the observed cosmic star formation history and reionization constraints from the Thomson scattering optical depth.

4.3 Diffusion in the IGM

As energetic particles are injected into the IGM, they will be subject to a random-walk through it. The travelled distance depends on the strength and structure of the intergalactic magnetic field, on which very little is known at high redshift. The standard assumption is to consider CR energy deposition as a uniform background (see, e.g., Lacki, 2015). The rationale behind this assumption is based on the following argument.

The slowest rate at which CRs can diffuse corresponds to the so-called Bohm diffusion. This assumes one scattering per gyroradius; the diffusion coefficient is $D_{B}=cr_{L}/3$, where $r_{L}$ is the particle Larmor radius. The maximum diffusion timescale between haloes is then given by

[TABLE]

where $\langle d\rangle$ is the average proper distance between them, and the Bohm diffusion coefficient for protons ($Z=1$) can be estimated as

[TABLE]

$B_{0}=10^{-16}$ G is the assumed IGM magnetic field strength at $z=20$, following Sazonov & Sunyaev (2015).

To estimate $\langle d\rangle$, we consider uniformly distributed haloes. Their average inter-distance is then

[TABLE]

where $dN/dM_{h}$ is the comoving halo density. At $z=20$, the haloes contributing mostly to the SFR are those with $M_{h}\sim M_{\rm min}(z=20)$, whose mean separation is $\langle d\rangle\sim 50$ kpc. From Eq. (41) and Eq. (29), we deduce that $t_{B}\lesssim t_{H}$ as long as the CR energy is $E\gtrsim 20$ keV. This result would support the idea of a uniform thermal deposition by CRs.

However, the above calculation is incomplete, as CRs may affect the environment in which they propagate. When CRs escape from the halo they produce an electric current to which the background plasma reacts by generating a return current that in turn leads to the development of small scale instabilities. The growth of such instabilities leads to large turbulent magnetic fields and to an enhanced particle scattering. In short, CRs may undergo self-confinement (Blasi et al., 2015). In this scenario, the particle diffusion timescale can be significantly larger than what we found in Eq. (41).

In order to get an estimate of the potential maximum effect associated with this mechanism, we generalize the formalism developed by Blasi et al. (2015) to the non-relativistic regime; moreover, we maximize the effect by assuming that all escaping CRs contribute to the self-generated magnetic field.

The differential number density (in momentum) of CRs escaping out from a halo at a distance $r$ from it can be written as777The injection spectrum assumed here is equivalent to Eq. (23), since $dN/d^{3}p\propto p^{-4}$ corresponds to $dN/dE\propto p^{-2}$.

[TABLE]

where $A(r)$ is obtained by imposing that the total pressure exerted on a surface $S=4\pi r^{2}$ by the CR source,

[TABLE]

is given by the CR pressure at the same distance

[TABLE]

where the source luminosity in CRs, $L_{\rm CR}$, for a typical star forming halo at $z=20$ is

[TABLE]

By equating Eqs. (45) and (46), one has

[TABLE]

The electric current $j_{\rm CR}$ associated with CRs streaming away from their sources can be written as

[TABLE]

where we have introduced

[TABLE]

Assuming that the non-resonant modes are able to grow on a timescale much shorter than $t_{H}$, the magnetic field saturates at a value $\delta B_{\rm s}$. The saturation level is set by equipartition between the energy density of the amplified field and the kinetic energy density of the CR current (see Eq. (49)):

[TABLE]

the last approximate equality holds since $pg(p)$ reaches a maximum at $\sim$GeV/c and remains constant for larger momenta. Numerically this yields $\delta B_{\rm s}\approx 0.01\mu$G at $r=1$ kpc.

The corresponding mean free path at a given epoch can be finally computed assuming Bohm diffusion, as in Eq. (42), $\lambda_{\rm CR}=\sqrt{t_{B}D_{B}}$. The magnetic field is given by Eq. (51) with $r$ representing now the mean free path. This results in

[TABLE]

If this is the case, it would imply that CRs heating is far from uniform; rather, it is highly patchy and clustered around the smallest star forming halos.

In practice, a number of neglected effects might reduce the efficiency of CR self-confinement. These are: (a) the presence of neutrals outside the fully ionized bubbles can quickly damp the waves generated through the CR streaming instability; (b) the B-field equipartition value in Eq. (51) might not be attained due to an inefficient CR-magnetic energy density conversion. Moreover, a sufficiently strong intergalactic B-field and/or a smaller galactic CR luminosity may result in a magnetic-CR energy density ratio that is too large for the development of the instability in the non-resonant regime.

We note that observations of the redshifted HI 21 cm line from these high redshifts would be very sensitive to the morphology of the pre-heated, neutral regions. We thus expect that the clustered heating scenario leaves unique imprints in the power spectrum of such radiation. Additionally, the analysis of the power spectrum should allow to discriminate between the case in which the heating source are X-rays or CRs. Finally, we could also gain precious information about the strength and structure of early intergalactic magnetic fields and the efficiency of CR acceleration by the first SNe. All these aspects are very hard to investigate with any other mean.

5 Conclusions

In this work we have shown that CRs can influence the temperature and ionization fraction of the IGM using a self-consistent model for galaxy formation and cosmic reionization. The model was designed to reproduce the observed SFR at redshift $z\lesssim 10$, and a cosmic history consistent the latest PLANCK results. Such data constrain the conversion efficiency of gas into stars ($f_{*}=0.04$) and the population-averaged escape fraction of ionizing photons into the IGM ($f_{\rm esc}=0.01$). From the supernova rate evolution given by the model we further derived the CR energy density.

CRs with energies $<1$ MeV lose energy predominantly by ionizations at redshifts $z>10$, and by Coulomb scatterings at lower redshifts (see Fig. 4). The energy lost via Coulomb collisions goes directly to heat, increasing the IGM temperature at $z\sim 10$ above the standard adiabatic thermal evolution $\Delta T\sim 10-200$ K, depending on the slope of the CR injection spectrum in the range $2<\alpha<2.5$, and on the transport efficiency of CRs out of the first star-forming structures. Such increase is comparable or higher than that produced by two other popular heating mechanisms, i.e. X-rays (Mesinger et al., 2013) and dark matter annihilations (Evoli et al., 2014). Plasma instabilities induced by blazar TeV photons can additionally heat up the IGM above the temperature induced by photo-heating (Chang et al., 2012; Chang et al., 2016). Compared to our results, this mechanism is relevant at lower redshifts, $z\lesssim 6$, and provides a more uniform background. The signal yielded by such contribution would therefore be easily distinguishable from the one expected from CRs.

Our model for the CR injection and transport in the IGM is based on two commonly accepted assumptions: (1) CRs escape from star-forming structures on a timescale much shorter than the energy loss timescale in the ISM, (2) CRs provide a spatially uniform energy density floor, i.e. a “background”.

The first assumption is certainly valid for Milky Way protons with $E>100$ MeV. Whether it holds for high-$z$ galaxies depends on poorly known quantities, such as the turbulent magnetic fields in these objects. On general grounds, however, weaker magnetic fields should correspond to a larger diffusion coefficient and a smaller size of the magnetic halo. Both factors lead to a shorter diffusion timescale than in the Galaxy.

The second assumption has been carefully investigated in § 4.3. We showed that CRs escaping from galaxies trigger streaming instabilities eventually amplifying the seed magnetic field up to equipartition. Under the most optimistic conditions for the development of the instabilities, such self-generated magnetic field might efficiently confine GeV particles around haloes for a time largely exceeding the Hubble time at $z\sim 20$. If true, this strongly clustered emission is expected to leave a specific imprint on the 21cm line power spectrum. Such detection would allow for the first time to study the structure and strength of magnetic fields in the Dark Ages. However, a testable prediction of this CR heating signature requires a more detailed model for the ejection and propagation of CRs in the pre-ionized bubbles and will be investigated in a forthcoming work.

Acknowledgments

N.L. thanks GSSI in L’Aquila for the warm hospitality during the preparation of the paper. We thank A. Mesinger for useful discussions. This work was partially supported by the ”Helmholtz Alliance for Astroparticle Physics (HAP)” funded by the Initiative and Networking Fund of the Helmholtz Association, and by the Deutsche Forschungsgemeinschaft (DFG) through the Collaborative Research Centre SFB 676 ”Particles, Strings and the Early Universe”.