Hadronic Interactions of Energetic Charged Particles in Protogalactic Outflow Environments and Implications for the Early Evolution of Galaxies

TL;DR

This paper models how energetic hadronic particles in protogalactic outflows heat the surrounding medium through pion production, affecting galaxy evolution and star formation, with implications for high-redshift galaxy observations.

Contribution

It introduces a detailed analysis of cosmic ray heating in protogalactic outflows, comparing advection and diffusion regimes and applying the model to high-redshift galaxy data.

Findings

Peak heating rate of 10^{-26} erg cm^{-3} s^{-1} at the outflow base in diffusion regime

Heating extends up to tens of kpc in advection regime

Approximately 10% of cosmic rays escape via advection

Abstract

We investigate the interactions of energetic hadronic particles with the media in outflows from star-forming protogalaxies. These particles undergo pion-producing interactions which can drive a heating effect in the outflow, while those advected by the outflow also transport energy beyond the galaxy, heating the circumgalactic medium. We investigate how this process evolves over the length of the outflow and calculate the corresponding heating rates in advection-dominated and diffusion-dominated cosmic ray transport regimes. In a purely diffusive transport scenario, we find the peak heating rate reaches $10^{-26}\;\! {\rm erg~cm}^{-3}\;\! {\rm s}^{-1}$ at the base of the outflow where the wind is driven by core-collapse supernovae at an event rate of 0.1 $\text{yr}^{-1}$, but does not extend beyond 2 kpc. In the advection limit, the peak heating rate is reduced to $10^{-28}\;\! {\rm erg~cm}^{-3}\;\! {\rm s}^{-1}$, but its extent can reach to tens of kpc. Around 10% of the cosmic rays injected into the system can escape by advection with the outflow wind, while the remaining cosmic rays deliver an important interstellar heating effect. We apply our cosmic ray heating model to the recent observation of the high-redshift galaxy MACS1149-JD1 and show that it could account for the quenching of a previous starburst inferred from spectroscopic observations. Re-ignition of later star-formation may be caused by the presence of filamentary circumgalactic inflows which are reinstated after cosmic ray heating has subsided.

Equations

$$\frac{1}{h^2}\frac{\mathrm{d}}{\mathrm{d}h}\left(\rho v{h}^2\right)$$

$$v\rho \frac{\mathrm{d}v}{\mathrm{d}h}$$

$$\frac{1}{h^2}\frac{\mathrm{d}}{\mathrm{d}h}\left[\rho v{h}^2\left(\frac{1}{2}{v}^2+\frac{{\gamma }_{\mathrm{g}}}{{\gamma }_{\mathrm{g}}-1}\frac{P}{\rho }\right)\right]$$

$$\frac{1}{h^2}\frac{\mathrm{d}}{\mathrm{d}h}\left[\frac{{\gamma }_{\mathrm{C}}}{{\gamma }_{\mathrm{C}}-1}{P}_{\mathrm{C}}{h}^2(v+v_{\mathrm{A}})\right]$$

$$B(h,{\mathcal{R}}_{\mathrm{SN}})=B_0{\left(\frac{{\mathcal{R}}_{\mathrm{SN}}}{{\mathcal{R}}_{\mathrm{SN},0}}\right)}^{1/2}{\left(1+{\left[\frac{h}{h_{\mathrm{B}}}\right]}^2\right)}^{-1},$$

$$I$$

$$\Phi (h)=-\frac{3\mathrm{G}{M}_{\mathrm{tot}}}{h}\ln \left\{h+\frac{h}{R_{\mathrm{s}}}\right\},$$

$$q=\mathcal{P}\left[{\mathcal{R}}_{\mathrm{SN}}{M}_{\mathrm{inj}}\right]=\mathcal{P}\left[\frac{{\alpha }_{\ast }{\mathcal{R}}_{\mathrm{SF}}{M}_{\mathrm{inj}}}{{\overline{M}}_{\mathrm{SN}}}\right],$$

$${\alpha }_{\ast }=\frac{{\int }_{M_{\mathrm{S}{\mathrm{N}}^{\ast }}}^{M_{\max }}\mathrm{d}M{M}^{-{\rm Y}}}{{\int }_{M_{\mathrm{cut}}}^{M_{\max }}\mathrm{d}M{M}^{-{\rm Y}}}.$$

$${\dot{ϵ}}_{\mathrm{th}}=\mathcal{Q}\left[\nu \xi {\mathcal{R}}_{\mathrm{SN}}{E}_{\mathrm{SN}}\right],$$

$${\dot{ϵ}}_{\mathrm{CR}}=\zeta \frac{{\Omega }_{\mathrm{A}}}{4\pi }\left[\nu \xi {\mathcal{R}}_{\mathrm{SN}}{E}_{\mathrm{SN}}\right]$$

$${\Omega }_{\mathrm{A}}=2\pi \left[1-\cos \left(\frac{\theta }{2}\right)\right].$$

$$\dot{ϵ}={\eta }_{\mathrm{SNe}}\left[\nu \xi {\mathcal{R}}_{\mathrm{SN}}{E}_{\mathrm{SN}}\right]={\eta }_{\mathrm{SNe}}\left[\frac{\nu \xi {\alpha }_{\ast }{\mathcal{R}}_{\mathrm{SF}}{E}_{\mathrm{SN}}}{{\overline{M}}_{\mathrm{SN}}}\right],$$

$${\eta }_{\mathrm{SNe}}=\left(\mathcal{Q}+\zeta \frac{{\Omega }_{\mathrm{A}}}{4\pi }\right).$$

$$Q=f\dot{ϵ}=f{\eta }_{\mathrm{SNe}}\left[\frac{\nu \xi {\alpha }_{\ast }{\mathcal{R}}_{\mathrm{SF}}{E}_{\mathrm{SN}}}{{\overline{M}}_{\mathrm{SN}}}\right],$$

$$\rho v\left[v\frac{\mathrm{d}v}{\mathrm{d}h}+\frac{{\gamma }_{\mathrm{g}}}{{\gamma }_{\mathrm{g}}-1}\frac{\mathrm{d}}{\mathrm{d}h}\left(\frac{P}{\rho }\right)\right]=-\rho v\frac{\mathrm{d}\Phi }{\mathrm{d}h}-(v_{\mathrm{A}}+v)\frac{\mathrm{d}{P}_{\mathrm{C}}}{\mathrm{d}h}.$$

$$\frac{\mathrm{d}P}{\mathrm{d}h}={\gamma }_{\mathrm{g}}\frac{P}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}h}-({\gamma }_{\mathrm{g}}-1)\left(\frac{v_{\mathrm{A}}}{v}\right)\frac{\mathrm{d}{P}_{\mathrm{C}}}{\mathrm{d}h}.$$

$$\frac{1}{v_{\mathrm{A}}}\frac{\mathrm{d}{v}_{\mathrm{A}}}{\mathrm{d}h}=-\frac{2h}{h_{\mathrm{B}}^2+h^2}-\frac{1}{2\rho }\frac{\mathrm{d}\rho }{\mathrm{d}h}.$$

$$\frac{1}{v_{\mathrm{A}}}\frac{\mathrm{d}{v}_{\mathrm{A}}}{\mathrm{d}h}\approx -\frac{2}{h}-\frac{1}{2\rho }\frac{\mathrm{d}\rho }{\mathrm{d}h},$$

$$\frac{1}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}h}+\frac{1}{v}\frac{\mathrm{d}v}{\mathrm{d}h}+\frac{2}{h}=0.$$

$$\frac{\mathrm{d}{P}_{\mathrm{C}}}{\mathrm{d}h}=\frac{{\gamma }_{\mathrm{C}}{P}_{\mathrm{C}}}{\rho }\left(\frac{v+v_{\mathrm{A}}/2}{v+v_{\mathrm{A}}}\right)\frac{\mathrm{d}\rho }{\mathrm{d}h},$$

$$\frac{\mathrm{d}P}{\mathrm{d}h}=\left\{{\gamma }_{\mathrm{g}}P-{\gamma }_{\mathrm{C}}({\gamma }_{\mathrm{g}}-1)P_{\mathrm{C}}\left[\frac{v+v_{\mathrm{A}}/2}{v+v_{\mathrm{A}}}\right]\left(\frac{v_{\mathrm{A}}}{v}\right)\right\}\frac{1}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}h}.$$

$$\rho v\frac{\mathrm{d}v}{\mathrm{d}h}+c_{\ast }^2\frac{\mathrm{d}\rho }{\mathrm{d}h}=-\rho \frac{\mathrm{d}\Phi }{\mathrm{d}h}$$

$$c_{\ast }^2=\frac{{\gamma }_{\mathrm{g}}P}{\rho }-\frac{{\gamma }_{\mathrm{C}}{P}_{\mathrm{C}}}{\rho }\left\{\left(\frac{v+v_{\mathrm{A}}/2}{v+v_{\mathrm{A}}}\right)\left({\gamma }_{\mathrm{C}}-{\gamma }_{\mathrm{C}}({\gamma }_{\mathrm{g}}-1)\left[\frac{v_{\mathrm{A}}}{v}\right]\right)\right\}.$$

$$\frac{\mathrm{d}v}{\mathrm{d}h}=\frac{2v}{h}\frac{\left(c_{\ast }^2-\frac{h}{2}\frac{\mathrm{d}\Phi }{\mathrm{d}h}\right)}{v^2-c_{\ast }^2}.$$

$$c_{\ast }^2-\frac{h}{2}\frac{\mathrm{d}\Phi }{\mathrm{d}h}=0.$$

$$c_{\ast }^2$$

$$=\frac{\mathrm{G}M(h_{\ast })}{2h_{\ast }}=\frac{F^2{V}_{\mathrm{c},\mathrm{vir}}^2}{2},$$

$$\mathrm{p}+\mathrm{p}\to \{\begin{array}{ll}& \mathrm{p}{\mathrm{\Delta }}^{+}\to \{\begin{array}{l}\mathrm{p}\mathrm{p}{\pi }^0{\xi }_0({\pi }^0){\xi }_{\pm }({\pi }^{+}{\pi }^{-})\\ \mathrm{p}\mathrm{p}{\pi }^{+}{\pi }^{-}{\xi }_0({\pi }^0){\xi }_{\pm }({\pi }^{+}{\pi }^{-})\\ \mathrm{p}\mathrm{n}{\pi }^{+}{\xi }_0({\pi }^0){\xi }_{\pm }({\pi }^{+}{\pi }^{-})\end{array}\\ & \mathrm{n}{\mathrm{\Delta }}^{++}\to \{\begin{array}{l}\mathrm{n}\mathrm{p}{\pi }^{+}{\xi }_0({\pi }^0){\xi }_{\pm }({\pi }^{+}{\pi }^{-})\\ \mathrm{n}\mathrm{n}2{\pi }^{+}{\xi }_0({\pi }^0){\xi }_{\pm }({\pi }^{+}{\pi }^{-})\end{array}\end{array},$$

$${\pi }^0$$

$${\pi }^{+}$$

$${\pi }^{-}$$

$${\dot{E}}_{\mathrm{rad}}=\frac{4}{3}{\sigma }_{\mathrm{T}}\mathrm{c}{\left(\frac{E_{\mathrm{e}}}{m_{\mathrm{e}}{\mathrm{c}}^2}\right)}^2{U}_{\mathrm{i}}$$

Full text

††pubyear: 2018††pagerange: Hadronic Interactions of Energetic Charged Particles in Protogalactic Outflow Environments and Implications for the Early Evolution of Galaxies–A.2

1 Introduction

Advection and diffusion are the two main mechanisms for the transportation of high-energy charged hadronic (cosmic ray) particles in galactic environments. These are evident in the Galactic interstellar medium (ISM) (see Schlickeiser, 2002; Strong et al., 2007; Gaggero et al., 2015b; Korsmeier & Cuoco, 2016; Yuan et al., 2017), and within the solar system (see Jokipii, 1966; Orlando & Strong, 2008; Abdo et al., 2011; Potgieter, 2013; Chhiber et al., 2017). The interplay between the two processes is determined by the extent of the bulk flows in the carrying medium (e.g. Dorfi, E. A. & Breitschwerdt, D., 2012; Uhlig et al., 2012; Heesen et al., 2016; Taylor & Giacinti, 2017; Farber et al., 2018), the structures and strengths of the local magnetic field (e.g. Parker, 1964; Jokipii, 1966; Berezinskii et al., 1990; Schlickeiser, 2002; Alvarez-Muniz et al., 2002; Aharonian et al., 2012; Gaggero, 2012; Snodin et al., 2016) and the amount of turbulence present in the system (e.g. Berezinskii et al., 1990; Schlickeiser, 2002; Candia & Roulet, 2004; Gaggero, 2012; Snodin et al., 2016).

The general consensus is that cosmic rays (CRs) are accelerated to high energies in violent environments, e.g. SN explosions, gamma-ray bursts, large-scale shocks in the ISM or galactic outflows, AGN jets, galaxy clusters and compact objects such as fast spinning neutron stars and accreting black holes (see Berezinsky et al., 2006; Brunetti et al., 2007; Pfrommer et al., 2007a; Dar & de Rújula, 2008; Reynoso, M. M. et al., 2011; Kotera & Olinto, 2011). Fermi processes (Fermi, 1949) have been suggested as viable mechanisms by which low-energy charged particles can be accelerated to attain relativistic energies. In astrophysical systems this might arise in shocks, such as those resulting from SN explosions. Systems such as starburst galaxies, which have frequent SN events, are therefore expected be abundant in energetic CRs (Karlsson, 2008; Lacki et al., 2011; Lacki & Thompson, 2012; Wang & Fields, 2014; Farber et al., 2018). Likewise protogalaxies, which have vibrant star forming activity and hence high SN event rates, should also be abundant in CRs. In a similar way to the shocks in the ISM generated by SN explosions, large-scale shocks in the intergalactic medium (IGM) and intracluster medium (ICM) can also be accelerators of CRs. There is evidence that energetic CRs are an important ingredient in galaxy clusters (Takami & Sato, 2008; Brunetti & Jones, 2014), providing pressure support to clusters’ structure. These CRs may play an important role in regulating the energy budget of the ICM through radiative losses, energy transportation and hadronic interactions.

While we have some understanding of the effects of CRs on the thermal and dynamical properties of the ISM and IGM in nearby astrophysical systems, our knowledge of the impacts of CR particles on the formation and evolution of structures on scales of galaxies or larger is very limited. The importance of CRs in protogalactic environments has gradually drawn more attention (e.g. Giammanco & Beckman, 2005; Stecker et al., 2006; Valdés et al., 2010; Sazonov & Sunyaev, 2015; Bartos & Marka, 2015; Leite et al., 2017; Owen et al., 2018). In particular, there are studies showing that CR heating of the ISM could lead to the distortion of the stellar initial mass function (IMF) and even quench the star formation process entirely (see Pfrommer et al., 2007b; Chen et al., 2016). CRs can also drive the large-scale galactic outflows (see Socrates et al., 2008; Weiner et al., 2009), which transfer energy and chemically enriched material into intergalactic space. The resulting pre-heating of the IGM in turn alters cosmological structural formation processes, and there are arguments that CRs might contribute to a certain degree of cosmological reionisation (see Nath & Biermann, 1993; Sazonov & Sunyaev, 2015; Leite et al., 2017).

On sub-galactic scales, the production of CRs is often attributed to supernovae (SNe), compact objects (such as spinning neutron stars) and accretion-powered sources, which are consequential of stellar evolution and hence star formation processes. However, the delivery of CR energy across a galaxy depends on strength and structure of the galactic magnetic field which, in turn, depends on the field evolution and hence the star-forming processes. SN explosions are energetic events. On the one hand, SN explosions would drive a large-scale galactic wind (Chevalier & Clegg, 1985; Socrates et al., 2008; Weiner et al., 2009), but on the other hand, they inject enormous amount of mechanical energy into the ISM within the galaxy which fuels the development of ISM turbulence (Dib et al., 2006; Joung et al., 2009; Gent et al., 2013; Martizzi et al., 2015, 2016). SN explosions also help the magnetisation of the entire galaxy (Zweibel & Heiles, 1997; Zweibel, 2003; Beck et al., 2012; Schober et al., 2013; Lacki & Beck, 2013).

In a magnetised medium with strong turbulence but no large-scale bulk flows, CR transport would be dominated by diffusion and this can lead to their effective containment within the host galaxy (see Owen et al., 2018) with their subsequent energy deposition into the media being regulated by hadronic, pion-producing interactions occurring above a threshold energy of 0.28 GeV (Kafexhiu et al., 2014). However, in the presence of a flow, CRs entangled into a magnetised medium can be advected along. This advection process would happen in large-scale galactic winds, where diffusion still takes place, but on time scales far longer than the flow time scale (see Berezinskii et al., 1990; Schlickeiser, 2002; Aharonian et al., 2012; Heesen et al., 2016). As such, CRs can be advected into intergalactic space causing heating of the circumgalactic medium. Imaging observations of nearby starburst galaxies have shown complex structural morphologies in which winds and outflows are “collimated" in a cone-like structure while the gases and stars beneath retain a planar galactic disk-like structure. Such structural complexity implies the coexistence of CR diffusion and advection — while in some regions the two processes would have comparable partitions in facilitating energy transport, in other regions one of them would dominate.

Here we further investigate the contribution of CRs to ISM and IGM heating via hadronic processes with a focus on the effects of CRs by galactic wind outflows. We organise the paper as follows. In § 2 we discuss the properties of galactic outflows driven by SNe and CRs and the outflow model used for our investigation of CR heating. In § 3, we present the formulation for CR propagation in the diffusion and advection dominated regimes. We also discuss the relevant mechanisms regulating the energy budget of CRs advected by bulk flows and the hadronic processes by which the CR energy is deposited into the ISM and/or IGM. In § 4, we show the results of our calculations of CR heating in protogalactic and outflow environments and demonstrate how the heating effect depends on model parameters, and discuss the astrophysical implications. An application to explain the inferred star-formation behaviour of the high-redshift galaxy MACS1149-JD1 (Hashimoto et al., 2018) is presented in § 5 and conclusions are given in § 6.

Our calculations assume that CRs are energetic protons. This assumption is based on the idea that CRs are produced and accelerated (see Berezinsky et al., 2006; Kotera & Olinto, 2011) in a similar manner in the distant Universe as they are in the nearby Universe, and that a substantial fraction of CRs detected on Earth (from the nearby Universe) are protons (Abbasi et al., 2010). This allows us to ignore the composition evolution of CRs as a first approximation. We also do not consider CR primary electrons explicitly as these charged, low-mass leptonic particles have considerably higher radiative loss rates than charged hadrons. This means that they would not be a major contributor to the global energy transportation picture. Hereafter, unless it is necessary, we do not differentiate between CR particle species, and CR protons are referred to as CRs.

2 Galactic Outflows

2.1 Observational aspects and phenomenology

Galactic-scale outflows have been observed in star-forming galaxies nearby, e.g. Arp 220 (Lockhart et al., 2015), and in the distant Universe (Frye et al., 2002; Ajiki et al., 2002; Benítez et al., 2002; Rupke et al., 2005a, b; Bordoloi et al., 2011; Arribas, S. et al., 2014). In active star-forming regions, the proximity of SNe allow the confluence of gas flows induced by the SN explosions to develop into a larger-scale wind. The build-up of these confluent winds eventually erupts as a large-scale galactic outflow, usually with a bi-conical structure along the minor axis of the host (Veilleux et al., 2005). The opening angles of the outflow cones are broad, of tens of degrees (Heckman et al., 1990; Veilleux et al., 2005), and their values vary between galaxies, e.g. around [math]$-$[math] in NGC 253 (Strickland et al., 2000; Bolatto et al., 2013) and approximately [math] in M82 (Heckman et al., 1990; Walter et al., 2002). Observations have shown that the hot X-ray emitting gas in an outflow can reach up to 3 kpc (Strickland et al., 2000; Cecil et al., 2002b, a), and the entire outflow structure could extend up to tens of kpc (see Veilleux et al., 2005; Bland-Hawthorn et al., 2007; Bordoloi et al., 2011; Martin et al., 2013; Rubin et al., 2014; Bordoloi et al., 2016). Galactic outflows are inhomogeneous, multicomponent, multiphase media, with warm, partially ionised gases intermingled with hot ionised bubbles and cooler, denser less ionised gas or neutral clumps. The outflow velocity has been measured from a few hundred km s*-1* in most cases, rising to a few thousand km s*-1* in a few extreme systems (Rupke et al., 2005b; Cecil et al., 2002a; Rubin et al., 2014), and the total mechanical power in the flow is estimated to be up to levels as high as $10^{43}$ erg s*-1* (see Cecil et al., 2002a). Galactic outflows play an active role in injecting mechanical energy into intergalactic space and in distributing chemically enriched matter into the environment (Aguirre et al., 2001b, a; Martin et al., 2002; Rupke et al., 2003; Adelberger et al., 2003; Aguirre et al., 2005; Bertone et al., 2005). There are also arguments that galactic outflows carry nG-strength magnetic fields from the within the galaxy to the surrounding IGM and ICM, where the seed fields would then be amplified by turbulence, dynamo mechanisms and/or cosmic-scale shears to the $\mu$G-levels inferred from observations (De Young, 1992; Goldshmidt & Rephaeli, 1993; Dolag et al., 1999; Dolag, K. et al., 2002; Bertone et al., 2006; Vazza et al., 2018). Of most relevance to this work, galactic outflows are efficient vehicles to transport CRs and the energy they carry across their source galaxy and to significant distances away from it (see Heesen et al., 2016).

2.2 Outflow Wind Structure

In our calculations, we adopt a working model that sufficiently captures the most essential microphysics and the associated global physics and astrophysics of the system. We focus on the redistribution of energy through the advective and diffusive transport of CRs and investigate the relative efficiency of CR heating in galactic outflows and the surrounding IGM between these two modes of CR transportation. Complexities such as the fine substructure of outflows, the multi-phase nature of the flow material and the re-acceleration of CR particles within the flow are worth separate further investigations in their own right, and so are not considered in detail in the present study, instead being left to future follow-up work.

Galactic outflows can be powered by different mechanisms. Most early models invoke thermally and/or SN-driven mechanisms (Larson, 1974; Chevalier & Clegg, 1985; Dekel & Silk, 1986; Nath & Trentham, 1997; Efstathiou, 2000; Madau et al., 2001; Furlanetto & Loeb, 2003; Scannapieco, 2005; Samui et al., 2008) while more recent studies have also considered radiatively-driven outflows (Dijkstra & Loeb, 2008; Nath & Silk, 2009; Thompson et al., 2015). CRs have also been regarded as a means by which galactic outflows may be driven: indeed, a CR driving mechanism would offer a good explanation of the observed soft X-ray diffuse emission from the Milky Way (Everett et al., 2008). At high redshift, actively star-forming galaxies would be abundant in CRs and thus are a clear candidate driver of outflows (see also Samui et al., 2010; Uhlig et al., 2012), while the high altitude winds which have the ability to inject CRs far into the circumgalactic medium are thought to be powered by CRs (Jacob et al., 2018). CRs influence the density and structure of the flow compared to other driving mechanisms (Girichidis et al., 2018) and also lose some of their energy in driving the outflow (e.g. Samui et al., 2010; Uhlig et al., 2012). These factors modify the heating effect that they are able to deliver when interacting with the wind fluid via hadronic interactions when compared to their role in winds driven by other mechanisms and, as such, mean that CRs must be self-consistently included in the modelling of the wind structure and dynamics. At high-redshift, the impact of smaller galaxies, of mass around $10^{9}~{}\text{M}_{\odot}$, on their environment is argued to be more important than more massive galaxies (see, e.g. Samui et al., 2008, 2009). As such we focus this study on these smaller systems, for which the impact will be greatest.

The direction of emergence of a galactic outflow is governed by the ‘path of least resistance’. In disk galaxies, a bi-conical flow pattern above and below the galactic plane is generally observed (see Veilleux et al., 2005), and this geometrical shape is expected because a spherically expanding wind from an actively star forming region at the core of the galaxy would be less obstructed by the the upper and lower edge of the galaxy that it encounters than the galactic plane. Emerging flows from spherical or near-spherical elliptical galaxies would have a less well-defined morphological pattern. An isotropic spherical outflow could arise if the outflowing wind from the galactic core region encounters all edges of the galaxy at a similar time, and if it is faced with similar inflowing pressures and resistances in all directions. We consider disk galaxies with bi-conical outflows. A schematic of the outflow model is shown in Fig. 1 (see also outflow wind morphology in e.g. Strickland et al., 2000; Ohyama et al., 2002; Veilleux et al., 2005; Cooper et al., 2008), where two distinct ‘zones’ are noted: Zone A is within the outflow cone where, we assume, CR transport is dominated by advection; Zone B is outside of the cone but within the galactic ISM and is the region in which CR transport is predominantly diffusive. We explore a range of opening angles, $\theta$, between [math] and [math] covering a similar range of values to a large subset of those observed in nearby starburst galaxy outflows. A reference value of [math] is chosen (if not otherwise specified) as a working representation111Hydrodynamical simulations suggest that, rather than remaining uniform throughout the extent of an outflow, the opening angles start at a low value of [math]$-$[math] near their base and then diverge to [math]$-$[math] well above and below the galactic plane (Mac Low & Ferrara, 1999; Martel & Shapiro, 2001; Pieri et al., 2007; Bordoloi et al., 2016), but this finer substructure is not accounted for in our model - instead we choose an opening angle which reflects that of the wider angle of the main part of the outflow..

Ipavich (1975) considered a 1-D numerical magnetohydrodynamic model for a CR-powered wind emerging from a galaxy with a point-like mass distribution. The model is parametrised with energy and mass injection, presumably provided by the SNe resulting from the starburst activity. A spherical geometry is assumed, with the wind emerging radially from a small active star-forming region enveloping the galactic core. Solving the associated magnetohydrodynamic equations yield several valid solutions depending on the boundary conditions adopted at a so-called critical point, at which the flow becomes supersonic. One such solution is that of an outflow wind with an asymptotic velocity for the outflow at distances sufficiently far from the star-forming galactic core. This idea was further developed by several authors since then, including  Breitschwerdt et al. (1987, 1991, 1993); Everett et al. (2008); Bustard et al. (2016); Recchia et al. (2016) and Samui et al. (2010). The latter of these accounts for a CR-driven outflow in the presence of an NFW density profile, with the intention of application to high-redshift starburst galaxies. We largely follow the Samui et al. (2010) outflow model here and use it to compute the density profile and other relevant conditions from which the advection of CRs and their subsequent hadronic interaction induced heating effect can be determined.

2.3 Outflow Model

We denote $h$ as a coordinate variable along a flow streamline and, in a spherical symmetric geometry, $h$ is the radial distance from the galactic core. We model the CRs and wind fluid as two separate but interacting components of an outflow wind in which the CR component has negligible mass density but non-negligible energy density. In the Samui et al. (2010) model the galactic outflow wind is a conic section of a spherical flow, with an asymptotic velocity arising at a sufficiently large distance from the galactic core region. This can be found by considering the steady-state spherically symmetric form of the fluid and CR equations (Ipavich, 1975; Breitschwerdt et al., 1991):

[TABLE]

where (1) is the mass continuity equation, (2) is the momentum equation and (3) is the energy equation for the wind fluid, while (4) is the evolution equation for the CR fluid component of the wind. $\rho$ is the density of the wind fluid, $v$ is the wind velocity, $P$ is the pressure of the wind fluid (gas), $P_{\rm C}$ is the CR pressure, $\Phi$ is the gravitational potential, $\gamma_{\rm g}=5/3$ is the adiabatic index for the gas component, and $\gamma_{\rm C}=4/3$ is the adiabatic index for the relativistic CR component. We specify the total mass injection rate into the wind as $\rho vh^{2}=q$, from equation 1, with $q$ as a mass injection rate due to SN mass-loading of the wind (see equation 8). $I$ is an energy exchange term between the CRs and baryonic wind fluid (see equation 6).

We adopt a magnetic field strength and morphology along the outflow cone according to:

[TABLE]

where $B_{\rm 0}=5\mu\text{G}$, $\mathcal{R}_{\rm SN,0}=0.1~{}\text{yr}^{-1}$, and where $h_{\rm B}$ is introduced as a characteristic scale over which the magnetic field does not strongly vary within the host galaxy of the outflow. Physically, the variation of the magnetic field with $h$ would only be expected in regions of the model that are well within the outflow cone. In regions which may better be regarded as interstellar environments, the magnetic field would vary less substantially with height. We find that a choice of $h_{\rm B}=1.5~{}\text{kpc}$ yields a relatively uniform magnetic field within a 0.5 kpc starburst region, falling only by around 10% from its peak value. Beyond this, $B(h)$ reverts to an inverse-square law behaviour with $h$ thus ensuring the conservation of magnetic flux along the outflow. The dependence of magnetic field strength at the base of the outflow (within the ISM of the host) on the square root of SN-rate follows from Schober et al. (2013), which models the development of magnetic fields in young starburst galaxies via turbulent dynamo amplification.

Our choice of $B_{\rm 0}$ is reflective of interstellar environments, where energy densities of CRs at the peak of their spectrum are comparable to that of the magnetic field. As CRs gyrate and stream along the magnetic field lines at speeds faster than the Alfvén velocity $v_{\rm A}=B(h)/\sqrt{4\pi\rho}$, they amplify interstellar Alfvén waves which have wavelengths comparable to the gyro-radii of the streaming CRs (Wentzel, 1968; Kulsrud & Pearce, 1969; Kulsrud & Cesarsky, 1971). In this process, known as the streaming instability (Wentzel, 1968; Kulsrud & Pearce, 1969), this leads to a resonant scattering effect, which slows the CRs and transfers momentum and energy from the CRs to the ambient medium after dampening of the waves at a rate given by $|v_{\rm A}\cdot\nabla P_{\rm C}|$ (e.g. Wentzel, 1971; Ipavich, 1975; Breitschwerdt et al., 1991; Uhlig et al., 2012). Further losses by the CRs to the wind fluid result from the work done by the CR pressure gradient in a bulk wind velocity $v$, arising at a rate of $|v\cdot\nabla P_{\rm C}|$ (Samui et al., 2010). Together, this allows us to define $I$ as the total energy exchange term between the baryonic wind fluid and the energetic CR component, given by:

[TABLE]

(Samui et al., 2010), where the minus sign is due to the energy exchange resulting in a loss by the CR component and a gain by the wind fluid.

Samui et al. (2010) solve this system of equations when adopting an NFW (Navarro–Frenk–White, Navarro et al., 1996) gravitational potential of the form

[TABLE]

for $M_{\rm tot}$ as the total galaxy mass and where $R_{\rm S}$ is the scale height, being the ratio of the virial radius of the galaxy and the concentration parameter, $R_{\rm vir}/c_{\rm g}$. This potential is relevant to galaxies like that which we also wish to model here. In the system of equations, $q$ is the volumetric mass injection (which is non-zero only within the starburst region, i.e. $h<h_{\rm inj}$). This may be quantified in terms of the SN event rate, $\mathcal{R}_{\rm SN}$, and the mass ejecta $M_{\rm inj}$ resulting from a SN explosion to estimate $\dot{M}$. The level of mass injection per event varies with SN types. Type II SNe offer a characteristic mass of around $10\;\!{\rm M}_{\odot}$, and Type I b/c of a few ${\rm M}_{\odot}$ (see, e.g. Branch, 2010; Perets et al., 2010). Note that lower mass stars take a longer time to complete their life cycles and so, at high-redshifts, only the very high-mass stars would have enough time to evolve to the SN stage within the host galaxy’s evolutionary timescale. Moreover, a low metallicity environment would yield a more top-heavy initial stellar mass distribution (e.g. Abel et al., 2002; Bromm et al., 2002; Bastian et al., 2010; Gargiulo et al., 2015), skewing the progenitor masses to favour the core-collapse SN channels more typical of massive stars. Thus, core-collapse SNe and hypernovae, with progenitors of masses ${{M}}_{{\rm SN}}\sim 8.5~{}{\rm M}_{\odot}$ or higher (see, e.g. Smartt et al., 2009; Smartt, 2009; Smith et al., 2011) would occur frequently in star-bursting protogalaxies. These core-collapse SN and hypernovae are extremely energetic, with $E_{\rm SN}\approx 10^{53}\;\!{\rm erg}$ per event (Smartt, 2009).

The mass injection rate may be parametrised as

[TABLE]

where ${\cal{R}}_{\rm SF}$ is the star formation rate and $\overline{M}_{\rm SN}$ is mean SN progenitor mass. The parameter $\mathcal{P}$ is the mass-loading factor, which is a scaling factor specifying the mass loaded into the wind for a given mass ejected from the progenitor star in a SN event. Although ${\mathcal{P}}$ could have a value above 1 and our limited knowledge of the ISM environment and the properties of SNe in protogalaxies prevents us from deriving a strong constraint for appropriate values of this parameter  (see e.g. Martin et al., 2002, which gives mass loading fractions of 10 and above in NGC 1569, among others), we conservatively adopt that ${\mathcal{P}}=0.1$. The parameter $\alpha_{*}$ is the fraction of stars that yield Type II SNe (and hypernovae), which can be estimated as

[TABLE]

As a conservative estimate, we adopt a Salpeter IMF index of $\Upsilon=2.35$222There is some evidence of a redshift-dependent IMF (Lacey et al., 2008; Davé, 2008; van Dokkum, 2008; Hayward et al., 2013), but this remains an open discussion (see e.g Bastian et al., 2010; Cen, 2010, for reviews). Given the lower metallicities, higher cosmic microwave background temperature which could influence molecular-cloud collapse, and the tendency for star-formation at high redshifts to arise in the ‘burst’ mode rather than more gradually (see Lacey et al., 2008), the IMF at high redshift may be more top-heavy – favouring the production of the very massive stars compared to the IMF observed in the current epoch, or the Salpeter IMF. A Top-heavy IMF has been claimed for some nearby starburst systems (Weidner et al., 2011; Bekki & Meurer, 2013; Chabrier et al., 2014), e.g. in M82 (Rieke et al., 1993; McCrady et al., 2003) and NGC 3603 (Harayama et al., 2008), and even for Galactic centre clusters (Stolte et al., 2005; Maness et al., 2007). In these cases, the Salpeter IMF index of 2.35 underestimates the number of high-mass stars and hence the SN events.. We set the the maximum stellar mass333Note this is at the upper end of likely progenitor masses to ensure that our calculation is conservative. In reality a greater proportion of massive stars are more likely to arise in the protogalactic environments. which could reasonably yield a SN explosion to be $M_{\rm max}=50\;\!{\rm M}_{\odot}$ (Fryer, 1999; Heger et al., 2003), the stellar mass cut-off $M_{\rm cut}\approx 1\;\!{\rm M}_{\odot}$ and the minimum mass required for a core-collapse SN event $M_{\rm SN^{*}}=8.5\;\!{\rm M}_{\odot}$ (Smith et al., 2011; Eldridge et al., 2013). This gives $\alpha_{*}\approx 0.05$, implying a scaling relation $\mathcal{R}_{\rm SF}\approx 160\;\!{\rm M}_{\odot}~{}{\mathcal{R}}_{\rm SN}$ between the star-formation rate $\mathcal{R}_{\rm SF}$ and the SN event rate $\mathcal{R}_{\rm SN}$ (see also Owen et al., 2018) . In determining $\dot{M}$, we set $M_{\rm inj}=2\;\!{\rm M}_{\odot}$. We also use $\overline{M}_{\rm SN}=10\;\!{\rm M}_{\odot}$ as a characteristic progenitor mass for core-collapse SNe (not strictly the mean, although we find our results are relatively insensitive to the exact choice of value here).

The energy budget is specified by the energy injection rate by CRs from SN events. This relates to our system of equations by the conservation law arising from combining and integrating equation 3 and 4 (see e.g. Breitschwerdt et al. 1991; Samui et al. 2010 for details). At large radii, it is clear that the constant of integration is $Q=qv_{\infty}^{2}/2$, the kinetic energy flow of the wind. This is the rate at which energy is taken out of the host system by the outflow. The rate of energy injected per unit volume may be expressed as the sum of that injected thermally and that injected via CRs. The thermal injection rate is given by:

[TABLE]

where the fraction of available energy which goes into driving the outflow is encoded by $\nu=0.1$ (Murray et al., 2005; Samui et al., 2010) and is also applicable to the CR component. The parameter $\mathcal{Q}$ is introduced as the thermalisation efficiency which implicitly accounts for the fraction of SN energy loss that is in radiative cooling and in transforming cool clumps into ionised gas. Observations of nearby systems, e.g. M82 (Watson et al., 1984; Chevalier & Clegg, 1985; Seaquist et al., 1985; Strickland & Heckman, 2009; Heckman & Thompson, 2017) suggest that both ${\mathcal{Q}}\sim 0.1-1$ and ${\mathcal{P}}\sim 0.1-1$. However, these values are not well constrained and conflicting values are assigned for the same system in come cases (cf. Bradamante et al., 1998; Strickland et al., 2000; Strickland & Heckman, 2009; Veilleux, 2008; Zhang et al., 2014). We thus consider a conservative benchmark model of $\mathcal{Q}=0.01$ for our calculations. The other parameter, $\xi$, is the fraction of the mechanical SN energy available in the presence of energy losses by neutrino emission. For core-collapse SNe, around 99% of the SN energy is carried away by streaming neutrinos (see Iwamoto & Kunugise, 2006; Smartt, 2009; Janka, 2012), and hence $\xi=0.01$. The energy injection rate via CRs is given as:

[TABLE]

where $\zeta$ is introduced as the fraction of SN energy passed to CR power, which is then available for transfer to the outflow wind and/or hadronic interactions. We adopt a characteristic value of $\zeta=0.1$ for this, which is slightly conservative (see Fields et al., 2001; Strong et al., 2010; Lemoine-Goumard et al., 2012; Caprioli, 2012; Morlino & Caprioli, 2012; Dermer & Powale, 2013; Wang & Fields, 2018). We note that, as CRs are initially radiated isotropically away from the source region (the starburst core of the host system), we must use the solid angle fraction ${\Omega_{\rm A}}/{4\pi}$ between the interfacing outflow regions (i.e. ‘Zone A’ in Fig. 1) and the core to properly account for the fraction of initially streaming CRs which are suitably directed to be able to drive the outflow wind (this stems from the distinct two zones of our model, between which CR transfer is taken to be negligible – see the discussion about our ‘Two-Zone’ approximation in section 4.1.3 for more details and our justification of this approach). For a single outflow cone, this is given by

[TABLE]

Combining equations 10 and 11 gives the total volumetric injection power by SNe as:

[TABLE]

in terms of SN event rate $\mathcal{R}_{\rm SN}$ or star-formation rate $\mathcal{R}_{\rm SF}$, where we introduce the combined SN efficiency term:

[TABLE]

In a CR-driven outflow, some amount of the injected energy from SNe is lost in driving the flow, leaving a fraction $f$ transferred into the wind kinetic energy. For the purposes of the volumetric energy injection term into the wind fluid, we may thus use

[TABLE]

where $f$ typically takes values of a few percent, with the rest of the ‘driving’ energy being lost as the wind climbs out of the gravitational potential of the host galaxy – see Samui et al. (2010) for an analytical expression for $f$ in an NFW profile which is adopted in the outflow model used here.

2.3.1 Velocity & Density Profile

We may manipulate equations (1) to (4) as follows. First, from equations (6) and (3), and specifying that for the flow model $h>h_{\rm inj}$, we may write

[TABLE]

Following Samui et al. (2010), we now multiply equation (2) by $v$ and subtract from equation (16) above:

[TABLE]

The Alfvén velocity is given by $v_{\rm A}=B(h)/\sqrt{4\pi\rho}$ and $B(h)\;\!h^{2}$ is conserved throughout the most part of the outflow cone (scaled to give $B_{0}$ as the interstellar magnetic field at the base of the outflow). Thus, we may differentiate $v_{\rm A}$ to give:

[TABLE]

In the limit where $h_{\rm B}<h$, this simplifies to

[TABLE]

which holds for all $h$ where the second term dominates, and so is a suitable approximation for our purposes. Furthermore, by differentiating equation (1),

[TABLE]

Combining equations (4) and (6) with these results (equations 19 and 20) gives

[TABLE]

while the gas pressure $P$ can be determined from this and equation (17) as:

[TABLE]

This, together with equation (21) can be substituted back into equation (2) to give

[TABLE]

where $c_{*}$ is introduced as an effective sound speed, defined by

[TABLE]

Finally, using equation (20) to substitute the density gradient in (23) allows the velocity gradient to be written as:

[TABLE]

The location at which the outflow velocity becomes supersonic is referred to as the critical point, $h_{*}$. At the critical point, the flow velocity will be equal to the effective sound speed, i.e. $v=c_{*}$, thus the denominator of equation (25) will vanish. For a smooth velocity through the critical point (as would be expected physically), we require the numerator to vanish, and it must go to zero more quickly at this point than the denominator to ensure a regular function through this point, i.e.

[TABLE]

This allows for a useful, alternative estimate for the value of $c_{*}$ (and hence $v$) to be made at the critical point: the gravitational potential gradient may be expressed in terms of the circular velocity of the system at the critical point $V_{\rm c,*}=\sqrt{{\rm G}M(h_{*})/h_{*}}$, i.e.

[TABLE]

where $M(h_{*})$ is the enclosed mass of the system up to the critical point. As galaxy rotation curves are approximately flat at large radii, the circular velocity at the virial radius would be comparable to that at the critical point. This allows us to use the full mass of the system in place of $M(h_{*})$ to enable easier parametrisation, and means that $V_{\rm c,vir}\approx V_{\rm c,*}$. We introduce $F$ in equation 27 to account for the small difference between $V_{\rm c,vir}$ and $V_{\rm c,*}$, with $F/\sqrt{2}$ typically being of order 1 for all plausible model parameter choices (see also Samui et al., 2010). For our reference model with mass $10^{9}~{}\text{M}_{\odot}$ and SN rate of $\mathcal{R}_{\rm SN}=0.1~{}\text{yr}^{-1}$, we find a mass outflow rate of $q=0.01~{}\text{M}_{\odot}~{}\text{yr}^{-1}$, energy flux of $\dot{\epsilon}=3.0\times 10^{38}~{}\text{erg}~{}\text{s}^{-1}$ per outflow cone, a critical point location at $h_{*}=0.32~{}\text{kpc}$ and a corrective factor of $F=1.05$.

The flow velocity at the critical point can then be used as a boundary condition from which equation 25 can be integrated. We adopt a numerical approach to do this, using a 4th order Runge-Kutta method (Press et al., 2007) to integrate both inwards and outwards from the critical point. To ensure a smooth solution over the critical point, we enforce a linear gradient across it locally using the method specified in Ipavich (1975). Fig 2 shows the resulting solution for our reference protogalaxy model (solid black line). This shows that the wind tends towards a terminal velocity of around $v_{\infty}\approx 290~{}\text{km}~{}\text{s}^{-1}$, and demonstrates the comparability between the flow velocity and circular velocity around the critical point.

An associated density profile can also be found numerically from equation 20. This is an important component of the model because, in section 3.1.1, we will show that the local density of a medium determines the level of CR heating that can arise via hadronic interactions. Fig. 3 shows the resulting density profile of the outflow when adopting the same reference model parameters used for the velocity profile. This corresponds to an ISM density (within the protogalaxy) of around $10~{}\text{cm}^{-3}$, and a temperature of around 105 K at the critical point. CR and gas pressure profiles can be similarly calculated, but are not important for the analysis in this paper.

3 Interactions of Energetic Cosmic Ray Particles in Protogalactic Outflows

3.1 Cosmic Ray Interactions

3.1.1 Hadronic Processes

Proton-proton (${\rm pp}$) interactions of CRs are expected to dominate over photo-hadronic interactions at ${\rm GeV}$ energies and above in most galactic and protogalactic systems (e.g. Mannheim & Schlickeiser, 1994; Owen et al., 2018). These ${\rm pp}$-interactions produce a shower of secondary particles which include hadrons, charged and neutral pions, leptons and neutrinos (see Pollack & Fazio, 1963; Gould & Burbidge, 1965; Stecker et al., 1968; Almeida et al., 1968; Skorodko et al., 2008; Dermer & Menon, 2009). The energies carried by the energetic protons (CRs) will be distributed among their descendant particles and through their subsequent interactions and decays. In particular, through the transfer of energy to the secondary charged pion particles and then to leptons (mainly electrons and positrons), the primary CR proton can deposit a fraction of its energy into the ambient medium.

The major channels of the ${\rm pp}$-interaction are

[TABLE]

where $\xi_{0}$ and $\xi_{\pm}$ are the multiplicities of the neutral and charged pions respectively while the $\Delta^{+}$ and $\Delta^{++}$ baryons are the resonances (see Almeida et al., 1968; Skorodko et al., 2008). The hadronic products will continue their interaction processes until their energies fall below the interaction threshold $E_{\rm p}^{\rm th}=0.28~{}{\rm GeV}$ (Kafexhiu et al., 2014)444This threshold is determined from the energy required for the production of a pair of neutral pions, being the lowest energy particle produced in the cascade, where $E_{\rm p}^{\rm th}=0.28~{}{\rm GeV}=2m_{\pi^{0}}+{m_{\pi^{0}}}^{2}/2m_{\rm p}$, for $m_{\pi^{0}}$ as the neutral pion rest mass and $m_{\rm p}$ as the proton rest mass., although it is expected that this would arise after just a few interaction events (see Owen et al., 2018).

The pionic products undergo decays, where most of the neutral pions ${\pi^{0}}$ will decay into two photons through an electromagnetic process,

[TABLE]

with a branching ratio of 98.8% (Patrignani et al., 2016) and on a timescale of $8.5\times 10^{-17}\;\!{\rm s}$. Charged pions $\pi^{\pm}$ will produce leptons and neutrinos via a weak interaction,

[TABLE]

with a branching ratio of 99.9% (Patrignani et al., 2016) and on timescales of roughly $2.6\times 10^{-8}\;\!{\rm s}$.

The $\gamma$-ray photons from the $\pi^{0}$ decay have minimal interactions with the surrounding medium and they effectively stream away from their production site. The neutrinos from the $\pi^{\pm}$ decay also interact minimally with the medium. While the $\gamma$-ray photons and neutrinos allow for a net energy escape from the interaction region, the leptons, which interact strongly with the ionised and magnetised ISM and/or outflow wind medium, play a key role in mediating the energy transfer process. Although some of the energy carried by the charged leptons is lost through inverse Compton and synchrotron processes, a non-negligible fraction can still be passed to the ISM in lepton-hadron coulomb scattering and collisions.

The energy transferred to the pions can be estimated from their production cross-sections. The parameterisation of the pion-production cross-sections proposed by Blattnig et al. (2000) gives a reasonable fit to the data (with only a minor discrepancy below 50-GeV, see Owen et al., 2018), when accounting for all pion-production branches. The ratios of the primary energy then passed to the different secondary species $\{\pi^{+},\pi^{-},\pi^{0}\}$ follows as $\{0.6,\;\!0.1,\;\!0.3\}$ at 1-GeV while this tends towards around $\{0.3,\;\!0.4,\;\!0.3\}$ at higher energies. Thus, the total fraction of CR primary energy passed to charged pion production is around 0.7, of which around 0.1 is lost to neutrinos (Dermer & Menon, 2009). On their decay to secondary electrons (and positrons - hereafter, we refer to all of the charged lepton secondaries as electrons for simplicity, without losing generality) and neutrinos, around 75% of the pion energy is passed to the neutrinos, while the electrons adopt around 25% (see, e.g. Aharonian & Atoyan, 2000; Loeb & Waxman, 2006; Dermer & Menon, 2009; Lacki & Thompson, 2012; Lacki & Beck, 2013). Overall gives the fraction of the CR primary energy ultimately passed to secondary electrons as around 0.15. This is equally split among each of the electrons produced. With a multiplicity of around 4 at GeV energies (see Albini et al. 1976; Fiete Grosse-Oetringhaus & Reygers 2010 for a fitted parametrisation of multiplicity data) which dominate the CR spectrum, this gives a typical secondary electron energy $E_{\rm e}$ of around 3.75% of that of the CR primary proton energy, $E$. Thus, a GeV CR proton undergoing a ${\rm pp}$ interaction would be expected to inject approximately 4 secondary electrons, each with an energy of around 40 MeV. We introduce the parameters $f_{\rm e}=0.0375$ as the fraction of primary energy passed to electrons, and $\xi_{\rm e}$ as their multiplicity. Both are technically functions of the primary proton energy, $E$ but we set them to be constant at their dominating $E=1~{}\text{GeV}$ value for our calculations here with no discernible impact on our results.

3.1.2 Leptonic Processes & Thermalisation

In protogalactic environments, the three main processes by which the secondary electrons release their energies are radiative cooling (via inverse Compton scattering with cosmological microwave background, CMB, photons and starlight photons, and/or via synchrotron emission when interacting with the ambient magnetic field), free-free cooling (mainly due to the electron-proton bremsstrahlung processes), and Coulomb collisions in the ISM. In the high-redshift galactic environments considered here, radiative losses are mainly caused by inverse-Compton scattering with the CMB and possibly starlight, if the host galaxy is able to sustain high star-formation rates (i.e. sufficient to yield a SN event rate of above $0.1~{}\text{yr}^{-1}$ – see Owen et al. 2018). Such loses arise at a rate of

[TABLE]

per particle (e.g. Rybicki & Lightman, 1986; Blumenthal, 1970), where ${\rm c}$ is the speed of light, $\sigma_{\rm T}$ is the Thomson scattering cross-section, $E_{\rm e}$ is the electron energy and $U_{\rm i}$ is the energy density in the radiation field (or magnetic field in the case of Synchrotron losses). The rate of free-free (bremsstrahlung) cooling per particle is

[TABLE]

where $\alpha_{\rm f}$ is the fine structure constant, and the energy loss of the electrons due to Coulomb interactions in the ionised ISM is

[TABLE]

where $\ln\Lambda\simeq 30$ is the Coulomb logarithm (see Dermer & Menon, 2009; Schleicher & Beck, 2013). Further losses arises due to adiabatic expansion of the CR fluid as it propagates along an outflow. This is quantified by

[TABLE]

(Longair, 2011), and applies equally to protons and electrons. Hence, the fraction of energy carried by the CR electrons that could be deposited into the ISM is simply

[TABLE]

where $\tau_{\rm C}$, $\tau_{\rm rad}$, $\tau_{\rm ff}$ and $\tau_{\rm ad}$ are timescales of Coulomb, radiative, free-free (bremsstrahlung) and adiabatic losses respectively. Overall, we can account for the branching ratios and cooling processes by introducing the term $f_{\rm therm}(E)=f_{\rm C}|_{E_{\rm e}}\left(\xi_{\rm e}\;\!f_{\rm e}\right)|_{E}$ into our calculations, which allows us to properly estimate the fraction of CR primary energy deposited that is ultimately thermalised.

This thermalisation of the CR electron secondaries does not occur immediately. Owen et al. (2018) shows the time-scale over which this arises can be estimated as

[TABLE]

which is shorter then the dynamical timescales estimated for the system in both the advective and diffusive CR transport regimes by at least an order of magnitude (as shown in Fig. 5) and so, for our purposes, we may assume that the thermalisation of the CR secondary electrons occurs rapidly and in the vicinity of the initial ${\rm pp}$ interaction. In line with this approximation, the rate at which energy is deposited into the ambient medium per unit volume is

[TABLE]

where $n_{\rm p}$ is the local number density of protons in the wind fluid, $E$ and $n(E)=E\;{{\rm d}N(E)}/{{\rm d}E\;{\rm d}V}$ are the energy and differential number density of CR protons respectively, and ${l^{*}_{\rm p\pi}}(=1/\hat{\sigma}_{\rm p\pi}n_{\rm p})$ is the mean-free-path of the interaction. The total inelastic cross-section of the ${\rm pp}$ interaction can be parametrised as

[TABLE]

(Kafexhiu et al., 2014), where $\chi=E/E^{\rm{th}}_{\rm p}$ and $E^{\rm{th}}_{\rm p}$ is the threshold energy, as introduced above. It follows that the rate of CR attenuation is

[TABLE]

and the corresponding heating rate of the medium is

[TABLE]

The energy limits in the integral above will be discussed in §3.2.

3.2 Cosmic Ray Energy Spectrum

3.2.1 Transportation & Spectral Evolution

The transport of CRs in a bulk flow is governed by

[TABLE]

(e.g. Schlickeiser, 2002) where $n=n(E,{\bf{x}})$ is the differential number density of CR protons (i.e. the number density of CR particles per unit energy) with an energy $E$ at a location ${\mathbf{x}}$. The $\nabla\cdot\left[D(E)\nabla n\right]$ term describes the diffusion process, specified by the diffusion coefficient $D(E)$ (see §3.3), while the $\nabla\cdot\left[{\bf{v}}n\right]$ term describes the advection of CRs in a bulk flow of velocity ${\mathbf{v}}$ (see §2.3). The mechanical and radiative cooling of the CR particles is specified by the energy loss rate, $b(E,{\bf{x}})$, the injection of CRs by the source term, $Q(E,{\mathbf{x}})$, and the attenuation of CRs by the sink term, $S(E,{\mathbf{x}})$. The radiative loss timescale of protons is generally longer than that of advection and diffusion in typical galactic environments. In some galactic-scale outflows, adiabatic cooling could be important, so we retain this in our calculations and simply set $b(E,{\mathbf{x}})=\dot{E}_{\rm ad}$ in this first study. However a thorough investigation of the adiabatic cooling effect of CRs and the competition between advection and diffusion in different astrophysical settings deserves a separate investigation.

We have adopted a scenario whereby CRs are produced in SN events, and these are expected to be most frequent in active star-forming regions at the base of outflows. As such, we consider a situation where CR protons are injected at the base of the outflow wind cone as an initial boundary condition for the solution of the CR transport equation. Since our current knowledge of the properties of CRs in protogalactic environments is very limited, the initial energy spectrum of the CRs is uncertain. Given that the acceleration of CRs is a universal process, e.g. via Fermi (1949) processes in SN remnants, we boldly assume a differential energy spectrum of the freshly injected CRs similar to that observed in the Milky Way, i.e. following a power-law

[TABLE]

Here, $\Omega$ is introduced as the solid opening angle (of the outflow cone) and we adopt a power-law index of $\Gamma=2.1$, in line with Milky Way observations of the galactic ridge – a region where abundant CR injection is likely to occur, and therefore is a reflection of the ‘fresh’ CR spectrum as required here (see, e.g. Allard et al. 2007; Kotera et al. 2010; Kotera & Olinto 2011, although slightly steeper indices of around 2.3–2.4 have been suggested in recent years for pure proton data in ‘fresh’ acceleration regions, e.g. Adrián-Martínez et al. 2016; H.E.S.S. Collaboration 2018). $\mathcal{N}(...)$ is the normalisation, given by

[TABLE]

with $E_{\rm max}$ as the maximum energy of interest, and $E_{0}$ (the lowest energy under consideration) as the reference energy. We adopt a minimum energy bound of $E_{0}=1\;\!{\rm GeV}$ which corresponds to the approximate energy above which hadronic interactions may arise (via the ${\rm pp}$ mechanism – see Kafexhiu et al. 2014), and a maximum energy of $E_{\rm max}=10^{6}\;\!{\rm GeV}$ (i.e. 1 PeV) being the maximum realistic energy that could be reached by CRs accelerated in SN remnants (Bell, 1978; Kotera & Olinto, 2011; Schure & Bell, 2013; Bell et al., 2013), while higher-energy particles would likely originate from outside the protogalaxy (Hillas, 1984; Becker, 2008; Kotera & Olinto, 2011; Blasi, 2014). With a power-law index of $\Gamma=2.1$, the $1-10^{6}~{}{\rm GeV}$ range harbours more than 99% of the total energy content of the CRs (see Benhabiles-Mezhoud et al., 2013). $\epsilon_{\rm CR}$ is the CR energy density – see section 3.2.2 for details, and $v_{\rm CR}$ is the characteristic velocity with which the CRs propagate macroscopically. In the case of free-streaming CRs, this is the speed of light, ${\rm c}$. In a diffusion-dominated system, $v_{\rm CR}$ would be the diffusive speed $D(E)/\ell_{\rm diff}$ while, in an advection-dominated scenario in which CRs are trapped in a local magnetic field, but transferred along by the flow of the fluid in which they are entrained, this would be the bulk flow velocity555The microscopic CR propagation speed would remain as ${\rm c}$ in all cases, however in the diffusion and advection scenarios their macroscopic propagation appears to be much less due to the small-scale deflections and scatterings with the local magnetic field, such that their propagation can no longer be approximated as streaming.. The differential CR flux can be used to write the CR differential number density as

[TABLE]

in line with the earlier definition in section 3.1.1.

3.2.2 Cosmic Ray Energy Densities

The CR energy density depends on whether the system is dominated by advection or diffusion, and is governed by the outflow velocity (in comparison to the diffusive speed). In an advection-dominated system, the CR energy density may be expressed as

[TABLE]

where we may approximate $v(h)$ with $v_{\infty}$, the terminal velocity of the outflow, for the purpose of modelling its large-scale redistribution of CR energy. In a diffusion-dominated system, it is given by

[TABLE]

where $\ell_{\rm adv,diff}$ are the characteristic length-scales of the system when dominated by advection or diffusion.

Here, the power of the CRs, $L_{\rm CR,eff}$, is related to the power of the SN explosions injecting them into the system via

[TABLE]

where factor $f$ was first introduced in equation 15 and accounts for the fraction of energy lost by the CRs in climbing out of the gravitational potential of the host galaxy (such that a fraction $1-f$ is retained by the CRs and so is available to undergo hadronic interactions). The other symbols retain their earlier definitions (see §2.3). When accounting for the flow solid angle, the factor $\Omega_{\rm A}/4\pi$ which appeared in equation 11, is also required.

In a galaxy harbouring CRs with limited bulk flows and advection, particles diffuse throughout the volume of the host on kpc scales (see, e.g. Owen et al., 2018). As such, we adopt $\ell_{\rm diff}=1~{}\text{kpc}$ as the characteristic diffusion length-scale of the system when particle transport is well-within the diffusive regime. Conversely, if transport is dominated by advection, advective outflows extend for tens of kpc (see Veilleux et al., 2005; Bland-Hawthorn et al., 2007; Bordoloi et al., 2011; Martin et al., 2013; Rubin et al., 2014; Bordoloi et al., 2016). As such, we adopt an advection length-scale for the propagating CRs of $\ell_{\rm adv}=10~{}\text{kpc}$. In the case of an outflow system with a SN rate of ${\cal R}_{\rm SN}=0.1~{}\text{yr}^{-1}$, outflow wind velocity of $v_{\infty}\approx 290~{}\text{km}~{}\text{s}^{-1}$ (see § 2.3) and a diffusion coefficient of $3.0\times 10^{28}~{}\text{cm}^{2}~{}\text{s}^{-1}$ (appropriate for a 1-GeV CR in a 5-$\mu$G ambient mean magnetic field - see § 3.3), equation 46 and 45 would suggest associated energy densities of $\epsilon_{\rm diff}\approx 170.0~{}\text{eV}~{}\text{cm}^{-3}$ and $\epsilon_{\rm adv}\approx 0.59~{}\text{eV}~{}\text{cm}^{-3}$ in the diffusive and advective regimes respectively, i.e. the advection of CRs reduces their energy density by almost 2 orders of magnitude at the base of the outflow cone.

These energy densities are largely consistent with the CR energy densities estimated from, e.g. M82, $\epsilon_{\rm CR}\approx 550~{}\text{eV}~{}\text{cm}^{-3}$, and NGC 253, $\epsilon_{\rm CR}\approx 260-350~{}\text{eV}~{}\text{cm}^{-3}$ (Yoast-Hull et al., 2016), starburst galaxies with similar SN rates to that adopted in the present model of $\mathcal{R}_{\rm SN}\approx 0.1~{}\text{yr}^{-1}$ (Lenc & Tingay, 2006; Fenech et al., 2010). While M82 hosts a clear outflow, its CR energy density would suggest that CR propagation is still diffusion-dominated in the system overall. NGC 253 also appears to be predominantly diffusive throughout much of the galaxy, with only advective transport dominating above a height of around 2 kpc (Heesen et al., 2007), being consistent with the relatively high measured CR energy density.

Advection-dominated outflow systems would have considerably more rapid outflow velocities compared to their diffusion-dominated counterparts (to ensure that the advecting flow has a faster velocity than the diffusing CRs). A clear example of a starburst with a rapid outflow is NGC 3079, and CR propagation in this system would therefore be expected to be predominantly advective. This galaxy is known to harbour a remarkably fast outflow wind, of central velocity of around $1,100~{}\text{km}~{}\text{s}^{-1}$ but perhaps rising to nearly $3,000~{}\text{km}~{}\text{s}^{-1}$ in some regions (Filippenko & Sargent, 1992; Veilleux et al., 1994b, a; Veilleux et al., 1999)666While NGC 3079 also hosts an active nucleus (AGN), analysis by Cecil et al. (2001) has shown that the outflow wind is driven by the nuclear starburst rather than by the AGN, and so is a valid comparison here.. Radio observations of NGC 3079 indicate average CR energy densities of around $8.0~{}\text{eV}~{}\text{cm}^{-3}$, with only a small variation throughout the host (Irwin & Saikia, 2003). Given that the SN rate $\mathcal{R}_{\rm SN}\approx 0.3-0.5~{}\text{yr}^{-1}$ (Irwin & Seaquist, 1988; Condon, 1992; Irwin & Saikia, 2003), if the system were fully diffusive, we would expect the CR energy density to be around $1,000-2,500~{}\text{eV}~{}\text{cm}^{-3}$, i.e. 3-5 times that of M82 or NGC 253. Instead, this estimated value is around 100 times less than the diffusive limit prediction and is therefore consistent with the energy density predicted by the advection-limit estimation.

3.2.3 Cosmic Ray and $\gamma$-Ray Spectral Comparisons

We may compare our CR injection spectral model defined in equations 42 and 43 with $\gamma$-ray observations of the Galactic Ridge (GR) – a region of abundant gas clouds and star-formation which is likely to be a useful tracer of CR interactions and their underlying spectrum from the resulting secondary $\pi^{0}$ decays. We model the expected CR spectral energy density in protogalaxies of characteristic SN rate $\mathcal{R}_{\rm SN}=0.01,0.1,1.0$ and $10.0~{}\text{yr}^{-1}$, along with our model prediction for the Milky Way with $\mathcal{R}_{\rm SN}\approx 0.015~{}\text{yr}^{-1}$ (e.g. Dragicevich et al., 1999; Diehl et al., 2006; Hakobyan et al., 2011; Adams et al., 2013) using equation 42 and the diffusive CR energy density given by equation 46. We acknowledge that the Milky Way case differs from the protogalaxy models in that the SN types are more likely to be dominated by those resulting from lower mass stars with longer lifetimes than in a starburst protogalaxy. Such SN events have a lower characteristic energy of around $E_{\rm SN}\approx 10^{51}~{}\text{erg}$ (instead of the $E_{\rm SN}=10^{53}~{}\text{erg}$ appropriate for core-collapse SNe with massive progenitors) with less energy loss to neutrinos – $\xi$ is taken to be 0.9 for the Milky Way model (see, e.g. models and simulations in Wright et al., 2017, which suggest neutrino losses are of around a few percent of the total Type 1a SN energy), rather than the $\xi=0.01$ value appropriate for Type II core-collapse SNe (e.g. Iwamoto & Kunugise, 2006; Smartt, 2009; Janka, 2012). The size of the system is also different, with $\ell_{\rm diff}\approx 30~{}\text{kpc}$ for the Milky Way (see, e.g. Xu et al., 2015), compared to $\ell_{\rm diff}=1~{}\text{kpc}$ adopted in our protogalaxy models.

$\gamma$ rays are produced by the decay of the $\pi^{0}$ secondaries produced in the CR pp interactions according to process 29. Since all of the $\pi^{0}$ energy is passed to the $\gamma$ rays, the relation between the CR spectral energy density and that of $\gamma$-rays is governed entirely by the inelastic cross-section for the production of $\pi^{0}$ secondaries – see section 3.1.1 for details. The CR energy flux is related to the $\gamma$-ray energy flux by

[TABLE]

where the cross-sections only show a weak energy-dependence (meaning that their values at $E_{0}$ are sufficient for our estimates)777This assumes local CR isotropy, and that the vast majority of CRs are attenuated by pp-interactions in the GR region. These assumptions should be assessed more carefully in future studies, and mean that the resulting estimates for CR number density from $\gamma$-ray emissions stated here are conservative and should be regarded as a lower limit.. Equation 48 combined with equation 44 can be rearranged to allow the CR spectral energy density in a $\gamma$-ray emitting region to be estimated. Applying this to $\gamma$-ray measurements of the GR in the region between l<|$$| in Galactic longitude, and b<|$$| in Galactic latitude above 1-GeV allows the local injected CR energy density driving this $\gamma$-ray emission to be estimated, as shown in Fig. 4. This indicates the injected spectral energy density for the four protogalaxy models considered in this study, with $\mathcal{R}_{\rm SN}=0.01,0.1,1.0$ and $10.0~{}\text{yr}^{-1}$ (the four black lines, solid, dashed-dotted, dashed and dotted respectively), and also that for the Milky Way GR region, with the model scaled for the Galactic parameters discussed above for reference. This GR line is compared to values derived from $\gamma$-ray data from Fermi-LAT (the black points, see Gaggero et al. 2017)888We directly use the results from the Fermi analysis undertaken in Gaggero et al. (2017) here. These points used the Fermi Science tools V10R0P4 with 422 weeks of PASS 8 data, and event class CLEAN. See Gaggero et al. (2017) for further details of the $\gamma$-ray data analysis. and H.E.S.S. (the higher energy grey points, from Aharonian et al. 2006, also shown in Gaggero et al. 2017), which are seen to be largely consistent with our scaled Galactic model. We note, however, that the uncertainties in our model parameters are likely to be much greater than the error bars indicated in the data points, so caution should be taken in drawing strong conclusions from this comparison.

From these data points, it is evident that there may be some motivation for a slightly steeper spectral index than that used in our protogalaxy models. However, it is not clear whether this results from differences between a true protogalaxy environment and the conditions in the GR region (which may not be truly comparable to this level), or whether this could be due to systematics in the data, or shortcomings of the crude conversion between $\gamma$-ray flux and CR spectral energy density which we invoke here. The analysis in Gaggero et al. (2017) suggests a best-fit power-law index of around $-2.29\pm 0.27$ is appropriate for the H.E.S.S. data, while the combined Fermi-LAT and H.E.S.S. analysis is consistent with a slightly steeper spectral index of $-2.49^{+0.09}_{-0.08}$ (with reduced $\chi^{2}$ of 3). These suggest that our adopted index of -2.1 is reasonable enough for our purposes, and we do not believe there is sufficient tension to adopt a steeper power-law that may not necessarily be any more or less physically motivated in a high-redshift protogalaxy – particularly as the choice index does not strongly impact our results for any sensible range of values.

3.3 Diffusion Coefficient

In a uniform magnetic field, the propagation of a charged particle describes a curve characterised by the Larmor radius $r_{\rm L}$, which is given by

[TABLE]

where $q$ is the magnitude of the charge of the particle. Propagation of CRs in a medium permeated by a turbulent, tangled magnetic field is more complicated. However, $r_{\rm L}$ can be used to derive a phenomenological prescription for the CR diffusion process. The diffusive speed of the particles is expressed in terms of the diffusion coefficient $D(E)$, which accounts for their scattering in the magnetic field and turbulence. This can be quantified in the direction parallel or perpendicular to the magnetic field lines, with that perpendicular to the lines typically being around two orders of magnitude smaller in an ISM environment (see, e.g. Shalchi et al., 2004, 2006; Hussein & Shalchi, 2014, among others). Radio observations suggest a principal magnetic field component is present in outflow cones directed perpendicularly to the disk of the host galaxy – see, for instance M83 (Sukumar & Allen, 1990), NGC 4565 (Sukumar & Allen, 1991), NGC 4569 (Chyży et al., 2006), NGC 5775 (Soida et al., 2011) and NGC 4631 (Hummel et al., 1988; Hummel et al., 1991; Brandenburg et al., 1993; Mora & Krause, 2013) among others.

Thus the diffusion along the magnetic field lines directed along the outflow cone dominates the macroscopic propagation of CR particles, and as a rough approximation we may parametrise the diffusion coefficient as a random walk process with mean-free-path characterised by the local Larmor radius, i.e. as

[TABLE]

where $\langle|B|\rangle|_{h}=|B(h)|$ is the characteristic magnetic field strength in the outflow at some position $h$, and the normalisation $D_{0}=3.0\times 10^{28}$ cm2 s*-1* is comparable to observations in the Milky Way ISM (Berezinskii et al., 1990; Aharonian et al., 2012; Gaggero, 2012) for a 1-GeV CR proton in a 5-$\mu$G magnetic field (corresponding to a reference Larmor radius $r_{{\rm L},0}$). The exponent $\delta$ encodes the effects of the interstellar turbulence, for which we adopt a value of $\delta=1/2$ (see also Berezinskii et al., 1990; Strong et al., 2007), i.e. corresponding to a Kraichnan turbulence spectrum which is considered a suitable model for the ISM (Yan & Lazarian, 2004; Strong et al., 2007) – we argue it is reasonable to expect the processes driving turbulence in high redshift protogalaxies are not unlike those in the Milky Way.

In diffusion-dominated systems, observations have not shown any strong evidence for large variations of the diffusion coefficient in galactic outflows in nearby galaxies, e.g. in NGC 7462 (Heesen et al., 2016). Despite the varying magnetic field in the presence of an outflow, diffusive propagation of CRs is not likely to extend far into the outflow cone. Thus, over the relevant length-scales, we argue that the expression for the coefficient above is effectively preserved along the flow such that $D(E,h)=D(E)$, with its temporal evolution and the spatial variation determined only by the temporal evolution and the spatial variation of the local characteristic magnetic field. In a system dominated by advection, magnetic field variations would presumably yield a more significant variation of the diffusion coefficient along the outflow cone. However, in such systems, diffusion is not important over large distances with advective flows and streaming instabilities taking precedence – so whether such variation of the diffusion coefficient is present is inconsequential to our analysis.

3.4 Cosmic Ray Transport

The diffusion timescale is given by

[TABLE]

and the advection timescale may be approximated as

[TABLE]

where $\ell\approx 5~{}\text{kpc}$ is the characteristic CR propagation length-scale (note that this is not necessarily the same as $\ell_{\rm adv}$ and $\ell_{\rm diff}$ used previously, which were specific to the nature of the system under consideration – here we instead chose a consistent length-scale over which processes can be compared, and which roughly corresponds to the distance over the bulk of CRs would be found in our outflow model – see sections 3.4.1 and 3.4.2). The timescale over which the CR particles deposit their energy is determined by the attenuation of the particles due to the pp-interaction, i.e it may be expressed in terms of pp-interaction mean-free-path, as

[TABLE]

Fig. 5 shows the comparison of the three timescales for CR protons with various energies, with a 5-$\mu{\rm G}$ mean magnetic field, a mean ISM number density $n_{\rm p}=10$ cm*-3* and a flow velocity of $290~{}\text{km}~{}\text{s}^{-1}$ (being that of the terminal flow velocity established for the CR-driven outflow in section 2.3.1). We note that this is intended to illustrate the relative importance of the processes at work in this system, with the adopted conditions comparable to those at the base of the outflow (i.e. the host galaxy ISM) where densities are highest and much of the CR attenuation would arise. The true outflow model has substantially different densities, with the profile falling to several orders of magnitude lower by 5 kpc (see Fig. 3), meaning that the CR attenuation timescale would be much greater at larger distances along the outflow cone. Nevertheless, the strong attenuation near the base of the outflow will dominate the timescales, meaning the estimate here remains a suitable approximation to illustrate the global picture of the system.

For all energies, $\tau_{\rm p\pi}<\tau_{\rm adv}$, which would imply that CR protons are substantially attenuated near the base of the outflow, where conditions are most similar to those assumed in the Fig. 5 approximation. However, over a length-scale comparable to the size of the host galaxy (of a radius of $\sim 1$ kpc), the time over which advection would arise is an order of magnitude less than that shown in Fig. 5. This means that absorption and advection would operate over comparable timescales, and a non-negligible fraction of CRs could be advected by the flow to reach distances beyond the galaxy from which they originate. Indeed, taking full account of the true density profile of the outflow means that attenuation would be substantially reduced compared to the situation indicated here – meaning that a substantial fraction of the CR energy density could be deposited outside of the host galaxy instead of within it. CR protons with $E_{\rm p}<10^{3}\;\!{\rm GeV}$ have a long diffusion time. With $\tau_{\rm diff}\gg\tau_{\rm p\pi}$, in the absence of advection, these CR protons will be contained within the galaxy and eventually release all their energy through the pp-interaction. CR protons with $E_{\rm p}\gg 10^{3}\;\!{\rm GeV}$, which have $\tau_{\rm diff}<\tau_{\rm p\pi}$ and $\tau_{\rm diff}<\tau_{\rm adv}$ (for outflows with $v\approx 290\;\!{\rm km~{}s^{-1}}$ or slower), could, however, diffuse out of the galaxy regardless of whether advection is present or not. However, only a very small fraction of the total CR energy density are harboured in this part of the CR spectrum, so their effects would not be of great astrophysical importance.

The timescale comparison gives a qualitative assessment of the relative importance of the advective and diffusive processes in the context of CR heating. A more quantitative analysis requires us to solve the transport equation (equation 41) explicitly, which we discuss in the remainder of this section. In our solution scheme, we consider that the system has settled into a steady state, which implies that we may set $\partial n/\partial t=0$. We adopt a numerical scheme in which CR protons are only injected at the base of the outflow cone (i.e. from the actively star-forming region), which practically transforms the source term $Q(E,{\rm x})$ into a boundary condition. However, observationally, the SNe sources of CRs can be distributed some way into an outflow. In extragalactic studies, the majority of SNe are found in galaxies up to around half of their estimated scale radius  (see, e.g. Hakobyan et al. 2012, 2014, 2016 which consider SNe in host galaxies up to 100 Mpc away). At higher redshifts, which would be most relevant to the starburst protogalaxies we model here, ISM conditions would presumably be more turbulent due to the higher SN activity and this may lead to a proportionally higher distribution of sources throughout the host system. Adopting a single boundary condition for the injection of CRs at the base of an outflow is thus insufficient to model the distribution of CR sources, particularly as we are intend to calculate the CR heating effect well within the ISM region of the outflow, down to 100 pc. We therefore calculate the outflow (both in the advection and diffusion cases) as the linear sum of scaled outflow solutions by a Monte-Carlo (MC) method, as outlined in section 4.1.2. In the next section we show calculations of two regimes: firstly, when the transport is dominated by advection and, secondly, when the transport is dominated by diffusion. We solve the transport equation explicitly in these two regimes, before accounting for the distribution of SN sources in the galactic core.

3.4.1 Advection Dominated Regime

In the advection dominated regime, we may drop the diffusion term. This reduces the transport equation (in the steady state) to

[TABLE]

(Here and hereafter, unless otherwise stated, we use the short-hand notation $n=n(E,{\bf x})$.) Suppose that the flow follows streamlines in the outflow cone. By symmetry, the flow is essentially 1-dimensional (specified by the co-ordinate $h$), and the transport equation, when the outflow has settled into a velocity profile $v(h)$ (see section 2.3), may be expressed as

[TABLE]

where the sink term now takes the form $S(E,h)={\rm c}~{}n~{}\hat{\sigma}_{\rm p\pi}(E)~{}n_{\rm p}(h)$. With the substitution $Z(E,h)=n(E,h)\;v(h)\;h^{2}$, the transport equation becomes

[TABLE]

The variable $b(E,h)=-{{\rm d}E}/{{\rm d}t}$ is known when the cooling processes are specified, and $\partial b(E,h)/\partial E$ can be found from this. When $b(E,h)$, the sink term and the boundary condition at the base of the outflow cone $h_{0}$ are set, the transport equation can be solved numerically using a finite difference method as described in Appendix A. In this work, we solve the equation for the case of the only non-negligible cooling process being that due to the adiabatic cooling of CRs propagating along the outflow cone, i.e. $b(E,h)=\dot{E}_{\rm ad}$, subject to the boundary condition that $Z(E_{0},h_{0})=n(E_{0},h_{0})~{}v(h_{0})~{}h_{0}^{2}$ with $h_{0}$ set to be $100~{}{\rm pc}$, the size of the starburst region (see Chevalier & Clegg, 1985; Tanner et al., 2016) and we use a reference energy $E_{0}$ at $1~{}{\rm GeV}$. Moreover, $Z(E,h)\propto E^{-\Gamma(h)}\;\!h^{2}\;\!v(h)$ with a power-law index $\Gamma(h_{0})=2.1$ at $h_{0}$. We invoke appropriate velocity profiles $v(h)$ modelled according to section 2.3.

By inspection of the transport equation, we may see that in the absence of energy-dependent CR cooling, softening (or hardening) of the CR energy spectrum over a large-scale galactic outflow will not occur in the advection dominated regime. This conclusion can also be reached in a qualitative analysis by comparing the attenuation and advection timescales and their energy dependencies. As shown in Fig 5, $\tau_{\rm adv}$ is independent of the CR energy and $\tau_{\rm p\pi}(E)$ has only very small variations across the energy range considered. Thus, without strong energy dependences in these two terms, there should not be significant evolution of the energy spectrum of the CRs that are advected by the flow. Fig. 6 shows the CR energy spectra (obtained by solving the transport equation numerically) over distances up to $h=50~{}{\rm kpc}$ for an outflow with a full opening conic angle of [math] and the velocity profile determined in section 2.3.1, which indicates negligible spectral evolution of the CRs along the flow.

3.4.2 Diffusion Dominated Regime

In the diffusion dominated regime, the transport equation takes the form

[TABLE]

The outflow cone is axi-symmetric and so the transport equation is 1-dimensional, specified by the coordinate $h$ (as in the advection dominated case). If the diffusion coefficient does not vary significantly along $h$, then we have

[TABLE]

Substituting $Z(E,h)=h^{2}n(E,h)$ into the equation yields

[TABLE]

After rearranging and expanding the energy derivative, we obtain

[TABLE]

The transport equation is solved numerically, with the scheme described in Appendix A. This requires two Neumann boundary conditions (i.e. step 1 in equation 75), and these are obtained directly from considerations of SN event rates, and the efficiency of CR production. Here we show how the two boundary conditions are constructed.

For the first one, we begin with

[TABLE]

where ${{\rm d}Z}/{{\rm d}h}$ at the boundary needs to be specified. As we lack a prescription that accounts for the acceleration of CRs in the starburst region with a transition to their transportation in the outflowing region, we adopt the assumption that ${{\rm d}Z}/{{\rm d}h}$ scales with that of CR electrons initially at the lower boundary of the wind cone. Hence, we have

[TABLE]

where $R$ is the scaling factor. Radiative cooling processes are generally inversely proportional to the fourth-order of the mass of the charged particles. Thus, we set $R=(m_{\rm e}/m_{\rm p})^{4}$, which implies that ${{\rm d}Z}/{{\rm d}h}$ is negligible for CR protons. This prescription is consistent with the CR proton flux being conserved at the boundary, i.e. $v(h_{0})Z(E)|_{h_{0}}=v(h_{0})\{n(E)h^{2}\}|_{h_{0}}={\rm constant}$. We may estimate the value of ${{\rm d}n_{\rm CRe}}/{{\rm d}h}$ (the gradient in CR electron number density) from radio observations. Observations of the nearby starburst outflows in NGC 7090 and NGC 7462 in the 6 cm and 22 cm bands (Heesen et al., 2016) indicate

[TABLE]

at the base of the wind. The second condition relates to the rate of change of the CR spectral index at the base of the outflow cone. Similarly, we assume a scaling with the CR electrons. Observations of the nearby starburst galaxies NGC 7090 and NGC 7462 (Heesen et al., 2016) suggest that $\Gamma_{\rm CR,e}^{\prime}={\rm d}\Gamma_{\rm CR,e}/{\rm d}h=-1.4$ at the base of the galactic outflow. Thus,

[TABLE]

where $R$ is as defined above. This requires that the cooling and spectral evolution at the lower boundary is insignificant, which, in turn, ensures a negligible variation in the spectral index of the CR protons.

We adopt the same numerical solution scheme used in solving the transport equation in the advection dominated regime. The solution $Z$ is obtained by integrating the transport equation with the two boundary conditions at the base of outflow cone using a Runge-Kutta method (as described in the Appendix A). The result is shown (for $h$ up to $2~{}{\rm kpc}$) in Fig. 6999The values shown in Fig. 6 are calculated for a protogalaxy model with $\mathcal{R}_{\rm SN}=0.1~{}\text{yr}^{-1}$. If we instead scale to a Milky Way-like model as described in section 3.2.3 then the CR differential fluxes are a little lower. Assuming an inelasticity of approximately 0.3 for the production of $\gamma$-ray producing neutral pions (see section 3.1.1), the resulting $\gamma$-ray flux can be estimated. Given that the Milky Way is predominately diffusive in terms of CR propagation, the CR differential fluxes according to the above spectral model approach would be around $1.0\times 10^{-5}\;{\rm ph}\;{\rm GeV}^{-1}\;{\rm cm}^{-1}\;{\rm s}^{-1}\;{\rm sr}^{-1}$ at 1-GeV, or $7.0\times 10^{-12}\;{\rm ph}\;{\rm GeV}^{-1}\;{\rm cm}^{-1}\;{\rm s}^{-1}\;{\rm sr}^{-1}$ at 1-TeV. This is consistent with $\gamma$-ray flux measurements of the Galactic Ridge above $E_{\rm\gamma}=10~{}\text{GeV}$ with Fermi-LAT, H.E.S.S. and VERITAS in e.g. Aharonian et al. (2006); Gaggero et al. (2015a); Archer et al. (2016); Gaggero et al. (2017).. The ratio of the normalised diffusion spectrum to the normalised advection spectrum gives an indication of the level of attenuation $\mathcal{A}_{\rm D/A}$ experienced by the CRs as they diffuse, as this is the cause of the turn-over in the spectrum – this is shown in Fig. 7, which illustrates the dominant effect of attenuation of the diffusing CRs up to around $5\times 10^{3}$ GeV, at which point the diffusion process becomes more important.

4 Results and Discussion

We consider starburst protogalaxies at high redshift, which studies have indicated could host substantial star-forming activity (e.g. Hashimoto et al., 2018; Watson et al., 2015) and strong galactic outflows (e.g. Steidel et al., 2010). Such star-forming activity makes these systems a likely host of abundant CRs – but here we primarily address whether they would predominantly operate as CR calorimeters (Thompson et al., 2007; Lacki et al., 2011), or whether strong outflow activity could provide a means by which CRs can escape and interfere with the circumgalactic and/or intergalactic environment. To help our understanding of how CR containment in young star-forming galaxies progresses and how the coexistence of diffusion and advection in different regions of the same host develops, we consider starburst protogalaxies that would be present at redshift $z\approx 7$, being among most distant objects that may be observed with current and near future deep optical/UV surveys, e.g. the Subaru HSC deep-field survey (The HSC Collaboration, 2012).

4.1 Hadronic Heating in Outflows

We have so far only considered outflow environments where the propagation of CRs is dominated by either advection or diffusion. In reality, both advection and diffusion would presumably operate simultaneously and a proper treatment of CR transport in an outflow would require a complete solution of the full transport equation. However, at any one location and particle energy, the dynamics would usually be dominated by only one of these processes. For instance, in regions where bulk velocities are low, diffusion would likely be more important than advection. This would be the case in regions of near the base of the outflow cone, where the low outflow velocity would not be able to compete against CR diffusion – indeed, such an effect is seen in numerical simulations (e.g. Farber et al., 2018). The opposite would be true at high altitudes, where the flow velocity is greater and could advect CRs faster than they would typically be able to diffuse. The relative importance of the contributions from each of these two process along an outflow would impact on the distribution of CRs and, by equation 40, would govern the location at which they deposit energy and thermalise.

4.1.1 Concurrent Advection & Diffusion

We may attain a reasonable approximation for the distribution of CRs in a system where both advection and diffusion operate by weighting the pure advection and pure diffusion limit solutions by their respective timescales at each position and energy, and summing these contributions together. Evaluating advection and diffusion timescales at each calculation increment accounts for both the variation of flow velocity over position as well as the variation of the diffusion coefficient over energy. The associated effective hadronic heating rate $\bar{H}$ in the concurrent advection/diffusion picture along the outflow then follows as:

[TABLE]

with

[TABLE]

where

[TABLE]

and

[TABLE]

Here, $\tau_{\rm adv}=h/v(h)$) is the position-dependent advection timescale, while $\tau_{\rm diff}=h^{2}/4D(E)$ is the energy-dependent diffusion timescale.

Individual advection-dominated and diffusion-dominated hadronic heating profiles are shown as the two black lines in Fig. 8, while the concurrent advection/diffusion heating power is indicated by the dashed red line. This demonstrates how diffusive propagation is important in the inner regions of the outflow, while advection dominates at higher altitudes above 0.4 kpc. Above this point, the outflow velocity is sufficiently greater compared to the typical diffusive speed of the CRs (see also Fig. 2). This is calculated when adopting a SN-event rate of $\mathcal{R}_{\rm SN}=0.1~{}\text{yr}^{-1}$, a conical galactic outflow of an opening angle of [math], an outflow thermalisation efficiency of $\mathcal{Q}=0.01$, and mass-loading factor of $\mathcal{P}=0.1$. For reference, these choices yield an outflow terminal velocity of $v_{\infty}\approx 290~{}\text{km}~{}\text{s}^{-1}$ and are the same as those used to produce the profile in Fig. 2.

4.1.2 Extended Starburst Region & Computational Scheme

Observations have indicated that SN events can arise throughout the disk of their host, at least up to around half of their estimated scale radius (Hakobyan et al., 2012, 2014, 2016). While sections 3.4.1 and 3.4.2 and the previous heating profile in Fig. 8 are suitable descriptions for CRs in outflows emerging from a single point-like region in the centre of their host galaxy, it is also necessary to consider the impact of a more physical extended core region. To do this, we adopt a Monte-Carlo (MC) scheme to simulate a spherical distribution of $N$ points up to 0.5 kpc from the protogalactic centre (being half the adopted scale-radius for the protogalaxy). We find that a choice of $N=1,000$ points yields a sufficient signal to noise ratio. The distribution of CRs calculated according to equation 56 or 3.4.2 is scaled by $1/N$. The scaled profile is then convolved with the MC spherical points distribution, with each point being taken to be a linearly independent $h_{0}$ boundary condition101010For this, we adopt a uniform spherical injection region as a first model. More detailed injection distributions are beyond the scope of the current paper but, e.g. a singular isothermal self-gravitating spherical injection weighting is proposed by Rodríguez-González et al. (2007) – or, see also Silich et al. (2011); Palouš et al. (2013) for other approaches. We found the choice of injection model, if reasonable, does not bear any strong influence on the results presented here.. The ensemble of individual CR profiles are then superposed to give a resulting total CR distribution in the outflow, and this accounts for the extended distribution of driving CR sources. This approach is applied to both the pure diffusion and advection calculations as well as in the combined case where both processes operate. This resulting summed CR distribution can then be used to determine the hadronic heating profile using the same method as per equation 65, and is shown in Fig. 9 where lines retain their earlier definitions, and where the starburst extended CR injection region is indicated in blue. This demonstrates the broadened profile for the extended injection, and shows how the transition between the advection and diffusion dominated transport zones is determined by the relative timescales over which they operate rather than the region in which the CRs are injected. A distinct picture of a lower ‘diffusion’ region in an outflow emerges in Fig. 9, with an ‘advection’ region at higher altitudes where the flow velocity is faster. This result follows the wind structure first introduced in Breitschwerdt et al. 1993 (see also Recchia et al. 2016).

4.1.3 Two-Zone Heating Rates

Our discussion up to this point has been predominantly concerned with the CR dynamics and heating distribution arising in the outflow cone, i.e. that labelled ‘Zone A’ in the schematic in Fig. 1. However, this only paints part of the picture: typically, a star-forming galaxy be unlikely to be enveloped entirely be a galactic outflow. In disk galaxies in particular, the outflow morphology would normally be bi-conical in nature (cf. section 1 and, e.g. Strickland et al. 2000; Ohyama et al. 2002; Veilleux et al. 2005; Cooper et al. 2008). Thus, as illustrated in Fig. 1, there would usually be a substantial region of the host galaxy that is not directly influenced by the outflow – and in this region (‘Zone B’), CR propagation would presumably operate predominantly by diffusion over kpc scales. It follows that CR heating in Zone B would therefore exhibit similar characteristics to the inner region of Fig. 9, where the advective flow is too slow to have any important effect on the redistribution of CRs.

We may broadly use the results calculated in Owen et al. (2018) to model the expected CR heating effect in Zone B (and the lower regions of Zone A), where purely diffusive CR transport was invoked for a protogalaxy of otherwise similar specifications to that considered in the present work. In the earlier study, the volumetric CR heating rate was calculated throughout a protogalactic ISM due to its irradiation by the entire CR emission from the galaxy. By contrast, in the current work we assume that CRs cannot readily propagate between the two zones (hereafter the ‘Two-Zone approximation’), which would mean that the CR irradiation experienced within each zone would be limited to the fraction of the total galactic CR emission which passes into that zone. We may evaluate this fraction by considering the solid angle subtended by each of the zones on the central starburst core, i.e. the CR power passing into Zone A would be

[TABLE]

(where the factor of 2 accounts for the bi-polar nature of the outflow) and that into Zone B follows as

[TABLE]

for $L_{\rm CR,eff}$ as the total CR power (see equation 47).

To properly attribute the results of Owen et al. 2018 at $\mathcal{R}_{\rm SN}=0.1~{}\text{yr}^{-1}$ to the CR heating levels expected in Zone B of the present study, a scaling by solid angle in both regions must be used. The resulting heating rate per steradian can then be used to compare heating power across the two zones, and is shown for Zone A in blue and Zone B in red in Fig. 10. In Zone B, the CR heating power is only calculated up to around 0.5 kpc, i.e. well within the ISM of the host – this is because the structure of the magnetic fields connecting the ISM to the circumgalactic medium is unclear and falls beyond the scope of this discussion. We further include lines for diffuse X-ray heating, heating from stellar radiation and heating from freely streaming CRs (which would arise if no magnetic field were present), as indicated by the relevant legend. These are intended for use as a comparison; for details regarding how these lines are calculated, we refer the reader to Owen et al. 2018. In Fig. 10, it is evident that in Zone B, the CR heating rate is dominant and even exceeds that due to the more conventional radiative processes. In Zone A, the CR heating rate is much reduced compared to the radiative heating lines, but still has the potential to play an important role in influencing the thermal properties of the ISM – particularly in the diffusion-dominated inner region, where it is maintained at a comparable level to Zone B.

We justify our use of the Two-Zone approximation as follows: we anticipate that the magnetic structure within the outflow would be perpendicular to the plane of the host galaxy and thus perpendicular to the magnetic field orientation within the Zone B ISM region. Indeed, such perpendicular magnetic structure in outflows is seen in simulation work where, e.g. the action of a CR-driven dynamo yields a perpendicular magnetic field configuration compared to the host galactic plane (Kulpa-Dybeł et al., 2011), or by the advection of the magnetic fields by the flows themselves (Bertone et al., 2005), by magnetic amplification via the CR streaming instability (Uhlig et al., 2012) along the outflow. This magnetic structure would also be consistent with polarised radio synchrotron emission above and below the planes of galaxies known to host outflows in the nearby Universe, with the polarisation direction aligned with the orientation of the outflow cone (see, e.g. Hummel et al., 1988; Sukumar & Allen, 1990, 1991; Hummel et al., 1991; Brandenburg et al., 1993; Chyży et al., 2006; Soida et al., 2011; Mora & Krause, 2013). We argue that the principle mechanism for CRs to permeate the Zone A/Zone B interface would be via diffusion. With magnetic field lines aligned in a direction parallel to the inter-zone boundary, diffusion across the interface would be severely hampered – the cross-boundary diffusion coefficient would effectively be perpendicular the the local magnetic field lines, and so would be around two orders of magnitude smaller than that along the field directions (e.g. Shalchi et al., 2004, 2006; Hussein & Shalchi, 2014), and substantially less than the effective ISM diffusion coefficient. The detailed substructure of the magnetic fields in these interfacing regions is not yet fully understood (Veilleux et al., 2005), but we argue that our prescription is consistent with existing work on relevant scales and that adopting an alternative model for CR transport across this boundary at this point would not imply an interpretation that is any more physical than that adopted here. We acknowledge that, in future studies, it will be critical to assess the magnetic fields in these interfacing regions across a range of length-scales to properly determine the permeability of the Zone A/Zone B interface to diffusing CRs.

4.1.4 Energy Deposition

For a SN event rate of $\mathcal{R}_{\rm SN}=0.1~{}\text{yr}^{-1}$, the total CR luminosity due to SN events in the model protogalaxy is $3.0\times 10^{41}~{}\text{erg}~{}\text{s}^{-1}$ (see equation 47), of which $1.1\times 10^{39}~{}\text{erg}~{}\text{s}^{-1}$ is passed into the outflow cone (when accounting for the geometry of the system and energy lost by the CRs in driving the outflow), being available to heat the ambient gases via hadronic interactions. In Fig. 11), it can be seen that the total integrated heating effect up to 50 kpc – the characteristic extent of a CR-driven galactic outflow (see, e.g. Veilleux et al., 2005; Bland-Hawthorn et al., 2007, for a discussion on the extent of outflows) – along the cone is around $2.0\times 10^{35}~{}\text{erg}~{}\text{s}^{-1}$, with much of this energy being deposited within the inner 0.3 kpc (even if the injection of CR energy is extended throughout a starburst region). We note that the thermalisation process of CR primaries is relatively inefficient, so their total heating power is much less than the energy released by the CR protons as they undergo hadronic interactions (see section 3.1.2 for a discussion on the respective loss channels compared to thermalisation). For instance, if we account for the power-law spectral energy distribution of CRs and adopt the mean energy of 8.2 GeV as a characteristic value, a CR would thermalise only a fraction of $1.5\times 10^{-3}$ of its initial energy on average111111This fraction would vary along the outflow cone due to the density profile of the wind fluid. The value quoted here, and those hereafter, is an average weighted by the heating profile to give a characteristic value for the hadronic heating efficiency along the outflow.. This means that the total power released by the CRs would be $1.4\times 10^{38}~{}\text{erg}~{}\text{s}^{-1}$ as they interact, being around 12.0% of the available CR power in the outflow cone. This fraction is consistent with the level found in the numerical models of Girichidis et al. (2018) (between 5 and 25%), which consider a combined transport scenario where outflow velocities close to the mid-plane are small, much like the diffusion limit and combined transport picture considered in the present work. Farber et al. (2018) also considers a pure advection scenario. They find substantial CR energy is harboured within a few kpc of the galactic plane because outflow velocities are low near $h=0$. The advected CRs therefore spend more time in the vicinity of the galactic plane and base of the outflow, and so are more likely to undergo hadronic interactions there. Moreover, the CR number density distributions resulting from the pure advection models of the Farber et al. (2018) study suggest a difference of around 4 orders of magnitude between the galactic disk and halo (when they invoke CR diffusion and coupling, the contrast falls to a little under 2 orders of magnitude). When further accounting for the density contrasts between the same locations (around 1 order of magnitude in the advection scenario), the CR hadronic heating rates in the plane and halo at around $h\approx 1~{}\text{kpc}$ would presumably differ by around 5 orders of magnitude – i.e. similar to the contrast in heating power found in comparable locations in this work. By contrast, Pakmor et al. (2016) find substantial CR energy (of order a few $\text{eV}~{}\text{cm}^{-3}$) up to heights of 5 kpc or more in similar systems when adopting even lower outflow velocities (around 100$~{}\text{km}~{}\text{s}^{-1}$) than those considered here – however, their calculations focus predominantly on CR transport, and the authors indicate that the hadronic interactions of the CRs had not been included in their simulations.

4.2 Model Parameters

We find that the critical quantities which govern the behaviour of the system are the SN event rate $\mathcal{R}_{\rm SN}$ and, to a lesser extent, the outflow velocity. Other parameters bear less importance – for example the influence of the outflow opening angle is demonstrated in Fig 12 where the black line is the result for the standard opening angle of [math] while the blue and red lines show a variation of +[math] and -[math] respectively (i.e. [math], blue and [math], red). While the impact here is minimal, it can be seen that, for larger opening angles, the heating rate is slightly higher: this is due to the greater fraction of the starburst core being subtended by the outflow cone and the associated higher injection rate of CR particles.

It is also of interest to consider the impact of directly modifying the outflow velocity. This is because different wind-driving mechanisms can yield vastly different outflow rates (Veilleux et al., 2005). To ensure a reasonably self-consistent estimate, the density profile must be scaled in conjunction with the velocity: from equation 1, the product $\rho\;\!v(h)\;\!h^{2}$ is conserved. So, for an equivalent mass injection rate, a velocity scaled by a factor $\Psi$ would require a corresponding scaling of the density profile by factor $\Psi^{-1}$. We note that this also results in a change in the associated mechanical energy of the wind, which is also scaled by a factor of $\Psi$ as follows from $\rho\;\!v_{\infty}^{2}$. The effect of scaling the system in this way is shown in Fig. 13, with the result essentially following directly from the initial condition for the CR energy density at the base of the outflow combined with the variation in the outflow density profile (which results from scaled velocity profile). In equation 45 it can be seen that the CR number density, and hence any heating effect the CRs may drive, is inversely proportional to $v_{\infty}$. Thus, if the velocity is reduced by an order of magnitude, the CR number density increases by an order of magnitude. However, the density of the outflow would also increase by an order of magnitude in such a scaling. Overall, this would increase the CR heating power in the advection limit by two orders of magnitude as it is proportional to both the CR number density and the density of the outflow wind fluid (which provides the target hadrons for the CR heating effect). In the diffusion limit, the picture is more straightforward: this time there is no dependence of the CR number density on the scaling, so the the CR heating power would only increase by one order of magnitude. We find that there is little bearing on the fraction of CRs absorbed in the outflow cone when a scaling is applied: when reducing the flow velocity by a factor of 10, the absorbed fraction falls slightly to 11.9% (compared to 12.0% in the baseline case – see section 4.1.4), while the fraction increases to 13.7% if increasing the outflow velocity by the same factor.

The most directly influential parameter is the SN event rate, $\mathcal{R}_{\rm SN}$, which essentially specifies the energy budget of the system. The CR heating rates scale with $\mathcal{R}_{\rm SN}^{2}$, as seen in Fig. 14, where a plausible range of values from $0.01~{}\text{yr}^{-1}$ (fairly quiescent, barely star-bursting system) to $10.0~{}\text{yr}^{-1}$ (extremely violent star-forming environment, perhaps possible in a system undergoing a major merger121212For example, some studies have suggested that Arp 220 could have a SN event rate of 2-4$~{}\text{yr}^{-1}$ (Lonsdale et al., 2006; Varenius et al., 2017). In Lonsdale et al. (2006), 4 new radio sources were observed in 12 months, and the authors argue that this is consistent with 4 new SN events and a corresponding SN rate, although this is based on limited statistics. The more recent study by Varenius et al. (2017) suggests a rate of $4\pm 2~{}\text{yr}^{-1}$ on the basis of the number of events observed, and by assuming only a small fraction of SN events are observable. So, while an event rate of 4$~{}\text{yr}^{-1}$ may be possible, the true level is likely to be lower than this (Page, 2018). We adopt the maximum value of 10$~{}\text{yr}^{-1}$ as an extreme case in the present work for completeness, but acknowledge this should be treated very much as an upper-limit.) are explored. Both the density of the outflow wind fluid (being the target for the hadronic interactions) and CR emission are dependent on $\mathcal{R}_{\rm SN}$. The CR heating power depends on the product of these two quantities, and so the dependence on $\mathcal{R}_{\rm SN}^{2}$ follows from this. The 12% fraction of CR proton energy absorbed in the outflow by hadronic interactions (as calculated in section 4.1.4) is independent of SN event rate, as is the outflow velocity profile. In both cases, this is because the contributing quantities which are dependent on $\mathcal{R}_{\rm SN}$ scale antagonistically.

Fig. 15 shows the timescales for various loss processes of secondary CR electrons near the base of the outflow at $h=0.1~{}\text{kpc}$ compared to the upper regions $h=50~{}\text{kpc}$. Since the microphysics which affect these secondary electrons is instrumental in governing the energy fraction of the CR primaries which can thermalise, it is useful to understand why the efficiency of CR hadronic heating may vary depending on the environment: the left panels (where lines are plotted in black) show that, when changing from a SN event rate of $\mathcal{R}_{\rm SN}=0.01~{}\text{yr}^{-1}$ to $\mathcal{R}_{\rm SN}=10.0~{}\text{yr}^{-1}$, the intensity of starlight should increase so as to proportionally reduce the inverse Compton timescales. In the absence of an outflow, this would lead to greater inverse Compton losses (due to scattering off a more intense stellar radiation field), a lower CR thermalisation efficiency and a resulting smaller increase in CR heating rate than a direct proportionality with $\mathcal{R}_{\rm SN}$ would suggest (see Owen et al., 2018). However, in the presence of an outflow, the density profile is enhanced by greater mass-loading rates from the SN winds, thus yielding an environment which is better able to thermalise CR secondaries. Moreover, this denser environment promotes hadronic interactions to boost CR heating. This is sufficient to roughly maintain the proportionality seen in Fig. 14, even though inverse Compton losses are more severe in more intense starlight (when $\mathcal{R}_{\rm SN}$ is greater). The underlying complexity of this behaviour also accounts for the slight deviation from a direct and consistent scaling between the lines in the advection limit at high altitudes in Fig. 14. While most of the processes in Fig. 15 operate at an appreciable rate (i.e. with timescales much shorter than the Hubble timescale) at the base of the outflow, only inverse Compton losses off the uniformly distributed CMB photons (and some starlight inverse Compton for the highest energy particles when the SN event rate is sufficiently high) are important at high altitudes. This demonstrates that the higher in altitude the CRs reach within an outflow, the less likely they are to thermalise. Indeed, above a few kpc, it can safely be assumed that the CRs have been fully advected out of their host galaxy and would instead impart their effects on the external environments (e.g. in the circumgalactic and inter-galactic medium).

4.3 Implications of Cosmic Ray Heating Power in and around Protogalaxies

The CR heating power due to hadronic interactions could be as high as $10^{-25}~{}\text{erg}~{}\text{cm}^{-3}~{}\text{s}^{-1}$ in a protogalaxy with a SN-rate $\mathcal{R}_{\rm SN}=0.1~{}\text{yr}^{-1}$, if outflows or galactic winds do not operate, such that strong CR containment within the host galaxy arises (Owen et al., 2018). We find that, if outflows are present, in the lower regions of the outflow cone (where CR diffusion still dominates over advection) a comparable heating power of $10^{-26}~{}\text{erg}~{}\text{cm}^{-3}~{}\text{s}^{-1}$ can also be attained. In the pure advection limit (i.e. when neglecting diffusion), a comparatively reduced CR heating rate is redistributed along an outflow cone compared to the diffusion limit, while a substantial fraction of the CRs do not engage with hadronic interactions in the outflow wind at all, instead escaping from the system entirely. These CRs may deposit their energy over much greater regions in the vicinity of the source galaxy and beyond. This could affect pre-heating and ionisation processes in the wider Universe (Nath & Biermann, 1993; Sazonov & Sunyaev, 2015; Leite et al., 2017), or alter the dynamics of the circumgalactic medium.

4.3.1 Cosmic Ray Containment

In the diffusion limit, CRs can become contained within the ISM of their host galaxy where they deposit a substantial fraction of their energy. This is consistent with $\gamma$-ray observations of nearby starbursts - e.g. M82 - which are $\gamma$-ray bright, suggesting that some non-negligible fraction of the CRs interact to produce pion decay $\gamma$-ray emission within the galaxy, rather than being advected away (see, e.g. Abdo, A. A. et al., 2010; Wang & Fields, 2014; Yoast-Hull et al., 2016; Heckman & Thompson, 2017; Wang & Fields, 2018). Since CR heating preferentially targets denser regions of the ISM (due to the particles of the interstellar gas being the principal targets in CR hadronic interactions), this could severely affect the ability of molecular clouds and cores to collapse into stars by raising their effective Jeans’ mass. This would push the initial mass function of forming stars to a more top-heavy form – or could even quench star-formation entirely, leading to a quiescent period until the ISM gas has sufficiently cooled to allow star-formation to resume. While the magnetic fields of star-forming regions would presumably also influence the local level of CR heating experienced, it is difficult to envisage a situation where magnetic shielding can act to such a level that there would be no impact by CR heating at all – indeed, the opposite effect may be true if magnetic field vectors are arranged in a way so as to preferentially direct diffusing CRs into denser, star-forming regions. We leave such detailed modelling to future work.

4.3.2 Cosmic Ray Escape

The picture is very different in the advection limit: in section 4.1.4, we find that around 12% of the total CR power within an outflow can be absorbed by hadronic processes, with the remaining 88% effectively being transported into the circumgalactic medium (and beyond). Although the exact fraction would vary depending on the details of the model adopted and parameter values used, this result would suggest that a substantial fraction of CR energy in an outflow is actually able to be transported away by advection if caught up by the flow. Moreover, our discussion in section 4.2 would indicate that this picture is not strongly sensitive to the choice of model parameters, with the majority of CRs in an outflow cone escaping for any reasonable parameter choices. In terms of energetics, equation 69 indicates that around 11.2% of the total CR luminosity of a host galaxy would pass into an outflow (with opening angle [math]). With the 12% hadronic absorption fraction of this arising within the cone, it follows that around 10% of the total CR luminosity of the source galaxy would escape into the circumgalactic and/or intergalactic medium – a CR power that would correspond to around $3.1\times 10^{40}~{}\text{erg}~{}\text{s}^{-1}$ at a SN event rate of $\mathcal{R}_{\rm SN}=0.1~{}\text{yr}^{-1}$. This represents a substantial contribution to the exterior energy budget around a starburst galaxy, contributing a power greater than the expected diffuse X-ray emission (at around $10^{38}~{}\text{erg}~{}\text{s}^{-1}$ – see, e.g. Watson et al. 1984) and is an intriguing result: while an outflow can facilitate the advection of substantial levels of CR energy into the extragalactic environment, this would only correspond to around 10% of the total CR energy being provided to the system by SN events. Thus, even in the presence of outflows, the CR calorimetric ability of starburst galaxies is likely to still remain relatively effective with a large fraction of CR energy still being contained and deposited within the ISM. This contained fraction would perhaps be slightly enhanced shortly after the onset of star formation before sustained outflow activity has had time to develop, or when star formation is bursty and distributed throughout the galactic disk so as to prevent the formation of a concentrated star-forming region in the central galaxy core (this would be required to drive a galactic-scale outflow). In terms of CR feedback effects, this would represent a ‘best of both worlds’ scenario with CRs being important both inside and outside their host galaxy when large-scale outflows begin to emerge.

In protogalaxies, CRs have a greater ability to couple with their ambient medium than radiation through mechanisms not open to photons, e.g. hadronic processes or magnetic scattering. We have seen that one of the implications of these advected CRs includes the elevated external heating effect calculated in this paper, up to a level of around $10^{-34}~{}\text{erg}~{}\text{cm}^{-3}~{}\text{s}^{-1}$, even when several of kpc away from their source galaxy. This opens up new questions about the impacts of such an effect on, e.g. pre-heating for cosmic reionisation (Sazonov & Sunyaev, 2015; Leite et al., 2017) and the ability for such a process to be maintained. Moreover, advected CRs may have a role in amplifying intra-cluster and intergalactic magnetic fields, e.g. through resistive generation (see Miniati & Bell, 2011, 2012; Beck et al., 2013; Lacki, 2015, for application of this process to escaping CRs from high-redshift protogalaxies and clusters at the cosmic dawn), and/or by driving the growth of non-resonant magnetohydrodynamical instabilities in weak, pre-existing seed magnetic fields (see Bell, 2004; Miniati & Bell, 2011; D’Angelo, M. et al., 2015; Samui et al., 2018). There may, however, be even more important impacts than this.

Clusters, proto-clusters, and groups/pairs of protogalaxies are supported by outward pressure gradients against their gravitational potentials, essentially in hydrostatic equilibrium (see, e.g. Suto et al., 2013; Biffi et al., 2016). In low-redshift systems this is thought to be dominated by gas pressure from hot, thermal intergalactic baryons. However, at higher redshifts (above $z\gtrsim 1$, see Lacki 2015) when star-formation rates (and hence SN rates) in the Universe were greater than in the current epoch, and substantial CRs may have been able to escape from their host environments, CR pressure could be begin to dominate (see, e.g. Ginzburg & Syrovatskii, 1963; Lacki, 2015; Butsky & Quinn, 2018). In a pair of protogalaxies, this could have the effect of pushing neighbouring galaxies away, and this could alter the distribution of highly star-forming protogalaxies hosting outflows (e.g. perhaps those observed as Lyman-$\alpha$ emitters) compared to their more quiescent counterparts. A CR-advecting outflow may also have the ability to heat external gas and thus reduce (or even halt) gas inflow to the host galaxy. In turn, this may hamper star-formation by depriving the system of the cool, inflowing gases potentially responsible for driving the starburst phase (Kereš et al., 2005; Dayal et al., 2013; Lu et al., 2015; Toyouchi & Chiba, 2015; Yabe et al., 2015; Dayal & Ferrara, 2018), and this would lead to quenching or a reduction of the star-formation rate of the system by ‘strangulation’ (e.g. Peng et al., 2015) – see also the following section 5, where we demonstrate the impact CR containment and subsequent outflow advection may have in the recently observed high-redshift galaxy MACS1149-JD1 (Zheng et al., 2012; Hashimoto et al., 2018).

5 Application to the High-Redshift Starburst Galaxy MACS1149-JD1

A spectroscopic analysis of the doubly-ionised oxygen forbidden O III line, at 88$\mu$m131313This line is found to be particularly strong in chemically un-evolved environments, being a useful tracer of star-formation in metal-poor young star-forming galaxies (e.g. De Looze et al., 2014; Cormier et al., 2015) including those observed at high-redshift – see Inoue et al. (2016); Rigopoulou et al. (2018)., of the high-redshift galaxy MACS1149-JD1 was conducted by Hashimoto et al. (2018), obtaining a redshift of $z=9.11$. This corresponds to a time when the Universe was around 550 Myr old (Wright, 2006; Hashimoto et al., 2018). From the O III emission line, Hashimoto et al. (2018) also deduced a star-formation rate of $\mathcal{R}_{\rm SF}=4.2^{+0.8}_{-1.1}~{}\text{M}_{\odot}~{}\text{yr}^{-1}$ at the observational epoch. The wider spectrum and the Balmer break revealed an additional, older population of stars reminiscent of strong starburst activity a few hundred Myr earlier, implying MACS1149-JD1 had experienced two distinguishable star-formation episodes by the time it was observed. A two-phase star-formation history has also been found in other nearer galaxies (Dressler & Gunn, 1983; Couch & Sharples, 1987; French et al., 2015). It requires a mechanism to rapidly switch the galaxy from a strong starburst state into a quiescent state. The quenching of starbursts in galaxies is believed to be caused by feedback, which halts the gas inflow and/or gravitational collapse of over-dense gases (e.g. via turbulence or heating of the ambient gas, as suggested in French et al. 2018). Given that abundant molecular gas reservoirs, which are able to fuel star-formation, are found to be present in over half of quenched post-starburst systems studied, a lack of dense gas cannot account for the quenching. Instead, star-formation must be prevented by some other mechanism (French et al., 2015; Rowlands et al., 2015; Alatalo et al., 2016).

Hashimoto et al. (2018) presented three scenarios for the star-formation history of MACS1149-JD1 to account for the observed spectra. These begin respectively at $z=17.0,15.4$ and $13.3$ with duration and rates corresponding to a total stellar mass of the galaxy of $1.1^{+0.5}_{-0.2}\times 10^{9}~{}\text{M}_{\odot}$ by $z=9.11$. All scenarios allow a quiescent period of at least 100 Myr between the end of the first starburst phase and the beginning of the later star-formation phase observed at $z=9.11$. The scenarios are distinguished by their star-formation rate and the duration of their star-formation episode: (1) a slow star-formation rate $\mathcal{R}_{\rm SF}\approx 8~{}\text{M}_{\odot}~{}\text{yr}^{-1}$ for a duration of 200 Myr; (2) a moderate star-formation rate $\mathcal{R}_{\rm SF}\approx 16~{}\text{M}_{\odot}~{}\text{yr}^{-1}$ for a duration of 100 Myr; and (3) a shorter, but more intense star-formation rate of $\mathcal{R}_{\rm SF}\approx 160~{}\text{M}_{\odot}~{}\text{yr}^{-1}$ for a duration of 10 Myr. The termination time of the star-formation episode was not firmly specified.

5.1 Termination of Star Formation

We propose that the first starburst episode revealed by the Hashimoto et al. (2018) observation was terminated by progressive CR heating, causing a feedback process. The rapid evolution of the massive metal-poor stars gave rise to SNe (see Abel et al., 2002; Bromm et al., 2002; Clark et al., 2011), which produced shocks and injected turbulence into the ISM, and allowed for the acceleration of CRs. The shocks and turbulence also facilitated the growth of the galactic magnetic field via a turbulent dynamo mechanism (Balsara et al., 2004; Balsara & Kim, 2005; Beck et al., 2012; Schober et al., 2013) in MACS1149-JD1. Initially, the magnetic field was not strong enough and energy carried by CRs was mainly transported to the IGM instead of deposited into the ISM. As the starburst progressed, magnetic field amplification continued, eventually reaching saturation. By that time, the magnetic field had not only attained suitable strength to trap sufficient amounts of CRs within the galaxy (despite a fraction being advected out), but could also sustain a prolonged period of CR heating of power above all other radiation sources of stellar origin (Owen et al., 2018). The sequence of these processes should operate on a timescale comparable to the duration of the star-formation episode, otherwise the star-formation would not have continued to proceed. This sets the requirement that the longest possible timescale would determine the progression of the heating process that eventually led to quenching. A heuristic argument points out that the timescale on which the magnetic field evolves to saturation and the dynamical response to the sustained CR heating would be comparable. Thus, the starburst episode would have two stages: (i) an initial magnetic field growth stage and (ii) a subsequent delayed quenching stage due to sustained CR heating. In the magnetic field growth stage, CRs are contained as the magnetic field increases in strength (even when it has not fully saturated). As such, a prolonged sustained CR heating effect is established. At this point, the infall of cold gas has not been shut down, so star formation proceeds. In turn, this continues to drive the magnetic field amplification.

The termination of star-formation can only arise when there is no remaining avenue by which it can be sustained at any appreciable rate. This is achieved by the internalised CR heating of the system, as well as the strangulation of the cold inflows (e.g. Peng et al., 2015). The internal CR heating stage commences when the magnetic field amplification nears its saturation strength, and the ISM is heated for a prolonged period. Although the internal thermodynamics of the galaxy are modified accordingly, the gas inflows fuelling star-formation will not cease instantaneously. To halt the inflows, we argue that CRs are advected by outflows resulting from bust-like star-formation throughout the host galaxy. If MACS1149-JD1 had a starburst phase star-formation rate of $\mathcal{R}_{\rm SF}\approx 16~{}\text{M}_{\odot}~{}\text{yr}^{-1}$, this would correspond to a SN event rate of $\mathcal{R}_{\rm SN}=0.1~{}\text{yr}^{-1}$ (see section 2.3), which would cause the magnetic field to saturate at $\mu$G levels within around 140 Myr (Schober et al., 2013; Owen et al., 2018). In section 4.3.2, we showed that the advected CR power would be of order $10^{40}~{}\text{erg}~{}\text{s}^{-1}$ for such a system, being typically around 10% of the available total CR luminosity from the SN events, and this would predominantly be injected into the circumgalactic environment where the CRs would preferentially interact with the high-density cold inflowing filaments to develop an appreciable CR heating power – possibly comparable in strength to that experienced in the ISM of the host galaxy (possibly as high as $10^{-25}~{}\text{erg}~{}\text{cm}^{-3}~{}\text{s}^{-1}$).

This additional heating effect would raise the temperature of the inflowing gases to a level where they cannot drive star formation effectively and begin to evaporate. We may use the virial theorem to estimate the timescale over which filamentary inflows may halted in this manner: the virial temperature is the temperature above which gravitational collapse of gas is halted and, presumably, any heating to temperatures above this level would lead to evaporation. If assuming inflows persist over filaments of lengths of up to 50 kpc (Dekel et al., 2009; Stewart et al., 2013; Goerdt & Ceverino, 2015; Dayal & Ferrara, 2018) and that they are the sole driver of the star-formation activity arising at a rate of $\mathcal{R}_{\rm SF}\approx 16~{}\text{M}_{\odot}~{}\text{yr}^{-1}$ with a 30% mass conversion efficiency (Turner et al., 2015; Meier et al., 2002; Behroozi & Silk, 2015; Sun & Furlanetto, 2016) with an inflow velocity of 400 km s*-1* (from the velocity offset of the Lyman-$\alpha$ line in MACS1149-JD1 – see Hashimoto et al. 2018), the steady-state mass of these inflows would be around $6.7\times 10^{7}~{}\text{M}_{\odot}$. If the typical diameter of these flows is similar to the galaxy which they feed, i.e. 1 kpc, this would give a virial temperature of $T_{\rm vir}\approx 5,000~{}\text{K}$ (Binney & Tremaine, 2008). If adopting a number density of $10~{}\text{cm}^{-3}$ for the filamentary inflows (i.e. comparable to the mean density of the ISM in the host galaxy), the estimated CR heating rate of $10^{-25}~{}\text{erg}~{}\text{cm}^{-3}~{}\text{s}^{-1}$ would suggest it would take of order only a few Myrs for the virial temperature of the inflows to be exceeded. At this point, the supply of gas to the galaxy would be strangulated, halting star-formation within a dynamical timescale $\tau_{\rm dyn}$ (being the time required for the system to respond to the strangulation and internal heating). Since the strangulation and ISM quenching timescales are comparatively short, the magnetic saturation and dynamical timescales alone specify the timescale over which star-formation would be quenched.

In Fig. 16 the three star-formation scenarios for the first starburst phase proposed by Hashimoto et al. (2018) are shown. Along with these, the corresponding evolution of the galactic magnetic field as in Owen et al. (2018), using the parameters obtained by Hashimoto et al. (2018) for MACS1149-JD1 and following the prescription of Schober et al. (2013) are also shown. The dynamical timescale for galaxies similar to MACS1149-JD1 is $\tau_{\rm dyn}\approx\sqrt{3\pi/(16\;\!G\rho)}\approx 23~{}\text{Myr}$141414This is calculated by assuming a characteristic ISM density of $10~{}\text{cm}^{-3}$ rather than using the dynamical mass in Hashimoto et al. (2018), which has large uncertainties through their lensing parameter $\mu$. Imposing the requirement of a delayed response to the quenching of star-formation after CR containment is attained (as described above) rules out the most intense starburst scenario among those proposed by Hashimoto et al. (2018). It also sets an upper limit for the star-formation rate of below $20~{}\text{M}_{\odot}~{}\text{yr}^{-1}$ and a latest limit of before 260 Myr (i.e. at redshift $z\approx 15.4$) for the initial star-formation episode in this galaxy.

5.2 Reinstatement of Star Forming Activity

The cause of the resumption of star-formation after a quiescent period of 100 Myr remains to be explained. One of the possibilities is the eventual cooling of the hot circumgalactic medium. If strong CR heating were sufficient to cause the evaporation of cold inflows, such filaments could only start to be re-instated once the CR emission for the galaxy is diminished. Presumably this would arise some time after the end of star-formation, once any outflow activity and substantial CR production had abated.

Prior to the end of star-formation, the large amount of advected CRs, ultraviolet (UV) radiation from the young stars and/or X-rays from the stars and their remnants in a starburst galaxy (see Hashimoto et al., 2018) could heat and ionise the circumgalactic medium and the IGM. Thus a hot, ionised bubble would be carved out around the galaxy. Without a supply of cold gas, star-formation in the galaxy would be quenched (Peng et al., 2015). Hashimoto et al. (2018) considered a uniform IGM and UV escape fraction of 20% and estimated that a Strömgren (1939) sphere up to a radius 0.4 Mpc from MACS1149-JD1 would be created by the stellar UV radiation alone. If the observed later star-formation episode were re-ignited by the inflow of cold, neutral IGM gas from outside this Strömgren sphere, a free-fall velocity of around $4,000~{}\text{km}~{}\text{s}^{-1}$ (i.e. about 1% of the speed of light) would be required in order to reach the galaxy in the required time (100 Myr). This speed is excessive, given that an inflow velocity of only a few hundred ${\rm km}~{}{\rm s}^{-1}$ was obtained from the measurement of a blueshift in the Lyman-$\alpha$ line compared to the O III 88-$\mu$m rest-frame (see also Dijkstra et al., 2006; Verhamme et al., 2015, for further details and applications of this technique to detect inflows).

However, if the gas within the Strömgren (1939) sphere around MACS1149-JD1 could be cooled sufficiently, its inflow would fuel the star-formation process after the quiescent period. The thermal free-free cooling timescale of a hot ionised gas is roughly given by

[TABLE]

where $n_{\rm e}$ is the electron number density and $T_{\rm e}$ is the electron temperature. Taking $n_{\rm e}=\langle n_{\rm H}\rangle\approx 10^{-2}\;\!{\rm cm}^{-3}$ (where $n_{\rm H}$ is the number density of the medium) the value inferred from the Strömgren (1939) sphere estimation in Hashimoto et al. (2018) and adopting an electron temperature $T_{\rm e}=10^{5}\;\!{\rm K}$ gives a cooling time equal to the duration of the quiescent period inferred for MACS1149-JD1. Although, in reality, $T_{\rm e}$ could be higher than $10^{5}\;\!{\rm K}$ in a gas photo-ionised by UV photons and/or keV X-rays, the cooling times would be shortened if the electron density increases. Ionised dense clumps and filaments embedded within the ionised circumgalactic medium or IGM could be cooled more quickly, and these could provide the gas reservoir required to refuel subsequent star-formation. While the Strömgren sphere prescription would hold in the low-density regions between filaments with the H ​II region still extending to Mpc scales, the extent of the H ​II region would be much reduced in the direction of the clumps and filaments. For an ionisation distance in the direction of the filaments of around 10% of the Strömgren radius of MACS1149-JD1, an inflow velocity of around 400 km s*-1*, consistent with the velocity inferred from the offset of the Lyman-$\alpha$ line, would be sufficient to account for the re-ignition of star-formation activity within 100 Myr.

6 Conclusions

The heating power derived from CRs undergoing hadronic interactions has the potential to exceed that due to radiative heating processes from stellar and diffuse X-ray emission by one or two orders of magnitude, with CR heating reaching powers of around $10^{-25}~{}\text{erg}~{}\text{cm}^{-3}~{}\text{s}^{-1}$ (when adopting a SN-rate $\mathcal{R}_{\rm SN}=0.1~{}\text{yr}^{-1}$). In this work, we have shown that this remains the case, even when a strong galactic outflow develops. A galactic outflow is able to remove around 10% of the CRs from the ISM overall in a model protogalaxy, although this advection of CRs is predominantly focussed within the cone of the outflow itself. This points towards a ‘Two-Zone’ picture: one region which is predominantly advective, where the CRs are transported by the bulk motion of the outflow wind in which they are entrained, and another in which CR transport is predominantly diffusive where CRs are contained and deposit much of their energy into their host galaxy’s ISM.

The enhanced heating power of CRs contained by a protogalactic magnetic field can have a range of important consequences for the future evolution of the host galaxy and its neighbours. This effect can particularly have implications for subsequent star formation and the initial mass function of stars in the host after CR containment, as well as an enhanced X-ray emission due to the radiative inverse Compton losses of secondary particles produced in CR interaction showers (see also Schober et al., 2015). The degree to which these effects arise and their relative importance can only be understood with more detailed modelling. While the advective transport of CRs modelled in the present paper does reduce some of the effects discussed above, we calculate that the heating rates and CR containment expected in the presence of strong galactic outflow activity are not vastly changed, and are only substantially reduced in the parts of the galaxy directly affected by the outflow wind. This means that the phenomenological picture painted above is largely unchanged.

Further to this, the level of CR heating observed and the fraction of CRs which may be transported by advection to heat the surroundings implies that the effects on the IGM in the vicinity of the host galaxy are non-negligible. This opens up new questions about the impacts of CR heating and ionisation on, e.g. pre-heating for cosmic reionisation (Sazonov & Sunyaev, 2015; Leite et al., 2017) or on amplifying intra-cluster and/or intergalactic magnetic fields (e.g. Miniati & Bell, 2011; Beck et al., 2013; Lacki, 2015), and the ability for such processes to be maintained. The influence of the advected CRs may be even more important than this. If around 10% of CRs are indeed able to escape from highly star-forming galaxies in the presence of outflows, the balance of hydrostatic equilibrium may be changed in the host galaxy’s surroundings due to the additional CR pressure this would introduce. This may be particularly important in larger scale groups of galaxies and clusters containing one or more actively star-forming protogalaxies if the additional CR pressure provided is non-negligible and maintained for a sufficient length of time.

Acknowledgements

ERO is supported by a UK Science and Technology Facilities Council PhD Studentship. XJ’s visit to UCL-MSSL was supported by a UCL MAPS Dean’s Summer Research Scholarship and a MSSL Summer Studentship. SC’s visit to UCL-MSSL was supported by a MSSL Summer Studentship. We thank Prof Chung-Ming Ko (NCU) for discussions regarding transport of cosmic rays, Prof Mat Page (UCL-MSSL) for SN rates in local starburst systems and Prof Daisuke Kawata (UCL-MSSL) for implications of cosmic-ray heating in galaxies. We thank the anonymous referee for their constructive and informative comments which led to refinements and improvements of our model and discussion. This research has made use of NASA’s Astrophysics Data Systems.

Appendix A Numerical Scheme for the Transport Equation

In § 3.4 we invoke a numerical scheme to solve the differential equations 56 and 3.4.2. Here we outline the numerical scheme used in each case in more detail.

A.1 Advection Regime

In this case, we aim to solve

[TABLE]

where $b(E,h)=-{{\rm d}E}/{{\rm d}t}$ is known analytically from the total cooling processes (see § 3.2) but, apart from adiabatic cooling, is negligible in the case of CR protons, and where $\partial b(E,h)/\partial E$ immediately follows from this analytically. We retain these terms in the following treatment so as to describe a general scheme applicable to other particles where these cooling terms may not be negligible.

The equation is discretised according to a numerical grid over 170 points in energy $E_{j}$ distributed linearly between $E_{0}=1~{}{\rm GeV}$ and $E_{\rm max}=10^{6}~{}{\rm GeV}$ and 10,000 points in position $h_{i}$ distributed linearly between $0.1~{}\text{kpc}$ and $100~{}\text{kpc}$, indexed by the notation $j$ and $i$ respectively. This gives the first-order difference equation

[TABLE]

Each calculated point $Z_{i,j}$ requires the adjacent points in the previous $h$-step over three energies. This allows the gradient $\partial Z/\partial E$ to be estimated, and the gradient at the central point may then be propagated forwards. In the ‘edge-cases’ we find it is sufficient to simply estimate the gradient from two of the three available points; the energy bins are of sufficient resolution that any small inaccuracies from taking fewer points to calculate the gradient at the edges of the grid are rapidly suppressed. Our tests with higher order difference schemes did not yield noticeably different results.

We use a 4th order Runge-Kutta (RK4) scheme (Press et al., 2007) with 5th order error estimation to arrive at the numerical solution at each grid point, subject to the boundary conditions discussed in the main part of the text.

A.2 Diffusion Regime

Here, we aim to solve the second order differential equation

[TABLE]

which can be done in a similar way to the advection equation above (again, subject to the earlier boundary conditions) and splitting into an ‘inner’ and ‘outer’ scheme. We may use the RK4 method to find the numerical solution ${\rm d}Z/{\rm d}h$ in the ‘inner’ scheme according to the difference equation

[TABLE]

and a further ‘outer’ step is required to arrive at the numerical result $Z$, with difference equation

[TABLE]

Again, tests with higher order difference schemes indicated the approach adopted here was adequate. We note that additional Neumann boundary conditions are needed for the inner numerical scheme, i.e. step 1 in equation 75, as detailed in the main text – see § 3.4.