Cosmic rays and non-thermal emission in simulated galaxies – III. Probing cosmic-ray calorimetry with radio spectra and the FIR–radio correlation

ABSTRACT

An extinction-free estimator of the star formation rate (SFR) of galaxies is critical for understanding the high-redshift universe. To this end, the nearly linear, tight correlation of far-infrared (FIR), and radio luminosity of star-forming galaxies is widely used. While the FIR is linked to massive star formation, which also generates shock-accelerated cosmic-ray (CR) electrons and radio synchrotron emission, a detailed understanding of the underlying physics is still lacking. Hence, we perform three-dimensional magnetohydrodynamical (MHD) simulations of isolated galaxies over a broad range of halo masses and SFRs using the moving-mesh code arepo, and evolve the CR proton energy density self-consistently. In post-processing, we calculate the steady-state spectra of primary, shock-accelerated and secondary CR electrons, which result from hadronic CR proton interactions with the interstellar medium. The resulting total radio luminosities correlate with the FIR luminosities as observed and are dominated by primary CR electrons if we account for anisotropic CR diffusion. The increasing contribution of secondary emission up to 30 per cent in starbursts is compensated by the larger bremsstrahlung and Coulomb losses. CR electrons are in the calorimetric limit and lose most of their energy through inverse Compton interactions with star light and cosmic microwave background (CMB) photons while less energy is converted into synchrotron emission. This implies steep steady-state synchrotron spectra in starbursts. Interestingly, we find that thermal free–free emission flattens the total radio spectra at high radio frequencies and reconciles calorimetric theory with observations while free–free absorption explains the observed low-frequency flattening towards the central regions of starbursts.

MHD, methods: numerical, cosmic rays, galaxies: magnetic fields, galaxies: starburst, radio continuum: galaxies

1 INTRODUCTION

The radio emission from star-forming galaxies is attributed to the interaction of CR electrons with an ambient interstellar magnetic field. These highly relativistic charged particles emit synchrotron radiation while gyrating around magnetic field lines, giving rise to a power-law spectrum at radio frequencies. Because CRs are accelerated at shocks of supernova (SN) remnants and because the lifetime of radio-emitting CR electrons in galaxies is shorter than tens of millions of years, their existence reveals ongoing star formation in a galaxy. Hence, a correlation between the radio synchrotron emission of star-forming galaxies and tracers of their star formation activity, such as the FIR luminosity, is expected. The latter traces ongoing star formation, because the ultraviolet (UV) light of a young stellar population is absorbed by dust and re-emitted in the FIR.

Indeed, a tight and nearly linear correlation has been observed between the radio and FIR luminosity of galaxies (van der Kruit 1971; Helou, Soifer & Rowan-Robinson 1985; Condon 1992; Yun, Reddy & Condon 2001; Bell 2003; Matthews et al. 2021; Molnár et al. 2021). Due to the tightness of this relation across different types of star-forming galaxies and across five orders of magnitude in luminosity, the FIR–radio correlation (FRC) is widely used to estimate the SFR of galaxies from their radio luminosities (see e.g. Heesen et al. 2014, 2019; Vollmer et al. 2020). It has the notable feature of not being affected by dust extinction, which is one of the main uncertainties of other SFR estimators, which employ, e.g. the Hα and/or UV emission. This makes the FRC a favourable method for estimating SFRs in particular for high redshifts galaxies, where dust properties are unknown, and where accurate SFR estimators are critical for deciphering the cosmic star formation history. However, it is still unclear, whether the FRC evolves with redshift (Lacki & Thompson 2010; Schleicher & Beck 2013; Schober, Schleicher & Klessen 2016).

A linear FRC naturally emerges if CR electrons calorimetrically lose most of their energy to synchrotron emission (Voelk 1989; Lisenfeld, Voelk & Xu 1996). This requires that the corresponding time-scale of synchrotron losses is shorter than all other loss processes that apply to CR electrons, i.e. losses due to inverse Compton (IC) and bremsstrahlung emission, as well as Coulomb losses and CR electron escape from the galaxy. The latter process arises due to advection and diffusion, and has been suggested to be relevant in galaxies with low SFRs (Thompson et al. 2006; Lacki, Thompson & Quataert 2010). But highly star-forming galaxies, such as NGC 253 and M82, are expected to be calorimetric, due to high photon energy densities and magnetic field strengths. However, in this picture, the fully cooled electron spectra would steepen by unity at high particle energies due to the energy dependence of IC and synchrotron cooling. This would imply steep radio spectra with spectral indices at around 1.1, which is in contrast to observed values of 0.5–0.8. Hence, these flat radio spectra of star-forming galaxies pose a challenge to the calorimetric theory.

To gain insights into these processes, the radio and gamma-ray spectra of famous starburst galaxies such as M82 or NGC 253 have been analyzed with so-called one-zone models in which a galaxy is represented as a single entity and parametrized by average quantities (Torres 2004; Lacki et al. 2010; Paglione & Abrahams 2012; Yoast-Hull et al. 2013; Eichmann & Becker Tjus 2016). One proposed solution for the problem of too steep radio spectra predicted from calorimetric theory has been bremsstrahlung and ionization losses in starburst galaxies (Thompson et al. 2006; Lacki et al. 2010; Basu et al. 2015), which could flatten CR electron spectra due to their shallower energy dependence in comparison to IC and synchrotron losses. However, at the same time, these additional processes open new energy loss channels for the radio-emitting CR electrons and thus diminish the resulting radio synchrotron luminosity. This would imply a sublinear FRC at high SFRs, which is in conflict with observations. To cure this problem, a second radio emission process has been proposed that fills in this missing radio emission in form of an increasing contribution of synchrotron emission from secondary electrons that are generated in hadronic CR proton interactions with the ambient interstellar medium (ISM, Thompson et al. 2006; Lacki et al. 2010), an efficient process in the dense centres of starbursts. This would successfully maintain the linear FRC, even at high SFRs. On the other hand, IC losses also compete with synchrotron losses that in turn might pose a challenge to use the FRC as SFR estimator, in particular at high redshifts, where IC losses due to scattering of CR electrons off CMB photons are expected to be larger. Thus, in addition to the question of calorimetry, understanding the relative importance of IC and synchrotron losses is essential.

To model these processes in more detail, we perform MHD simulations with the moving-mesh code arepo, simulate the CR proton energy density and model the non-thermal CR spectra. This is the third and last paper of a series in which we present our results. The first paper (Werhahn et al. 2021a, hereafter Paper I) details the steady-state modelling of CRs in our simulations and compares the results to CR data. We find very good agreement for our CR electron and proton spectra as well as the positron fraction as a function of particle energy when compared to Voyager-1 and AMS-02 data. The resulting gamma-ray emission of star-forming galaxies is analysed in the second paper (Werhahn et al. 2021b, hereafter Paper II), where we successfully match the observed FIR–gamma-ray relation as well as observed gamma-ray spectra of star-forming galaxies. Here, we apply the same formalism as described in Paper I and study the resulting radio synchrotron emission from primary and secondary CR electrons. In a companion paper, Pfrommer et al. (2021) investigate the origin of the global FRC. In particular, the role of the turbulent small-scale dynamo is identified and the processes leading to the observed scatter in the FRC are analyzed.

This paper is structured in the following way. In Section 2, we describe our simulations of isolated galaxies, the steady-state modelling of the CR proton, primary and secondary electron spectra, as well as the calculation of the resulting radio emission and absorption processes. Section 3 provides an overview of the time-scales of the non-radiative and radiative losses of CR protons and electrons. We furthermore describe the modelling of the FRC, in Section 4, where we analyse the primary and secondary contribution to the total radio luminosity, quantify the calorimetry of CR electrons and assess the processes that have been proposed to ‘conspire’ to maintain an almost linear FRC at high SFRs. Eventually, we scrutinize three possible mechanisms, that could be responsible for the observed flat radio spectra of star-forming galaxies such as NGC 253, M82, and NGC 2146 in Section 5, before we summarize our results in Section 6.

2 DESCRIPTION OF THE METHODS

2.1 Simulations

We perform simulations with the moving mesh code arepo(Springel 2010; Pakmor et al. 2016), as described in Paper I and Paper II. Our simulations start from the gravitational collapse of a gas cloud embedded in a dark matter (DM) halo prescribed by an NFW (Navarro, Frenk & White 1997) profile, with masses ranging from M₂₀₀ = 10¹⁰ to |$10^{12}\, \mathrm{M_{\odot }}$|⁠, in order to resemble realistic halo masses from dwarf to Milky-Way (MW) sized galaxies. Subsequently, we simulate the formation of a rotationally supported disc that forms stars stochastically above a critical density threshold (Springel & Hernquist 2003). After an initial burst of star formation, the SFR decreases exponentially over time. At SNe, CRs are instantaneously injected with an energy fraction ζ_SN of the kinetic energy of the SN explosion, which we vary from 5 to 10 per cent. The lower range of ζ_SN is inferred from a combination of kinetic plasma simulations at oblique shocks (Caprioli & Spitkovsky 2014) and three-dimensional MHD simulations of CR acceleration at SN remnant shocks (Pais et al. 2018), which is followed by a detailed comparison of simulated radio, X-ray and gamma-ray emission maps and spectra to observational data (Pais et al. 2020; Pais & Pfrommer 2020; Winner et al. 2020).

The ideal MHD prescription (Pakmor & Springel 2013) governs the evolution of the magnetic field. Starting from an initial seed magnetic field B₀ permeating the gas cloud before the collapse, the magnetic field is exponentially amplified by a small-scale dynamo, before it grows further via an inverse cascade until saturation (Pfrommer et al. 2021). We chose two different values for the initial magnetic field of B₀ = 10⁻¹⁰ and |$10^{-12}\, \mathrm{G}$|⁠, which represent the pre-amplified magnetic field in a protogalactic environment. CRs are described as a relativistic fluid (Pakmor et al. 2016; Pfrommer et al. 2017a) with adiabatic index of 4/3 and we account for adiabatic changes of the CR energy density as well as Coulomb and hadronic loses due to interactions with the ISM. Furthermore, CRs are advected with the gas (in our ‘CR adv’ models) and are additionally allowed to anisotropically diffuse along magnetic field lines (in our ‘CR diff’ models) with a parallel diffusion coefficient along magnetic field lines of |$D_0=10^{28}\mathrm{cm^2\, s^{-1}}$|⁠.

A summary of all simulations analysed in this paper is presented in Table 1. The concentration parameter of the NFW profile is given by c₂₀₀ = r₂₀₀/r_s = 12 in all simulations, where the characteristic scale radius of the NFW profile is denoted by r_s and r₂₀₀ is the radius enclosing a mean density that is equal to 200 times the critical density of the universe. We refer to the simulation in the ‘CR diff’ model with a halo mass of |$M_{200}=10^{12}\, \mathrm{M_\odot }$|⁠, B₀ = 10⁻¹⁰ G, and ζ_SN = 0.05 as ‘fiducial halo’ in the following.

Table 1.

Overview of the simulations in this paper. Shown are (1) the halo mass M₂₀₀, (2) the CR transport model: In the ‘CR adv’ model, we only account for CR advection with the gas whereas the ‘CR diff’ model additionally allows for anisotropic diffusion, (3) the initial magnetic field B₀, (4) the injection efficiency of CRs at SNRs, ζ_SN, and (5) the referenced figures.

M200 (M⊙)	CR model	B0 (G)	ζSN	Figures
1010	CR adv	10−10	0.05	3, 6, 7, 12
1011	CR adv/CR diff	10−10	0.05	3, 6, 7, 12
3 × 1011	CR adv/CR diff	10−10	0.05	3, 6, 7, 12
1012	CR adv/CR diff	10−10	0.05	1–13
1010	CR adv	10−12	0.05	A2
1011	CR adv/CR diff	10−12	0.05	A2
3 × 1011	CR adv/CR diff	10−12	0.05	A2
1012	CR adv/CR diff	10−12	0.05	A2
1010	CR adv	10−12	0.10	A2
1011	CR adv/CR diff	10−12	0.10	A2
3 × 1011	CR adv/CR diff	10−12	0.10	A2
1012	CR adv/CR diff	10−12	0.10	A2

2.2 CR steady-state spectra

In addition to the primary CR electron population that is accelerated at SNRs together with CR protons, hadronic interactions of CR protons with the ISM lead to the production of neutral and charged pions. Pions decay further into muons and eventually to neutrinos and secondary electrons and positrons (hereafter referred to as secondary electrons).

In order to model the spectral distribution of CRs, we assume a steady state and solve the diffusion-loss equation (Ginzburg & Syrovatskii 1964; Torres 2004) in each cell for protons, primary and secondary electrons (see Paper I), which reads

$$\begin{eqnarray*} \frac{\mathrm{}f_i(E_i)}{\tau _{\mathrm{esc}}}-\frac{\mathrm{d}}{\mathrm{d}E_i}\left[f_i(E_i)b(E_i)\right]=q_i(E_i). \end{eqnarray*}$$

(1)

Here, E_i denotes the CR energy, the subscript i = e, p specifies the CR species and f_i(E_i) = dN_i/(dVdE_i) is the resulting equilibrium spectral density for either protons (denoted by p), primary or secondary electrons (denoted by e). We define the source function q_i(E_i) = q_i(p_i)dp_i/dE_i in terms of a power law in momentum with an exponential cut-off given by

$$\begin{eqnarray*} q_{i}(p_{i})\mathrm{d}p_{i} = C_{i} p_{i}^{-\alpha _{\mathrm{inj}}} \exp [-(p_i/p_{\mathrm{cut},i})^{n}]\mathrm{d}p_{i}, \end{eqnarray*}$$

(2)

where p_i = P_i/(m_ic) denote normalized momenta, m_i is the proton/electron rest mass, and c is the speed of light. We adapt n = 1 for protons and n = 2 for primary electrons (Zirakashvili & Aharonian 2007; Blasi 2010) in equation (2) and assume cutoff momenta for protons |$p_{\mathrm{cut,p}}=1\, \mathrm{PeV}/(m_{\mathrm{p}}c^2)$| (Gaisser 1990) and for electrons |$p_{\mathrm{cut,e}}=20\, \mathrm{TeV}/(m_{\mathrm{e}}c^2)$| (Vink 2012). Both, primary electrons and protons are injected with a spectral index of α_inj = 2.2. We discuss the source function of secondary electrons in Section 2.3.

The resulting steady-state distribution of CR protons is ensured to match the CR energy density in each cell, after all cooling processes have been taken into account. The energy loss processes are given by the cooling rate b(E_i) = −dE_i/dt, that comprise hadronic and Coulomb losses in the case of CR protons. Additionally, losses due to particle escape are quantified by the time-scale |$\tau _{\mathrm{esc}}=(\tau _{\mathrm{adv}}^{-1}+\tau _{\mathrm{diff}}^{-1})^{-1}$|⁠, where advection and diffusion losses are included. These are estimated via

$$\begin{eqnarray*} \tau _{\mathrm{diff}}=\frac{L_{\mathrm{CR}}^{2}}{D} \end{eqnarray*}$$

(3)

and

$$\begin{eqnarray*} \tau _{\mathrm{adv}}=\frac{L_{\mathrm{CR}}}{v _{\mathrm{z}}}. \end{eqnarray*}$$

(4)

Here, L_CR = ε_CR/|∇ε_CR| is an estimate for the diffusion length. For diffusion losses, we assume an energy dependence of the diffusion coefficient D = D₀(E/E₀)^δ, where |$D_0=10^{28}\mathrm{cm^2\, s^{-1}}$|⁠, |$E_0=3\, \mathrm{GeV}$|⁠, and δ = 0.5, which can be inferred from observations of beryllium isotope ratios (Evoli et al. 2020). However, we find in Paper II, that gamma-ray spectra of highly star-forming galaxies provide a better match for a shallower energy dependence of the diffusion coefficient with δ = 0.3 and hence, we also adopt both values of δ here. For advection losses, only the verical velocity component |$v _z$| pointing away from the disc is taken into account, because we show in Paper I, that azimuthal fluxes in and out of a computational cell nearly compensate each other, so that advection predominantly happens in vertical direction.

The assumption of a cell-based steady state has been analyzed in Paper I, where we compare the time-scale of the change in the total energy density of CRs over a global simulation time-step τ_CR to the loss time-scale that includes cooling and diffusion losses τ_all. We find that the requirement for a steady-state configuration, i.e. τ_all/τ_CR < 1 is reached in cells that contribute most to the non-thermal emission processes, i.e. both in the radio and gamma-ray regime. Hence, we consider this a good assumption for our study concerning the non-thermal emission processes. However, because the steady-state assumption breaks down in low-density regions, in outflows and near SNRs, we aim towards a more accurate treatment of the time evolution of CR protons and electrons (Winner et al. 2019; Girichidis et al. 2020) in future work.

2.3 Primary and secondary CR electrons

To obtain the steady-state distribution f_e of primary CR electrons, we adapt equation (2) as the source function and solve the steady-state equation, i.e. equation  (1). The normalization of f_e is set by requiring the primary CR electrons to reproduce the observed ratio of electrons to protons, |$K_{\mathrm{ep}}^{\mathrm{obs}}$|⁠, at 10 GeV in a snapshot of our simulations resembling the MW in terms of halo mass and SFR. To this end, we average over the cell-based steady-state spectra around the solar galactocentric radius and re-normalize the spectra according to the observed electron to proton ratio. From this, we can infer the injected ratio of primary electrons to protons before cooling, |$K_{\mathrm{ep}}^{\mathrm{inj}}\approx 0.02$|⁠, that is needed, in order to obtain the observed value of |$K_{\mathrm{ep}}^{\mathrm{obs}}\approx 0.01$| at 10 GeV (Cummings et al. 2016) after cooling. We assume that the injected ratio of electrons to protons is universal and adapt it to the rest of the galaxy as well as all other simulated galaxies.

To calculate the production of secondary CR electrons and positrons, we adopt the parametrization of Kelner, Aharonian & Bugayov (2006) for large kinetic proton energies |$T_{\mathrm{p}}\gt 100\, \mathrm{GeV}$|⁠. Because Yang, Kafexhiu & Aharonian (2018) provide a more detailed modelling of the differential cross-section of pion production near the threshold of pion production up to |$10\, \mathrm{GeV}$|⁠, we use their parametrization for small proton energies and perform a cubic spline interpolation in the intermediate energy range, combining it with our own fit of the total cross-section of pion production to the data (see equations B1, B5, and B6 in Paper I).

In addition to Coulomb and escape losses, CR electrons also suffer losses due to radiative processes, including synchrotron, bremsstrahlung and IC emission (see Section 3). In order to model the latter, we assume that the energy density of the incident radiation field is composed of the CMB and stellar radiation, i.e. ε_ph = ε_CMB + ε_⋆. The stellar contribution can be approximated by the FIR emission that results from dust-reprocessed UV light from young stellar populations, which we approximate with a blackbody distribution characterized by a temperature T = 20 K (Calzetti et al. 2000). The resulting photon energy density in each cell is calculated as

$$\begin{eqnarray*} \varepsilon _{\star }=\sum _{i}\frac{L_{\mathrm{FIR},i}}{4 \pi R_{i}^{2}c}, \end{eqnarray*}$$

(5)

where we infer the FIR–luminosity L_{FIR, i} from the SFR in each cell (Kennicutt 1998) and sum over the fluxes arriving from all other cells at a distance of R_i by using a tree-code. If a cell is itself actively star forming, we equate the distance with the cells’ radius, which we derive from its volume V_i via R_i = [3V_i/(4|$\pi$|⁠)]^1/3.

Of particular interest, in this work, is the relative importance of radio emission due to primary versus secondary CR electrons. The ratio of their distribution functions at 10 GeV has been derived in Paper I using a simplified analytical approximation:

$$\begin{eqnarray*} \frac{f_{\mathrm{e}}^{\mathrm{prim}}}{f_{\mathrm{e}}^{\mathrm{sec}}} & =K_{\mathrm{ep}}^{\mathrm{inj}}\frac{3}{128}16^{\alpha _{\mathrm{p}}}\left(1+\frac{\tau _{\pi }}{\tau _{\mathrm{esc}}}\right), \end{eqnarray*}$$

(6)

where |$\tau_{\pi}$| denotes the time-scale of pion production via hadronic CR proton interactions with the ISM, α_p is the spectral index of the CR proton distribution, and |$f_{\mathrm{e}}^{\mathrm{sec}} = f_{\mathrm{e}^+}^{\mathrm{sec}} + f_{\mathrm{e}^-}^{\mathrm{sec}}$| is the total steady-state distribution of secondary electrons and positrons. Adapting |$K_{\mathrm{ep}}^{\mathrm{inj}}=0.02$| and calorimetric conditions for secondary electron production, i.e. τ_|$\pi$| ≪ τ_esc (⁠|$\tau_{\pi}$| ≃ τ_esc), which consequently implies α_p = α_inj = 2.2, equation (6) yields |$f_{\mathrm{e}}^{\mathrm{prim}}/f_{\mathrm{e}}^{\mathrm{sec}}\approx 0.2$| (0.4). In fact, this can be interpreted as a lower limit, because as soon as escape losses become more important than hadronic losses, i.e. |$\tau_{\pi}$| ≳ τ_esc, secondary production ceases to be efficient and primary electrons will dominate over secondary electrons and positrons. Thus, we expect for simulations that only account for CR advection and neglect CR diffusion secondary electrons to dominate by a factor of ∼5 at 10 GeV in the calorimetric limit. By contrast, if we also include energy-dependent diffusion, we will obtain steeper CR proton spectra with α_p ≳ 2.2 as soon as diffusion losses become important. This increases the ratio of primary to secondaries at 10 GeV significantly and we can expect primary electrons to dominate over secondaries.

2.4 Radio emission and absorption processes

Radio emission from star-forming galaxies is due to non-thermal synchrotron as well as thermal bremsstrahlung emission, i.e free–free emission. The modelling of these emission processes, together with the corresponding absorption processes, is summarized in the following.

2.4.1 Radio synchrotron emission

Using our steady-state CR electron population f_e, we calculate the radio synchrotron emissivity, |$j_{\nu }=E\, \mathrm{d}N_{\gamma }/(\mathrm{d}\nu \mathrm{d}V \mathrm{d}t)$|⁠, in each computational cell following Rybicki & Lightman (1986), via

$$\begin{eqnarray*} j_{\nu }=\frac{\sqrt{3}e^{3}B_{\perp }}{m_{\rm e}c^{2}}\intop _{0}^{\infty }f_{\mathrm{e}}(p_{\mathrm{\rm e}})F(\nu /\nu _\rm {c})\mathrm{d}p_{\mathrm{e}}, \end{eqnarray*}$$

(7)

where B_⊥ is the component of the magnetic field perpendicular to the line of sight, e denotes the elementary charge, and we use an analytical approximation provided by Aharonian, Kelner & Prosekin (2010) for the dimensionless synchrotron kernel |$F(\nu /\nu _\rm {c})$| (equation A2), where |$\nu _\rm {c}$| is the critical frequency (see Appendix A1 for details). The typical synchrotron emission frequency ν_syn is related to the electron Lorentz factor γ_e via¹

$$\begin{eqnarray*} \nu _\rm {syn}=2 \nu _c =\frac{3eB_\perp }{2\pi \, m_{\rm {e}}c}\, \gamma _\rm {e}^2\simeq 1~\rm {GHz}\, \frac{B_\perp }{1\, \mu \mbox{G}}\, \left(\frac{\gamma _\mathrm{e}}{10^4}\right)^2. \end{eqnarray*}$$

(8)

This indicates that observations of synchrotron emission at a fixed frequency typically probe electrons with lower (higher) Lorentz factors in higher (lower) magnetic field strengths, which we refer to as the ‘ν_c-effect’ in the following. The specific radio synchrotron intensity I_ν at frequency ν (in units of |$\mathrm{erg\, cm^{-2}\, s^{-1}\, Hz^{-1}\, sterad^{-1}}$|⁠) is obtained by integrating j_ν along the line of sight (denoted by s) and the specific radio luminosity (in units of |$\mathrm{erg\, s^{-1}\, Hz^{-1}}$|⁠) is obtained by integrating j_ν over the entire galaxy, yielding

$$\begin{eqnarray*} I_\nu = \frac{1}{4\pi } \intop _0^{\infty } j_\nu \mathrm{d}s \qquad \mbox{and}\qquad L_{\mathrm{\nu }}=\intop j_\nu \mathrm{d}V. \end{eqnarray*}$$

(9)

2.4.2 Thermal free–free emission and absorption processes

In addition to non-thermal synchrotron emission, we also expect a contribution from thermal free–free emission in the radio band, which predominantly depends on the electron density and temperature. To model this accurately, a multiphase model of the ISM with radiative transfer would be required. Since this is beyond the scope of this work, we adapt a simplified model here, where we assume a temperature of T = 8000 K for the warm-ionized medium, which dominates this emission component. To estimate the corresponding electron density, we take a fixed fraction |$\xi _\rm {e}$| of the electron density |$n_\mathrm{e}=x_\mathrm{e}\, n_\mathrm{H}$| that we calculate from the fraction of free electrons x_e provided by our simplified pressurized ISM model (Springel & Hernquist 2003), that models the ISM with an effective equation of state, and leave |$\xi _\rm {e}$| as a free parameter of order unity.

In addition, we account for free–free absorption and synchrotron self-absorption (SSA) by means of the radiative transfer equation (see equation A13). To this end, we (i) rotate the simulation into a desired inclination of the disc with the line of sight, (ii) construct thinly spaced slices through the simulation box perpendicular to the line of sight, from the front of the simulation volume to its back, and (iii) cumulatively add the optical depth, that is a sum of the optical depth of free–free absorption, τ_ff, and SSA, τ_SSA. Because the electron density along the line of sight depends on galactic inclination, the resulting amount of free–free emission and absorption inherits this dependence and so does the shape of the radio spectrum. Whereas free–free absorption predominantly affects the low-frequency part of radio spectra (see Section 5.1), we find that SSA has a negligible effect at radio frequencies studied here.

3 TIME-SCALES

This section will provide an overview of the time-scales of all relevant processes, which we define as |$\tau =E/\dot{E}$|⁠. Because we aim at explaining the FRC at 1.4 GHz, we compute the time-scales at a fixed observational frequency. According to equation (8), electrons that emit synchrotron radiation at GHz frequencies in |$\mu$|G magnetic fields have a typical energy of |$E_{\mathrm{e}}=10^4\, m_{\mathrm{e}}c^2\approx 5$| GeV. These electrons can be primary and secondary electrons while the latter have been hadronically created by CR protons with typical energies of E_p ≈ 16E_e ≈ 80 GeV.

3.1 Non-radiative processes

The injection time-scale of secondaries corresponds to the time-scale of hadronic interactions of CR protons with the ISM, and is given by

$$\begin{eqnarray*} \tau _{\pi } = \frac{1}{c n_{\mathrm{N}}K_{\mathrm{p}}\sigma _{\mathrm{pp}}} \approx 240 \, \left(\frac{n_{\mathrm{N}}}{0.1~\rm {cm}^{-3}}\right)^{-1}\mathrm{Myr}, \end{eqnarray*}$$

(10)

where the cross-section of hadronic interactions is |$\sigma _{\mathrm{pp}}\approx 44 \, \mathrm{mbarn}$| for α_p = 2.2 (Pfrommer et al. 2017b), the inelasticity of hadronic interactions is given by K_p = 1/2 (Mannheim & Schlickeiser 1994), n_N = n_H + 4n_He = (X_H + 1 − X_H)ρ/m_p = ρ/m_p is the number density of target nucleons in the ISM, where X_H = 0.76 denotes the hydrogen fraction and ρ is the gas density. Consequently, the hadronic time-scale directly traces the gas density, as can bee seen in the upper right-hand panels in Fig. 1, which shows |$\tau _{\pi}$| for our M82-like galaxy, see Table 2.

Figure 1.

The upper panels show from the left- to right-hand maps (slices) of the characteristic cooling time-scales of synchrotron and IC cooling for CR electrons, as well as the time-scale of hadronic interactions of CR protons, that is relevant for the production of secondary electrons. In the lower panels, we show the ratios of IC to synchrotron cooling, as well as the ratio of escape to synchrotron and to IC losses, respectively, where escape losses include losses due to CR advection and diffusion. The time-scale ratios are averaged over thin slices with a thickness of 500 pc and the time-scales are calculated at an energy of 10 GeV. All maps are shown for a snapshot with a halo mass of |$M_{200}=10^{12}\, \mathrm{M_{\odot }}$| at t = 2.3 Gyr.

Coulomb interactions of CRs with the ambient medium of an electron number density n_e act on a time-scale (Gould 1972)

$$\begin{eqnarray*} \tau _\mathrm{Coul} &=& \frac{2E_\mathrm{e} \beta _{\mathrm{e}}}{{3\sigma _\rm {T}n_{\mathrm{e}}m_{\mathrm{e}}c^3}} \left[\ln \left(\frac{m_{\mathrm{e}}c^{2}\beta _{\mathrm{e}}\sqrt{\gamma _{\mathrm{e}}-1}}{\hbar \omega _{\mathrm{pl}}}\right)-\ln (2)\right.\nonumber \\ && \times \,\left.\left(\frac{\beta _{\mathrm{e}}^{2}}{2}+\frac{1}{\gamma _{\mathrm{e}}}\right)+\frac{1}{2}+\left(\frac{\gamma _{\mathrm{e}}-1}{4\gamma _{\mathrm{e}}}\right)^{2}\right]^{-1} \propto \frac{E_\mathrm{e}}{n_\mathrm{e}}\propto \frac{1}{B^{1/2} n_\mathrm{e}}, \nonumber \\ \end{eqnarray*}$$

(11)

where the normalized electron velocity is denoted by |$\beta _\mathrm{e}=v _\mathrm{e}/c$|⁠, σ_T is the Thompson cross-section, the plasma frequency is defined by |$\omega _\mathrm{pl}=\sqrt{4\pi e^{2}n_{\mathrm{e}}/m_{\mathrm{e}}}$| and we used in the last step the relativistic limit and equation (8) so that the last expression is only valid at fixed synchrotron emission frequency. Adapting an electron number density of |$n_\mathrm{e}=0.1\, \mathrm{cm^{-3}}$|⁠, Coulomb losses act on time-scales of ∼2 Gyr for highly relativistic electrons with γ_e = 5 GeV/(m_ec²) ≈ 10⁴, whereas middly relativistic electrons with γ_e ∼ 10 cool on time-scales of |$\tau \approx 3\, \mathrm{Myr}$|⁠. Thus, Coulomb losses typically remove the low-energy part of the CR electron spectrum. Consequently, in order for radio synchrotron-emitting electrons (with typical Lorentz factors γ_e ∼ 10⁴ for |$B=1~\mu$|G and |$\nu _\rm {syn}=1.4$| GHz) to be affected by Coulomb losses, very high densities are required.

3.2 Radiative processes

At high electron energies, radiative losses due to synchrotron, IC or bremsstrahlung emission are typically dominant. The loss time-scale due to synchrotron emission is given by

$$\begin{eqnarray*} \tau _\mathrm{syn} =\frac{3E_\mathrm{e}}{\sigma _{\mathrm{T}} c \beta _{\mathrm{e}}^{2} \gamma _{\mathrm{e}}^{2}\varepsilon _{B} } \propto \frac{1}{E_\mathrm{e} B^2}\propto B^{-3/2}, \end{eqnarray*}$$

(12)

where |$\varepsilon _B=B^2/(8\pi)$| and we used equation (8) in the last step so that this dependence is only valid at fixed synchrotron emission frequency. For instance, in central regions of starburst galaxies with magnetic field strengths of |$\approx 10\, \mu$|G, where electrons with typical energies of E_e ≈ 2 GeV are responsible for radio synchrotron emission at GHz frequencies, the synchrotron cooling time-scale is |$\tau _{\mathrm{syn}}\approx 70\, \mathrm{Myr}$|⁠. In starburst nuclei with up to |$B\approx 50~\mu$|G, synchrotron cooling acts on even shorter time-scales of |$\tau _\mathrm{syn}\approx 6\, \mathrm{Myr}$|⁠.

In an ambient radiation field with a photon energy density ε_ph, CR electrons scatter off of these photons and lose energy on a time-scale

$$\begin{eqnarray*} \tau _\mathrm{IC} =\frac{3E_\mathrm{e}}{\sigma _{\mathrm{T}} c \beta _{\mathrm{e}}^{2} \gamma _{\mathrm{e}}^{2}\varepsilon _{\mathrm{ph}} } \propto \frac{1}{E_\mathrm{e} \varepsilon _\mathrm{ph}}\propto \frac{B^{1/2}}{\varepsilon _\mathrm{ph}}, \end{eqnarray*}$$

(13)

where the last step is only valid at fixed synchrotron emission frequency. Note that ε_ph is a sum of a stellar contribution (ε_⋆, modelled by equation 5) and the CMB (⁠|$\varepsilon _{\mathrm{CMB}}\approx 4.16\times 10^{-13}\, (1+z)^2\, \mathrm{erg\, cm^{-3}}$|⁠, where z is the cosmic redshift). The CMB is ubiquitous and independent of the local SFR (or ε_⋆). For a photon energy density of ε_ph = 5ε_CMB, the IC cooling time-scale ranges from ∼40 to ∼140 Myr in regions with magnetic field strengths of 1 –|$10\, \mu$|G, where we again adapted equation (8). Because |$\tau _\mathrm{syn}/\tau _{\mathrm{IC}}\propto B_\mathrm{ph}^2/B^2$|⁠, where |$B_{\mathrm{ph}}=(8\pi \varepsilon _{\mathrm{ph}})^{1/2}$| is the equivalent magnetic field strength, IC losses will always dominate over synchrotron losses if B_ph > B. Because the photon energy density of the CMB corresponds to |$B_{\mathrm{CMB}}\approx 3~\mu \mathrm{G}$| today, IC losses will be relevant as soon as the magnetic field drops below |$3~\mu \mathrm{G}$|⁠, and correspondingly at larger magnetic field strengths if we additionally account for a stellar radiation field.

In a fully ionized medium with a proton number density n_p, CR electrons lose energy due to the emission of bremsstrahlung on a time-scale (Blumenthal & Gould 1970)

$$\begin{eqnarray*} \tau _\mathrm{brems} = \frac{E_\mathrm{e}}{4\alpha r_{0}^{2}cn_{\mathrm{p}}\beta _{\mathrm{e}}\gamma _{\mathrm{e}}m_\mathrm{e}c^2} \left[\ln (2\gamma _{\mathrm{e}})-\frac{1}{3}\right] ^{-1} \propto \frac{1}{n_{\mathrm{p}} \ln E_\mathrm{e}}, \end{eqnarray*}$$

(14)

where α is the fine-structure constant and r₀ denotes the electron radius. Again, adapting electrons with energies of 5 GeV, bremsstrahlung losses act in a medium with |$n_{\mathrm{p}}=0.1\, \mathrm{cm^{-3}}$| on time-scales of τ_brems ≈ 480 Myr. Consequently, high densities are required in order for bremsstrahlung losses to be able to compete with synchrotron or IC losses, as we will further discuss in Section 5.2.

3.3 Large dynamic range

The SFR-gas surface density relation, |$\dot{\Sigma }_\star \propto \Sigma _\rm {gas}^{1.4}$|⁠, of star-forming and starburst galaxies suggested by Kennicutt (1998) is valid over a large dynamic range of five orders of magnitude in gas surface density. As a result, the CR cooling time-scales that depend on gas density are expected to vary on a similarly large range of scales. Because the photon energy density is related to star formation, IC cooling strongly depends on the gas density via the Kennicutt (1998) relation, too. Similarly, the magnetic field strengths of star-forming galaxies are found to scale with gas surface density (Robishaw, Quataert & Heiles 2008), resulting in a strong variation of the synchrotron cooling time-scale with gas density as well. However, these cooling time-scales are not only expected to significantly vary among different types of galaxies, i.e. from dwarfs to starbursts, but also within a galaxy. In particular, the relative importance of the cooling and loss time-scales is of relevance for shaping their spectra and determining their non-thermal emission properties.

This is exemplified in Fig. 1, where we show maps of the time-scale of synchrotron and IC cooling of CR electrons and the hadronic time-scale of CR protons, both at an energy of 10 GeV (upper panels), for our fiducial halo at |$t=2.3\, \mathrm{Gyr}$|⁠. As expected, the synchrotron cooling time-scale traces the magnetic field (see Fig. 2) and varies from a few Myr in the very central regions with strong magnetic fields up to a few tens of Myr in the disc at around 5 kpc from the centre. At larger radii, the synchrotron time-scale increases with decreasing magnetic field as τ_syn ∝ B⁻². By contrast, the IC cooling scales as |$\tau _\mathrm{IC}\propto \varepsilon _\mathrm{ph}^{-1}$| and consequently, IC cooling remains important at larger radii, where the photon energy density is still high and approaches ε_CMB. Consequently, IC dominates over synchrotron cooling beyond 5 kpc (see lower left-hand panel in Fig. 1). Note that this ratio is independent of the electron energy, due to the identical energy dependence of IC and synchrotron losses. In those regions where synchrotron cooling dominates over IC cooling, the former is also faster than advection and diffusion losses, i.e. escape losses (lower middle panels in Fig. 1). However, in the outskirts of the galactic disc, IC losses dominate over escape losses and only within the outflows, advection losses become relevant.

Figure 2.

Face-on (upper panels) and edge-on (lower panels) maps of the CR electron distribution f_e at 10 GeV (left-hand panels), the magnetic field strength (middle panels) and the total (primary and secondary) synchrotron intensity at 1.4 GHz (right-hand panels) of a simulation with |$M_{200}=10^{12}\, \mathrm{M_{\odot }}$|⁠, |$B_0=10^{-10}\,$| G and ζ_SN = 0.05 at |$t=2.3\,$| Gyr (i.e. the same snapshot as shown in Fig. 1).

4 THE FIR–RADIO CORRELATION

We will now describe our modelling of the FRC from our simulations, which allows us to dissect the contribution from primary and secondary electrons to the total radio synchrotron luminosity. Furthermore, we will assess calorimetry and analyse possible deviations from the calorimetric picture in starburst galaxies.

4.1 Modelling the FRC

Post-processing our simulations, we derive the steady-state CR electron distributions f_e (equation 1) in every resolution element. We show the sum of the primary and secondary electron equilibrium distribution at 10 GeV in the left-hand panels of Fig. 2 for our fiducial halo, i.e. a simulation with halo mass |$M_{200}=10^{12}\, \mathrm{M_\odot }$|⁠, after 2.3 Gyr. Together with the magnetic field of our MHD simulation (shown in the middle panels of Fig. 2), we can directly calculate the resulting radio synchrotron emission of the steady-state electron distribution, using equation (7). This is shown in the right-hand panels of Fig. 2 in terms of the specific radio synchrotron intensity I_ν at a frequency of ν = 1.4 GHz as it would be observed in a face-on (upper panel) or edge-on (lower panel) configuration. The radio emission in our projected maps is clearly dominated by the central region up to radii of ∼10 kpc, where the magnetic field is strong. Particular filamentary features are visible in the edge-on view of the radio emission. Those trace the morphology of the magnetic field, which reaches up more than 10 kpc above and below the disc, and which is shaped by a strong central outflow driven by the CR pressure gradient.

Because the radio synchrotron emissivity j_ν in equation (7) depends on the component of the magnetic field perpendicular to the line of sight and hence on the viewing angle, the calculation of the observed radio luminosity also depends on the observed orientation. For the FRC, we chose to calculate the radio luminosity for a face-on configuration of our galaxies. The corresponding edge-on luminosities are typically a factor of about two smaller, which introduces a natural scatter in the FRC (Pfrommer et al. 2021). We correlate the face-on radio synchrotron luminosities against the SFR (⁠|$\dot{M}_{\star }$|⁠) of our simulated galaxies with an initial magnetic field of B₀ = 10⁻¹⁰ G and ζ_SN = 0.05 in Fig. 3. The corresponding FIR luminosities as derived from the SFRs using the Kennicutt (1998) relation are shown at the upper horizontal axis. For each simulated galaxy, we take snapshots at the peak of the SFR and at times when the SFR has successively decreased by an e-folding, which yields snapshots that are equally spaced in |$\log \dot{M}_{\star }$|⁠. The effect of varying the initial magnetic field B₀ and the injection efficiency of CRs at SNRs, ζ_SN, on the FRC is discussed in Appendix  B. Note that we only show the radio synchrotron luminosities here and discuss possible contributions from thermal free–free emission in Section 5.1.

Figure 3.

Top panels: we compare the observed FRC (Bell 2003, open circles and a fit to the data, grey line) to our simulated FRC (black symbols for the total radio synchrotron luminosity) in our model that accounts for CR advection and diffusion (‘CR diff’, left-hand panels) as well as our CR advection-only model (‘CR adv’, right-hand panels) and assume face-on configurations for all galaxies. Colour coded are contributions from primary (blue and orange) and secondary (green and red) CR electrons and different symbols correspond to simulations with different halo masses (indicated in the legends). All simulations assume |$B_0 = 10^{-10}\, \mathrm{G}$| and ζ_SN = 0.05. Bottom panels: we show the corresponding relative contributions to the total radio luminosity due to primary and secondary electrons, respectively. Note that the secondary emission dominates in our ‘CR adv’ model while the primary emission dominates in our ‘CR diff’ mode l.

Our simulation models which only allow for CR advection (‘CR adv’, right-hand panel of Fig. 3) and where we additionally account for anisotropic CR diffusion (‘CR diff’, left-hand panel) match the observed FRC (Bell 2003) over a broad range of SFRs. However, the radio luminosities of the smallest haloes with |$M_{200}=10^{10}\, \mathrm{M_\odot }$| fall short of the FRC in our ‘CR diff’ model and are not visible within the range of radio luminosities shown in Fig. 3. As discussed in Pfrommer et al. (2021), these dwarf galaxies show a slower dynamo growth in comparison to the more massive haloes. Additionally, they generate strong outflows that are launched due to the inclusion of anisotropic diffusion of CRs, enabling CRs to escape from the galactic disc, which further reduces the resulting radio emission. By contrast, in the ‘CR adv’ model, the small-scale dynamo efficiently amplifies the magnetic field and CRs accumulate in the galaxy because they cannot diffuse out of the disc by construction. Hence, our dwarf galaxies with |$M_{200}=10^{10}\, \mathrm{M_\odot }$| are able to reach the FRC in the ‘CR adv’ model (see Fig. 3, right-hand panel).

In addition to the total radio luminosity, Fig. 3 also shows the contributions arising from primary and secondary electrons separately, which are further analysed in the following section.

4.2 Secondary versus primary synchrotron emission

Complementary to Pfrommer et al. (2021), here we aim to understand the individual contributions of primary and secondary electrons to the total radio synchrotron luminosity, which enables us to dissect the reasons for our successful reproduction of the observed FRC. First of all, in our approach, the injected CR proton luminosity is given by

$$\begin{eqnarray*} L_{\mathrm{p}}=\zeta _{\mathrm{SN}} \dot{M}_{\star } \epsilon _{\mathrm{SN}} \ , \end{eqnarray*}$$

(15)

where |$\epsilon _{\mathrm{SN}}=E_{\mathrm{SN}}/M_{\star } = 10^{51}\, \mathrm{erg}/ (100~\mathrm{M_\odot })=10^{49}\, \mathrm{erg\, M_{\odot }^{-1}}$| is the SN energy released per unit mass (see Pfrommer et al. 2017a), where we assume a Chabrier (2003) initial mass function. Adopting our fiducial injection efficiency of CR energy at SN remnants, ζ_SN = 0.05 (Pais et al. 2018), and using the Kennicutt (1998) relation between SFR and FIR luminosity, we obtain

$$\begin{eqnarray*} L_{\mathrm{p}}\approx 5.5\times 10^{-4} \, \left(\frac{\zeta _{\mathrm{SN}}}{0.05}\right)\, L_{\mathrm{FIR}}. \end{eqnarray*}$$

(16)

The primary electron population obtains a fraction ζ_prim of the total luminosity of CR protons in our modelling that is given by (see Appendix A in Paper I)

$$\begin{eqnarray*} L_{\mathrm{prim-e}}=\zeta _{\mathrm{prim}}\, L_{\mathrm{p}}=\left(\frac{m_{\mathrm{p}}}{m_{\mathrm{e}}} \right)^{\alpha _{\mathrm{p}} -2}K_{\mathrm{ep}}^{\mathrm{inj}}\, L_{\mathrm{p}}. \end{eqnarray*}$$

(17)

Using |$K_{\mathrm{ep}}^{\mathrm{inj}}=0.02$| and α_p = 2.2, this yields ζ_prim ≈ 0.09.

On the other hand, hadronic interactions of CR protons with the ISM lead to the production of secondary electrons and positrons. Assuming that CR protons exclusively cool via hadronic interactions, the fraction of total proton luminosity injected into secondary electrons and positrons, i.e. ζ_sec = L_sec-e/L_p can be estimated from the following consideration. Charged pions are produced in hadronic CR proton–proton interactions with a multiplicity 2/3 and on average, electrons/positrons obtain 1/4 of the pion energy, which yields a secondary fraction of about 2/3 × 1/4 = 1/6 of the total proton luminosity. Considering additionally the nuclear enhancement factor of 1.4 to 1.6, which accounts for heavier nuclei in the composition of CRs and the ISM (Białłas, Bleszyński & Czyż 1976; Stephens & Badhwar 1981), we arrive at an energy fraction of ζ_sec ≈ 0.25. Hence, secondary electrons and positrons obtain a luminosity of

$$\begin{eqnarray*} L_{\mathrm{sec-e}}&=\eta _{\mathrm{cal,p}}\, \zeta _{\mathrm{sec}}\, L_{\mathrm{p}}, \end{eqnarray*}$$

(18)

where we additionally accounted for the calorimetric fraction of protons η_{cal, p}. This factor quantifies the fraction of CR proton luminosity that is used for pion production, which reaches values up to η_{cal, p} ≈ 0.8 for highly star-forming galaxies (see Paper II). We assume a simple power-law momentum spectrum for the source function of CR electrons with a low-momentum cutoff p_min and calculate the volume-integrated injected electron source function via:

$$\begin{eqnarray*} Q_\mathrm{e} (p_\mathrm{e})= \frac{\mathrm{d} N_\mathrm{e}}{\mathrm{d} p_\mathrm{e}\mathrm{d} t} = \int q_\mathrm{e} (p_\mathrm{e})\mathrm{d} V = \mathcal {C}_\mathrm{e} p_\mathrm{e}^{-\alpha _{\mathrm{e,inj}}}\theta (p_\mathrm{e}-p_\mathrm{min}), \end{eqnarray*}$$

(19)

where θ(p) denotes the Heaviside step function and |$\mathcal {C}_{\mathrm{e}}$| is the normalization. Assuming furthermore α_{e, inj} > 2, this enables us to define the total injected CR electron luminosity via

$$\begin{eqnarray*} L_\mathrm{e}&=&\int _0^{\infty }Q_\mathrm{e}(p_\mathrm{e})\, T_\mathrm{e}({\rm p}_\mathrm{e})dp_\mathrm{e}=\frac{\mathcal {C}_\mathrm{e}\, m_{\mathrm{e}}c^2}{\alpha _\mathrm{e,inj}-1}\nonumber\\ &&\times \left[\frac{1}{2}\, \mathcal {B}_\frac{1}{1+p_\mathrm{min}^2}\left(\frac{\alpha _\mathrm{e,inj}-2}{2},\frac{3-\alpha _\mathrm{e,inj}}{2}\right)+ p_\mathrm{min}^{1-\alpha _\mathrm{e,inj}}\, \frac{T_\mathrm{e}(p_\mathrm{min})}{m_{\mathrm{e}}c^2}\right],\nonumber\\ \end{eqnarray*}$$

(20)

$$\begin{eqnarray*} \equiv \mathcal {C}_\mathrm{e}\, m_{\mathrm{e}}c^2 A_\mathrm{bol} (p_\mathrm{min},\alpha _\mathrm{e,inj}), \end{eqnarray*}$$

(21)

where |$T_\mathrm{e}(p_\mathrm{e})=\left(\sqrt{1+p_\mathrm{e}^2}-1\right)\, m_{\mathrm{e}}c^2$| is the kinetic electron energy and |$\mathcal {B}_y(a,b)$| denotes the incomplete beta function (Abramowitz & Stegun 1965). Assuming that the synchrotron cooling time of CR electrons is shorter than their escape time (Voelk 1989), the emitted synchrotron luminosity from a steady-state electron distribution f_e can be calculated as

$$\begin{eqnarray*} \nu L_{\nu } (\mathrm{GHz}) &=& \displaystyle\frac{E_\gamma \mathrm{d} N_\gamma }{\mathrm{d}\ln \nu \, \mathrm{d} t} =\eta _{\mathrm{syn}}\, \displaystyle\frac{E_\mathrm{e}\mathrm{d} N_\mathrm{e}}{2\, \mathrm{d}\ln \gamma _\mathrm{e}\, \mathrm{d} t} \approx \eta _{\mathrm{syn}} \, \displaystyle\frac{\gamma _\mathrm{e}^{2-\alpha _\mathrm{e,inj}}}{2A_\mathrm{bol}}\, L_\mathrm{e} \end{eqnarray*}$$

(22)

$$\begin{eqnarray*} \equiv \eta _{\mathrm{syn}} \, \zeta _\mathrm{bol}(p_\mathrm{min},\alpha _\mathrm{e,inj},\gamma _\mathrm{e}) \, L_\mathrm{e}. \end{eqnarray*}$$

(23)

Here, we use dln ν = 2dln γ_e (see equation 8) and equation (21), while assuming the relativistic limit dp_e ≈ dγ_e in equation (19). Furthermore, we define the bolometric electron fraction |$\zeta _\mathrm{bol} = \gamma _\mathrm{e}^{2-\alpha _\mathrm{e,inj}} / (2A_\mathrm{bol})$|⁠, that accounts for the fraction of total electron luminosity that could potentially be converted into a specific synchrotron luminosity at frequency ν, given a magnetic field B and electrons with Lorentz factor γ_e(ν, B) (equation 8). In addition, we introduce the calorimetric synchrotron fraction η_syn, that quantifies the fraction of the available CR electron luminosity ζ_bolL_e that is actually converted into synchrotron emission, which will be discussed in more detail in Section 4.3.

Equations (19)–(23) can be separably adapted to primary and secondary electrons. Hence, we arrive at an FRC that reads for the specific luminosity as a function of FIR luminosity for each CR electron population as

$$\begin{eqnarray*} \nu L_{\nu ,\mathrm{prim}}\mathrm{(GHz)}\approx 1.0 \times 10^{-7}\, \left(\frac{\eta _{\mathrm{syn,prim}}}{0.1}\right)\left(\frac{\zeta _\mathrm{bol}}{0.02}\right) L_{\mathrm{FIR}}, \end{eqnarray*}$$

(24)

and

$$\begin{eqnarray*} \nu L_{\nu ,\mathrm{sec}}\mathrm{(GHz)}\approx 1.3\times 10^{-8}\, \left(\frac{\eta _{\mathrm{cal,p}}}{0.8}\right) \left(\frac{\eta _{\mathrm{syn,sec}}}{0.1} \right) \left(\frac{\zeta _\mathrm{bol}}{0.001}\right) L_{\mathrm{FIR}}, \nonumber \\ \end{eqnarray*}$$

(25)

where we use ζ_SN = 0.05. The 1.4 GHz radio luminosity as a function of SFR reads (using again Kennicutt 1998),

$$\begin{eqnarray*} L_{1.4\mathrm{GHz,prim}}\approx 2.0\times 10^{27}\, \left(\frac{\eta _{\mathrm{syn,prim}}}{0.1}\right) \left(\frac{\zeta _\mathrm{bol}}{0.02}\right) \left(\frac{\dot{M}_{\star }}{\mathrm{M_{\odot }\, yr^{-1}}}\right), \end{eqnarray*}$$

(26)

and

$$\begin{eqnarray*} L_{1.4\mathrm{GHz,sec}}&\approx& 2.6\times 10^{26}\, \left(\frac{\eta _{\mathrm{cal,p}}}{0.8}\right) \left(\frac{\eta _{\mathrm{syn,sec}}}{0.1}\right) \nonumber\\ &&\times \, \left(\frac{\zeta _\mathrm{bol}}{0.001} \right) \left(\frac{\dot{M}_{\star }}{\mathrm{M_{\odot }\, yr^{-1}}}\right), \end{eqnarray*}$$

(27)

where we adopted γ_e = 6 × 10³, which is the typical Lorentz factor of CR electrons emitting synchrotron radiation at 1.4 GHz in magnetic fields of |$2\, \mu \mathrm{G}$| (see equation 8) and ζ_bol(p_min = 1, α_{e, inj} = 2.2, γ_e = 6 × 10³) ≈ 0.02 for primary electrons, as well as ζ_bol(p_min = 1, α_{e, inj} = 2.7, γ_e = 6 × 10³) ≈ 0.001 secondary electrons. This analytical estimate yields a synchrotron luminosity from secondary electrons that is a factor of ∼10 smaller than the luminosity arising from primary electrons, if the assumptions η_{cal, p} = 0.8, η_{syn, prim} = η_{syn, sec} = 0.1 and ζ_bol = 0.02 hold for primary electrons, whereas ζ_bol = 0.001 for secondary electrons. We will discuss the choice of these parameters in the following.

The primary and secondary contributions to the total synchrotron luminosity |$L_{1.4\, \mathrm{GHz}}$| in our simulations (Fig. 3, dark-grey symbols) are separately shown in the upper panels of Fig. 3 with different colours, as indicated in the legend. The lower panels show the fraction of the primary and secondary luminosities, respectively. Whereas the primary luminosities fall short of the secondary emission in the ‘CR adv’ model (right-hand panels), the radio luminosity is strongly dominated by primary electrons in our ‘CR diff’ model (left-hand panels), with an increasing contribution from secondaries towards higher SFRs. This is due to the difference in both models in the parameters entering equation (27).

First, the calorimetric proton fraction η_{cal, p}, representing the efficiency of secondary electron production due to hadronic interactions, has been shown in Paper II to decrease with decreasing SFR in the ‘CR diff’ model from 0.7 in starburst galaxies, down to 0.1 in our dwarf galaxies. In the ‘CR adv model’, η_{cal, p} reaches values up to 0.8 and only typically varies by a factor of ∼2. The reduced efficiency of pion production due to CR diffusion can be understood by the fact that CR diffusion is another CR proton loss process that competes with hadronic losses and thus lowers the secondary radio luminosity. However, the difference of η_{cal, p} for star-forming galaxies with |$\dot{M}_{\star } \gtrsim 1~\mathrm{M_\odot ~yr^{-1}}$| is not significant enough to fully explain the discrepancy in L_prim versus L_sec in our different CR transport models.

A more relevant effect for explaining the sub-dominant role of secondary electrons for the total radio luminosity in the ‘CR diff’ model is the effect of the steepening of the CR proton spectra due to energy-dependent diffusion. As an example, we show in Fig. 4 CR spectra of CR protons, primary and secondary electrons of a simulation with |$M_{200}=10^{12}\, \mathrm{M_\odot }$| for both CR transport models at t = 1.3 Gyr. The averaged CR spectra (shown with colours) are compared to momentum power-law spectra that are typical for the various transport and cooling processes (shown with grey-dashed lines), in which we account for the full proton/electron dispersion relation. This causes the mild downturn of the proton spectrum at energies |$E_\rm {kin}\lesssim 1$| GeV. We can clearly see the effect of Coulomb-cooling, which suppresses the spectra in comparison to the momentum power-law spectrum and which is stronger for CR protons in comparison to primary electrons, due to the dependence of Coulomb-cooling on |$\beta =v /c$| (see equation 11 and Paper I).

Figure 4.

From the left- to right-hand panels, we show CR proton, primary, and secondary electron spectra of simulations with |$M_{200}=10^{12}\, \mathrm{M_{\odot }}$|⁠, ζ_SN = 0.05, |$B_0=10^{-10}\, \mathrm{G}$|⁠. We contrast the CR spectra of snapshots at |$t=1.3\, \mathrm{Gyr}$| of a simulation that includes advection and anisotropic diffusion of CRs (‘CR diff’, blue) to a simulation that only accounts for CR advection (‘CR adv’, yellow). In both cases, we average the spectra over the radial range, that includes 99 per cent of the total radio luminosity at 1.4 GHz (i.e. |$R\approx 7\, \mathrm{kpc}$| for both snapshots) and the height, where the magnetic field has decreased by an e-folding, as indicated in the legend. The grey dashed lines show momentum power-law spectra (with arbitrary normalization to exemplify scaling properties) that are typical for the various transport and cooling processes (with indices indicated in the panels). These spectra account for the full proton/electron dispersion relation, which causes the mild downturn of the proton spectrum at energies |$E_\rm {kin}\lesssim 1$| GeV.

At high energies, the CR proton spectrum stays flat in the ‘CR adv’ model, i.e. α_{p, adv} = α_inj = 2.2, whereas the spectral index approaches α_{p, diff} = 2.7 in the ‘CR diff’ model, which accounts for energy-dependent diffusion, i.e. D ∝ E^0.5. Because the steady-state proton distribution determines the source function of secondary electrons, α_{sec-e, inj} = α_p, the latter also exhibits a steeper steady-state spectrum in the model accounting for energy-dependent diffusion losses after accounting for the radiative steepening of the electron spectra due to IC and synchrotron interactions, α_sec-e = α_{p, diff} + 1 = 3.7. By contrast, the steady-state spectrum of primary electrons has a similar spectral shape in both models, with α_prim-e = 3.2 at high energies, where IC and synchrotron cooling dominate. Note that the spectra in Fig. 4 are averaged over the radii that include 99 per cent of the total radio luminosity in both snapshots, respectively, and a characteristic scale height of the magnetic field (as indicated in the legend) so that we show representative CR spectra for explaining the radio synchrotron emission at 1.4 GHz.

This effect of steeper secondary CR electron spectra is entering equation (27) in the form of ζ_bol, that is a function of p_min, γ_e(ν, B) and α_e (see equation 23). We exemplify these dependencies in Fig. 5, that shows the strong decrease of ζ_bol with increasing spectral index α_e of the CR electron source function, which affects the available luminosity of primary and secondary electrons, that can be converted to radio synchrotron emission. As expected, steeper electron spectra imply a smaller amount of electrons available for synchrotron emission at a fixed frequency. While the low-momentum cut-off of the CR electron spectra does not have a significant impact on ζ_bol, the electrons’ Lorentz factor γ_e naturally becomes increasingly more relevant for increasing spectral indices. As a result, ζ_bol is affected by the ν_c-effect.

Figure 5.

We show the bolometric electron fraction ζ_bol as implicitely defined in equation (23). It quantifies the fraction of the total electron luminosity (from a CR electron source function approximated by a power law with injection spectral index α_{e, inj} and low-momentum cutoff p_min, see equations 19 and 20) that is available for synchrotron emission at a frequency ν(γ_e, B). For instance, in an ambient magnetic field of |$B=2~\mu \mathrm{G}$|⁠, electrons with a Lorentz factor γ_e = 10⁴ radiate synchrotron emission at GHz frequencies (see equation 8).

Consequently, we identify two main effects that can potentially be responsible for suppressing the radio luminosity via ζ_bol: (i) the steeper secondary CR electron source functions with the limiting value of α_{sec-e, inj} = 2.7 due to CR diffusion decrease L_{1.4GHz, sec} (in comparison to primaries with α_{prim-e, inj} = 2.2 in both models and also secondaries with α_{sec-e, inj} = 2.2 in the ‘CR adv’ model), and (ii) the ν_c-effect, i.e. synchrotron emission at a fixed frequency in a lower magnetic field strength is generated by higher energetic electrons, which are less abundant. This implies a smaller value of ζ_bol, which leads to a lower observed radio luminosity. The second effect is subdominant for primary CR electrons in both CR transport models and for secondary CR electrons in the ‘CR adv’ model because ζ_bol does not significantly depend on γ_e for α_{e, inj} = 2.2. However, with increasingly larger α_{e, inj}, ζ_bol is reduced by up to a factor of 5, when we increase γ_e by an order of magnitude. This is in fact the case within our simulations: As pointed out by Pfrommer et al. (2021), our simulations exhibit saturated magnetic field strengths of 0.1 to |$14~\mu$|G for SFRs of 0.01 to 30 M_⊙ yr⁻¹ in our ‘CR diff’ model. This implies typical Lorentz factors of electrons emitting at 1.4 GHz that range from γ_e ≈ 3.5 × 10⁴ to ≈ 3.5 × 10³. As a result the ν_c-effect indeed plays an important role and is particularly strong for secondary electrons: as we move from starburst to dwarf galaxies the saturated magnetic field strengths decrease and imply increasing Lorentz factors of CR electrons, which diminishes the total radio luminosity (if it is observed at a fixed frequency). But in the case of our ‘CR diff’ model, diffusive losses become increasingly more important towards lower star-forming galaxies (Paper II), which steepens α_p = α_{inj, sec.e} and hence, the secondary contribution to the total radio luminosity experiences an additional suppression, which explains the decreasing fraction of secondary radio luminosity with decreasing SFR, as shown in the lower left-hand panel in Fig. 3 for our ‘CR diff’ simulations.

4.3 Testing electron calorimetry

In order to quantify electron calorimetry of our simulated galaxies as a function of SFR and halo mass, we first define the calorimetric synchrotron fraction following equation (23) as

$$\begin{eqnarray*} \eta _{\mathrm{syn,prim/sec}}=\frac{\nu L_{\nu ,\mathrm{prim/sec}}}{\zeta _\mathrm{bol} L_\mathrm{e}} =\frac{ \sum _i \nu L_{\nu ,\mathrm{prim/sec},i}}{\sum _i\zeta _{\mathrm{bol},i} L_{\mathrm{e},i}}, \end{eqnarray*}$$

(28)

where in each cell i, we define the fraction of the total CR electron luminosity that radiates synchrotron emission at a frequency ν, given a magnetic field B_i within a cell i,

$$\begin{eqnarray*} \zeta _{\mathrm{bol},i} = \frac{\left[\gamma _{\mathrm{e}}(B_i,\nu)\right]^{2-\alpha _\mathrm{e,inj}}}{2 A_\mathrm{bol}}, \end{eqnarray*}$$

(29)

and L_e = L_prim-e + L_sec-e denotes the electron luminosity. For simplicity, the calorimetric synchrotron fractions of primary and secondary electrons η_{syn, prim/sec}, are defined in such a way that they add up to a total calorimetric synchrotron fraction

$$\begin{eqnarray*} \eta _{\mathrm{syn}}=\eta _{\mathrm{syn,prim}} + \eta _{\mathrm{syn,sec}} = \frac{\nu L_\nu }{\zeta _\mathrm{bol} L_\mathrm{e}}. \end{eqnarray*}$$

(30)

To fulfill this condition of additivity, we assume α_{e, inj} = 2.2 in all cells for both primary and secondary electrons. This implies that η_{syn, sec} represents in our definition a lower limit for the calorimetric synchrotron fraction of secondary electrons, since they can exhibit steeper injected spectral indices α_{sec.e, inj} > 2.2, if energy-dependent CR diffusion is included, as discussed in Section 4.2.

To calculate these calorimetric synchrotron fractions of our simulations, we estimate the available proton luminosity from the SFR using equation (15) in each cell,² and use this to compute the electron luminosities from equation (17) and (18). In the upper panels of Fig. 6, we show the calorimetric synchrotron fractions of our simulated galaxies (the same snapshots as shown in Fig. 3) for primary and secondary electrons, as well as the total calorimetric fraction. These range from 0.03 to 0.65 in our ‘CR diff’ models and from 0.06 to 0.67 in our ‘CR adv’ models. This partly explains the large scatter that we obtain in our FRC in Fig. 3. In particular, there is a difference in the calorimetric synchrotron fraction for different halo masses at the same SFR, that leads to a spread in radio luminosities. This is a result of the growth of the magnetic dynamo that continues to amplify magnetic fields at large disc radii of the galaxies, after saturation of the small-scale magnetic dynamo (Pfrommer et al. 2021).

Figure 6.

Top panels: we show the calorimetric synchrotron fraction (see equation 28) for our CR diffusion (left-hand panel) and advection model (right-hand panel) as a function of SFR. In addition to the total fraction including primary and secondary electrons (black), we also show the fractions for primaries and secondaries separately (colour coded as indicated in the legend). Bottom panels: here, we show the bolometric fractions of synchrotron (filled symbols) and IC emission (open symbols) defined in equations (31) and (32) which add up to unity by construction. The different symbols correspond to simulations with different halo masses as indicated in the top left-hand panel.

Furthermore, we observe a global trend towards higher calorimetric fractions with higher SFRs, which is particularly strong in the ‘CR diff’ models. It remains to be seen whether an improved model for the ISM in global (dwarf) galaxy simulations with explicit SN-driven turbulence (Semenov, Kravtsov & Gnedin 2018; Gutcke et al. 2021) in a cosmological setting can further amplify the magnetic field in the outskirts of the disc via a small-scale dynamo. This would increase the values of |$\eta _\rm {syn}$| at low SFRs in dwarfs over what is found in our models, where fast CR diffusion in the ‘CR diff’ model in low-mass haloes efficiently quenches the small-scale dynamo and hence suppresses the calorimetric synchrotron fractions.

So far, we only considered the fraction of the injected electron luminosity that can potentially be converted to synchrotron emission at a fixed frequency, i.e. we only selected electrons with a certain Lorentz factor γ_e(ν, B) that emits into a given frequency window. But in addition to that, we are also interested in comparing the total (bolometric) amount of energy that is radiated via synchrotron emission in comparison to IC scattering. To extend these considerations over the whole energy range, we define bolometric fractions of synchrotron and IC emission according to

$$\begin{eqnarray*} \eta _{\mathrm{syn}}^{\mathrm{bol}} = \frac{\sum _i B_i^2V_i}{\sum _i \left(B_i^2+B_{\mathrm{ph},i}^2 \right)V_i} \end{eqnarray*}$$

(31)

and

$$\begin{eqnarray*} \eta _{\mathrm{IC}}^{\mathrm{bol}} = \frac{\sum _i B_{\mathrm{ph},i}^2V_i}{\sum _i \left(B_i^2+B_{\mathrm{ph},i}^2 \right)V_i}. \end{eqnarray*}$$

(32)

These definitions are derived from the synchrotron and IC loss rates, which depend on the energy density of the magnetic field ε_B ∝ B² and the photon energy density |$\varepsilon _\mathrm{ph}\propto B_\mathrm{ph}^2$|⁠, respectively (see equations 12 and 13). Per definition, |$\eta _{\mathrm{IC}}^{\mathrm{bol}}+\eta _{\mathrm{syn}}^{\mathrm{bol}}=1$|⁠, and hence, these fractions provide a measure of the fraction of the total electron luminosity that is lost to IC or synchrotron emission (at all possible frequencies) if we only consider these two radiative loss processes. To be consistent with the calculation of η_syn, we again only sum over all cells with SFRs |$\dot{M}_\star \gt 0$| and show the resulting bolometric fractions of IC and synchrotron emission in the lower panels of Fig. 6. In almost all analysed snapshots, we find that |$\eta _{\mathrm{IC}}^{\mathrm{bol}} \gt \eta _{\mathrm{syn}}^{\mathrm{bol}}$|⁠, which implies that IC losses are usually dominating synchrotron losses for CR electrons. For our small haloes with |$M_{200}=10^{10}$| and |$10^{11}\, \mathrm{M_\odot }$|⁠, this is not surprising: Pfrommer et al. (2021) found that the magnetic field in these small haloes saturate below the energy density of the CMB and hence, IC-cooling via scattering off of the CMB alone already dominates over synchrotron cooling. There are only three cases for which |$\eta _{\mathrm{IC}}^{\mathrm{bol}} \le \eta _{\mathrm{syn}}^{\mathrm{bol}}$|⁠, which represent highly star-forming galaxies with |$M_{200}=10^{12}\, \mathrm{M_\odot }$| in the ‘CR adv’ model. They saturate at the highest magnetic field strengths in comparison to all other shown simulations (Pfrommer et al. 2021) and hence, result in the highest bolometric fractions of synchrotron emission, that manage to dominate over IC emission.

4.4 ‘Conspiracy’ at high gas densities

We found in Section 4.2, that there is an increasing contribution of secondary emission to the total radio luminosity toward larger SFRs in our ‘CR diff’ model. If this trend holds, the question arises of how to maintain the almost linear behaviour of the FRC. As proposed by Lacki et al. (2010), in addition to the larger contribution of secondary electrons to the radio luminosity in high-density starburst galaxies, one would also expect more losses due to Coulomb interactions and bremsstrahlung emission due to their dependence on gas density, see equations (11) and (14). Furthermore, IC losses are expected to increase with higher SFRs due to the higher photon energy densities from young stellar populations. Consequently, in order to maintain an almost linear FRC, those effects have to ‘conspire’ to (approximately) cancel each other (Lacki et al. 2010).

To quantify deviations from the FRC and to test the relative contributions of the proposed effects, we utilize the ratio of total IR (TIR) to radio synchrotron emission as defined by Helou et al. (1985):

$$\begin{eqnarray*} q=\log _{10}\left(\frac{L_{\mathrm{TIR}}}{L_{1.4\, \mathrm{GHz}}} \right) - 3.67, \end{eqnarray*}$$

(33)

where the TIR luminosity is given by L_TIR ≈ 1.75L_FIR (Calzetti et al. 2000). We define the deviation from a given FRC via the difference Δq that results from comparing the full and a modified FRC for which we exclude different contributions to the total synchrotron emission. This enables us to study the impact of individual processes in modifying the slope of the FRC and whether these processes pose a challenge to the observed quasi-linearity of the FRC.

First, we assess the effect of neglecting the contribution of secondary electrons to the total radio luminosity |$L_{1.4\, \mathrm{GHz}}$|⁠, and denote the corresponding parameter with q_no sec. This is expected to be larger than q, where primary and secondary contributions are included. Hence, the difference

$$\begin{eqnarray*} \Delta q_{\mathrm{sec}}= q - q_{\mathrm{no\ sec}}\le 0 \end{eqnarray*}$$

(34)

enables us to quantify the contribution of secondary radio emission to the FRC.

Second, we consider the case of disregarding losses due to bremsstrahlung and Coulomb interactions to infer their effect on the radio emission and define the resulting difference as

$$\begin{eqnarray*} \Delta q_{\mathrm{brems+Coul}} = q - q_{\mathrm{no\ brems+Coul}}. \end{eqnarray*}$$

(35)

To calculate q_{no brems + Coul}, we compute the radio luminosity in equation (33) from steady-state spectra that we obtain by setting b_Coul = b_brems = 0. If these losses are significant in shaping the steady-state distribution, we expect |$q_{\mathrm{no\, brems+Coul}}\le q$| because fewer losses due to processes other than synchrotron emission potentially give rise to a larger radio luminosity. Similarly, we calculate the change in q from neglecting all non-synchrotron cooling processes, i.e. bremsstrahlung, Coulomb and IC cooling, which is denoted by

$$\begin{eqnarray*} \Delta q_{\mathrm{brems+Coul+IC}} = q - q_{\mathrm{no\ brems+Coul+IC}}. \end{eqnarray*}$$

(36)

Fig. 7 shows the deviation from a linear FRC when excluding secondary radio emission or bremsstrahlung and Coulomb losses. In our ‘CR diff’ runs, both effects become relevant in galaxies with SFRs |$\dot{M}_\star \gtrsim 1~\mathrm{M_\odot \, yr^{-1}}$|⁠. In these star-forming and star-bursting galaxies, the missing radio emission from excluding secondary CR electrons is compensated by the additional radio emission, resulting from the neglect of bremsstrahlung and Coulomb losses. As a consequence, the deviation from calorimetry that would be expected in starbursts, due to the increasing relevance of bremsstrahlung and Coulomb interactions in these high-density galaxies, seems to be counterbalanced by including the contribution of secondary radio emission.

Figure 7.

We show the change in the FIR-to-radio luminosity ratio, quantified by q (see equation 33), for our ‘CR diff’ (left-hand panel) and ‘CR adv’ models (right-hand panel). Here, we study how individual cooling processes modify the FRC: the effect of neglecting bremsstrahlung and Coulomb losses, that both depend on gas density, is shown via Δq_{brems + Coul} (equation 35) in light blue (orange) with open symbols in our ‘CR diff’ (‘CR adv’) models, whereas the effect of additionally neglecting IC losses, Δq_{brems + Coul + IC}, is shown by the corresponding full symbols. Changes in q that result from disregarding the secondary radio emission are denoted by Δq_sec (equation 34) and are shown in dark blue (red). Different symbols correspond to different halo masses as indicated in the left-hand panel. Note that both panels have different ranges on their vertical axes.

IC losses elevate q across all SFRs and show an almost constant shift of q for our ‘CR diff model’, indicating that the IC emission does not significantly modify the slope of the FRC in this model. While IC losses become increasingly more relevant toward high SFRs for every individual galaxy simulation in our ‘CR adv’ model, on average the distribution of Δq_{brems + Coul + IC} is an elevated version of Δq_{brems + Coul} with an increased scatter. Hence, in spite of most of CR electrons lose their energy through IC interactions (Fig. 6), this shows that the subdominant synchrotron cooling process nevertheless appears to be a quasi-calorimetric tracer of the SFR, implying an almost linear slope of the FRC. Note that the overall change in q is comparably small in the ‘CR diff’ model and Δq ≤ 0.37 for all considered effects. For a mean observed q of 2.21 in starburst galaxies (Helou et al. 1985), this is at most a 17 per cent effect.

In contrast to our ‘CR diff’ model, galaxies of our ‘CR adv’ model (right-hand panel in Fig. 7) exhibit a contribution of secondaries that is almost independent of SFR. This is because neglecting diffusion losses increases the importance of hadronic losses, which are faster than escape losses |$\tau _{\pi}$| ≪ τ_esc. As a result, the ratio of primaries to secondaries does not depend on the gas density or the SFR, but only depends on |$K_\mathrm{ep}^{\mathrm{inj}}$| and α_p (see equation 6), and is hence approximately constant across the range of SFRs probed by our set of simulations that only accounts for CR advection. Furthermore, as discussed in Section 4.2, the ν_c-effect can only significantly affect secondary electrons in the model accounting for energy-dependent CR diffusion, which results in an SFR-dependence of the secondary contribution to the total radio synchrotron emission and hence explains the decrease of Δq_sec with increasing SFR in the ‘CR diff’ model.

5 RADIO SPECTRA

As pointed out by Thompson et al. (2006), the calorimetric assumption must hold for starburst galaxies like Arp 220. Here, the large photon energy density implies a short IC cooling time-scale:

$$\begin{eqnarray*} \tau _\mathrm{IC}\sim 4\times 10^3\, \mathrm{yr}\, \left(\frac{\varepsilon _\rm {ph}}{10^{-6}~\rm {erg~cm}^{-3}}\right)^{-1}\, \left(\frac{B}{3~\rm {mG}}\right)^{1/2}. \end{eqnarray*}$$

(37)

In order for escape losses to be faster, unphysical wind velocities of order |$\sim 20\,000\, \mathrm{km\, s^{-1}}$| would be required for CR electrons to be advected from the compact star-bursting region of size ∼100 pc. Consequently, losses due to IC emission must be larger than escape losses and we would expect steep radio spectra. This is because an IC-cooled steady-state CR electron spectrum with α_e = 3.2 implies a synchrotron spectral index of α_ν = (α_e − 1)/2 = 1.1, which is larger than the observed spectral indices of ∼0.5 to 0.8. The same argument holds, if synchrotron losses dominate the cooling of CR electrons, due to the identical energy dependence of synchrotron and IC losses.

In order to solve this tension, we analyse in the following three possible mechanisms that could be responsible for flattening the radio spectra, in comparison to the steep radio spectra that are expected if IC and synchrotron cooling dominate. First, a non-negligible contribution of thermal free–free emission to the radio emission could have an influence on the radio spectral shape. Second, because bremsstrahlung losses have a weaker energy dependence, i.e. b_brems ∝ Eln E (equation 14), in comparison to IC and synchrotron losses with b_IC/syn ∝ E², they could lead to electron spectra that reflect the injected spectral index and hence lead to flatter radio spectra, if they are faster than IC and synchrotron losses, i.e. |$\tau _{\mathrm{brems}}\lt (\tau _{\mathrm{IC}}^{-1}+\tau _{\mathrm{syn}}^{-1})^{-1}$|⁠. Similarly, Coulomb cooling could affect the radio spectral shape because b_Coul ≈ const. for γ_e ≫ 1. A third possibility for maintaining flat electron spectra and hence yielding flat emitted synchrotron spectra is advection. Because advection losses do not depend on energy, the injected electron spectral index of 2.2 would be maintained in the regions where those losses are dominant. We test these possibilities in the following sections.

So far, in our fiducial model, we adopted an energy dependence of the diffusion coefficient D ∝ E^δ where δ = 0.5. In order to be consistent with our findings in Paper II, where we found that a shallower dependence of δ = 0.3 enables better matches to the observed gamma-ray spectra of NGC 253, M82, and NGC 2146, we adopt in the following δ = 0.3 for the calculation of the radio spectra.³

5.1 Thermal and non-thermal radio emissions

Figs 8 and 9 show the radio spectra of the three star-forming galaxies NGC 253, M82, and NGC 2146. We compare observations by Kapińska et al. (2017), Adebahr et al. (2013) and Klein et al. (2018) to radio spectra derived from our simulated galaxies that resemble the observed ones in terms of their SFRs and radio luminosities, respectively (see Table 2). We show the observed and simulated total spectra (top panels) and the spectra from the central regions (middle panels). The flux density is inferred from our modelling of the intensity I_ν, as it would be observed from a galaxy with an inclination angle ϕ at a luminosity distance d (see equations A13 and A14). Below the total and central radio spectra, we show the spectral indices of the differently modelled components that are computed via

$$\begin{eqnarray*} \alpha _\nu = -\frac{\mathrm{d}\log F_\nu }{\mathrm{d}\log \nu }. \end{eqnarray*}$$

(38)

We find that the total radio synchrotron spectra exhibit spectral indices from 0.6 to 1.1 for radio frequencies ranging from 0.1 to 100 GHz. Those are significantly steeper than the observed radio spectra at frequencies larger than a few GHz. However, our modelling of the thermal free–free emission (see Section 2.4.2) is able to flatten the total spectra so that they approach the observed values also at frequencies ν ≳ 1 GHz. This has also been found for the starburst galaxies NGC 253, M82, and Arp220 in Peretti et al. (2019), although in contrast to our results, they find that the radio emission is dominated by secondary emission.

Figure 8.

Spectra of synchrotron emission (blue lines) and thermal free–free emission (magenta lines) of our simulations that are similar to NGC 253 and M82 in terms of SFR and total radio luminosity (see Table 2). The spectra are computed for inclination angles ϕ = 74° for NGC 253 (Iodice et al. 2014) and ϕ = 80° for M82 (Lynds & Sandage 1963; McKeith et al. 1995). The upper panels show the spectra for the whole galaxies, whereas for the middle panels, we cut out a region around the galactic centre with a radius R_central (see Table 2; dashed lines), resembling the observed central regions, respectively. To allow a comparison with the spectral shapes of the observed spectra, the simulated spectra are re-normalized by L_{1.4GHz, obs}/L_{1.4GHz, sim} (see Table 2). Thermal free–free absorption affects the synchrotron spectra at low frequencies (green lines) and is strongest for the central spectra, due to the high central gas density. For NGC 253 and M82, we plot the total (black points) and central (grey points) observed radio spectrum by Kapińska et al. (2017), Klein, Lisenfeld & Verley (2018) and Adebahr et al. (2013), as indicated in the legend. The light grey points in the middle right-hand panel correspond to the spectrum of the halo, i.e. the difference between the total and central spectrum, which matches the spectral shape of our synchrotron component. Bottom panels: The solid (dashed) lines show the spectral indices of our total (central) simulated radio spectra for synchrotron, thermal free–free and total emission, respectively. For a variation of R_central, see Appendix A3.

Figure 9.

Top panel: we compare the observed total radio spectrum of NGC 2164 (Klein et al. 2018) to simulated spectra of synchrotron and thermal free–free emission of a simulation that matches NGC 2164 in terms of SFR and total radio luminosity (see Table 2). Here, we adopt an inclination angle of ϕ = 63° (Della Ceca et al. 1999). The middle panel shows the spectra in the central region (within |$R_\mathrm{central}=0.3\, \mathrm{kpc}$|⁠) while the bottom panel shows the spectral index of our simulated radio spectra for synchrotron, thermal free–free and total emission, respectively.

Table 2.

This table summarizes the properties of our simulated galaxies, that resemble the observed galaxies NGC 253, M82 and NGC 2146 in terms of their SFRs and total radio luminosities. The observed radio luminosity is derived from the observed flux density at the frequency bin that is closest to |$\nu =1.4\, \mathrm{GHz}$| using the data shown in Figs 8 and 9 and from the distances summarized below. We use the simulated radio luminosities with the observationally determined inclination angle ϕ and adapt δ = 0.3. The central radius R_central and the constant fraction ξ_e of the electron density provided in our pressurized ISM model (Springel & Hernquist 2003) are parameters used to construct the radio emission spectra from the simulations, including thermal free–free emission and absorption.

Galaxy	SFR (observationa	Distanceb	|$L_{1.4\, \mathrm{GHz}}$| (observation/	Rcentral(observation/	ϕ(observationc /	ξe	Simulation
	/simulation)		simulation)
	( M⊙ yr−1)	(Mpc)	(erg s Hz−1)	simulation) (kpc)	simulation)		|$M_{200}\, (\mathrm{M_{\odot }}),~t\, (\mathrm{Gyr})$|
NGC 253	5.03 ± 0.76	3.3	7.67 × 1028	0.15–0.25	74°	–	–
	4.110	–	3.87 × 1028	0.15	74°	0.8	M200 = 3 × 1011, t = 1.1
M82d	10.4 ± 1.6	3.7	1.21 × 1029	0.45	80°	–	–
	6.457	–	1.62 × 1029	1.50	80°	1.0	M200 = 1012, t = 2.3
NGC 2146	14.0 ± 0.5	15.2	3.02 × 1029	–	63°	–	–
	25.520	–	5.55 × 1029	0.30	63°	0.5	M200 = 1012, t = 0.7

a

Kornecki et al. (2020).

b

NGC 253: Mouhcine et al. (2005); M82: Vacca et al. (2015); NGC 2146: Gao & Solomon (2004).

c

NGC 253: Iodice et al. (2014); M82: Lynds & Sandage (1963), McKeith et al. (1995); NGC 2146: Della Ceca et al. (1999).

d

This snapshot is shown in Figs 1, 2, 10, 11, and 13.

The spectrum of the central region of a simulated galaxy is obtained by computing the intensity |$I_\nu (\boldsymbol {r}_\perp)$| within a circle around the centre with the radius R_central (see Table 2) that approximately corresponds to the radius of the observed central region. Only for M82, we adapt a larger radius to match the total observed flux from the central region with our model. Interestingly, we find that the central data both of NGC 253 and M82 show an even flatter radio spectrum in comparison to the total spectrum, with particularly strong free–free absorption at small frequencies. This is convincingly explained by the higher electron density in central regions of galaxies, as it is also manifested within our model of NGC 253. It also becomes clear from the observations of M82, that the flattening of the radio spectra is a feature that is predominantly originating from the central region of the galaxy. The observations by Adebahr et al. (2013) of the central ∼450 pc of M82 clearly show a flat radio spectrum, whereas the halo spectrum (i.e. the total spectrum minus the central data) is similarly steep in comparison to our radio synchrotron component. Hence, the flat component of the radio spectrum must result from the central region and can be reconciled with the steep radio synchrotron spectra.

In Fig. 10, we show modelled free–free emission maps at 1 GHz of our M82-like galaxy to illustrate the dependence of the thermal emission on the observed region and also on the viewing angle. Observing a galaxy edge-on leads to a higher column density along the line of sight, which enhances the thermal free–free emission in the mid-plane. Furthermore, the ratios of radio synchrotron to thermal free–free emission at two radio frequencies shown in the middle and right-hand panels of Fig. 10 highlight our conclusion drawn from the radio spectra as discussed above: while the synchrotron emission dominates at low frequencies, thermal free–free emission prevails towards higher frequencies.

Figure 10.

Face-on and edge-on maps of the projected thermal free–free emission at 1 GHz (left-hand panels) of our simulation with |$M_{200}=10^{12}\, \mathrm{M_\odot }$| at t = 2.3 Gyr. The middle and right-hand panels show the ratio of the synchrotron emission to the thermal free–free emission at 1 and 30 GHz, respectively.

Similarly to NGC 253 and M82, the observed spectrum of NGC 2146 is flatter than the radio synchrotron emission of our simulated galaxy. However, adding a thermal free–free component enables us to reproduce the flat observed spectrum at frequencies above ∼10 GHz. In the case of NGC 2146, there are no published data of the central region. We chose to show the central spectrum of NGC 2146 from a region with radius of 0.3 kpc, where we find a pronounced flattening of the central spectrum at low frequencies due to free-fee absorption.

We furthermore note that the radio luminosities given in Table 2 for our simulations include radio synchrotron emission and thermal free–free emission, whereas the luminosities adapted in our analysis of the FRC are only the non-thermal radio synchrotron luminosities. We determine the contribution of the thermal free–free emission to the total luminosity at 1.4 GHz to 24.7, 6.3, and 5.7 per cent for NGC 253, M82, and NGC 2146, respectively. Except from NGC 253, this is in agreement with thermal fractions at 1.5 GHz derived from observations of 61 nearby galaxies in Tabatabaei et al. (2017).

5.2 Can bremsstrahlung or Coulomb losses yield flat radio spectra?

In order to assess the possibility that bremsstrahlung cooling can lead to flatter radio spectra in comparison to the case when IC and synchrotron losses dominate, we first show in Fig. 11 a slice through the disc of the bremsstrahlung cooling time-scale at an energy of 10 GeV, which is typically acting on time-scales of 50−200 Myr and is shortest in high-density regions (see equation 14). However, synchrotron cooling is faster in the central few kpc and in the outflow, where the magnetic field is strong. In the regions where bremsstrahlung losses are faster than synchrotron cooling, IC cooling kicks in (right-hand panel of Fig. 11). Therefore, bremsstrahlung losses are subdominant in most regions of this galaxy and are not able to significantly shape the electron spectra at the considered electron energy. Consequently, they cannot be responsible for flattening these at the corresponding frequencies. Instead, IC and synchrotron losses steepen the injected electron spectra by unity, yielding α_ν = 1.1 in the fully cooled limit.

Figure 11.

Face-on and edge-on maps of the time-scale of bremsstrahlung cooling, |$\tau _\rm {brems}$| (left-hand panels), as well as the ratio of |$\tau _\rm {brems}$| to the cooling time-scale due to synchrotron and IC emission (middle and right-hand panels, respectively) at an electron energy of 10 GeV. The maps are shown for our fiducial halo with |$M_{200}=10^{12}\, \mathrm{M_\odot }$| at 2.3 Gyr and the ratios are averaged over slices with a thickness of 0.5 kpc.

However, as can be seen from the radio spectral index shown in the bottom panels of Fig. 8 and 9, α_ν only approaches this limit of very high frequencies ν ≳ 10 GHz. One reason for this is the fact that we include diffusive losses in our steady state equation,⁴ which only steepen the spectra by δ = 0.3 (or 0.5). This implies that we obtain radio spectral indices ranging from 0.75 (or 0.85) if diffusion losses dominate, to 1.1 if IC or synchrotron losses prevail. Furthermore, in regions with particularly large magnetic field strengths, we observe lower energetic electrons (due to the ν_c-effect), which are more affected by Coulomb and bremsstrahlung cooling in comparison to high-energy electrons, where IC and synchrotron losses steepen the spectra.

To quantify this effect and to assess the importance of bremsstrahlung cooling in flattening radio spectra for galaxies along the SFR sequence, we show in Fig. 12 the ratio of the combined time-scale of bremsstrahlung and Coulomb cooling

$$\begin{eqnarray*} \tau _{\mathrm{brems,Coul}}= \left(\tau _{\mathrm{brems}}^{-1} + \tau _\mathrm{Coul}^{-1} \right)^{-1} \end{eqnarray*}$$

(39)

to the combined synchrotron and IC time-scale, i.e.

$$\begin{eqnarray*} \tau _{\mathrm{syn,IC}}= \left(\tau _{\mathrm{syn}}^{-1} + \tau _\mathrm{IC}^{-1} \right)^{-1}, \end{eqnarray*}$$

(40)

as well as to the combined time-scale of synchrotron, IC and escape losses,

$$\begin{eqnarray*} \tau _{\mathrm{syn,IC,esc}}= \left(\tau _{\mathrm{syn}}^{-1} + \tau _\mathrm{IC}^{-1} + \tau _\mathrm{esc}^{-1}\right)^{-1}. \end{eqnarray*}$$

(41)

In order to show the cooling time-scale ratios in regions that are relevant for synchrotron emission, we weight the corresponding cooling rates (τ⁻¹) with the synchrotron luminosity of each cell when averaging them over the whole galaxy. In both, our ‘CR diff’ and ‘CR adv’ models, we find cooling via bremsstrahlung and Coulomb losses to be slower than IC and synchrotron cooling for electron with energies E_e = γ_e(B)m_ec²|_10 GHz, i.e. electrons that are typically responsible for the emission of synchrotron radiation at 10 GHz, depending on the ambient magnetic field (see equation 8).⁵

Figure 12.

Ratios of the synchrotron-weighted time-scales of bremsstrahlung and Coulomb cooling τ_{brems, Coul} (equation 39) to the time-scale of synchrotron and IC cooling τ_{syn, IC} (as well as to that of radiative cooling and escape losses τ_{syn, IC, esc}, open symbols; see equations 40 and 41) for our ‘CR diff’ model (left-hand panels) and ‘CR adv’ model (right-hand panels). The time-scales of each snapshot are calculated by averaging over the corresponding cooling rates (i.e. τ⁻¹) from all cells while weighting them with the synchrotron luminosity. The time-scales are computed at the typical electron energy responsible for radio synchrotron emission at |$\nu =10\, \mathrm{GHz}$| (top panels) and |$\nu =1.4\, \mathrm{GHz}$| (lower panels), given the magnetic field in each cell (see equation 8).

However, if we consider synchrotron-emitting electrons at ν = 1.4 GHz we find that bremsstrahlung and Coulomb losses are not at all negligible in comparison to IC and synchrotron losses in most of our snapshots, but instead, τ_{brems, Coul} ≲ τ_{syn, IC} at electron energies E_e = γ_e(B)m_ec²|_1.4 GHz. This is due to the fact that we typically probe lower energetic electrons when observing radio synchrotron emission at a lower frequency, given the same magnetic field. The mean Lorentz factor of electrons emitting at 1.4 GHz averaged over all cells ranges between γ_e = 10⁴ and 2 × 10⁴ for all SFRs, corresponding to E_e = 5–10 GeV. By contrast, the average γ_e weighted with the synchrotron luminosity in each cell decreases from 10⁴ at small SFRs down to 10³ in our starburst galaxies. Hence, the synchrotron emission at ν = 1.4 GHz in starburst galaxies with their high magnetic field strengths probes electrons with energies down to E_e = 0.5 GeV.

The additional inclusion of escape losses (open symbols in Fig. 12) that become increasingly more important for decreasing SFRs in our ‘CR diff’ model shows that only highly star-forming galaxies are able to maintain τ_{brems, Coul} ≲ τ_{syn, IC, esc} at the relevant energies, determined by γ_e(B) (equation 8). This is in accordance with our results from Section 4.4, where we found that neglecting bremsstrahlung and Coulomb losses increases the radio luminosity at 1.4 GHz for an increasing SFR. In the ‘CR adv’ model, we observe a similar trend in comparison to the ‘CR diff’ model: At high electron energies (e.g. at 10 GeV), bremsstrahlung and Coulomb losses are subdominant while they become relevant if we account for the ν_c-effect, where we typically probe lower energetic electrons. If we additionally include escape losses (i.e. only advection losses in this model), they do not significantly influence the considered time-scale ratios.

In summary, we conclude that the synchrotron emission at 1.4 GHz in highly star-forming galaxies probes electrons with γ_e ∼ 10³ which are mostly affected by bremsstrahlung and Coulomb losses. By contrast, electrons at higher energies are dominated by synchrotron and IC cooling and consequently, at higher frequencies, the radio synchrotron spectra are steepened by those cooling processes. Only the modelling of a thermal free–free component as discussed in Section 5.1 is thus able to flatten the radio spectra towards high frequencies in starbursts. We furthermore find an increasing relevance of energy-dependent diffusion losses of CR electrons with decreasing SFRs. At the same time, observations at 1.4 GHz probe higher energetic electrons (γ_e ∼ 10⁴) at low SFRs so that bremsstrahlung and Coulomb losses do not strongly affect the radio synchrotron luminosity towards low SFRs (see Fig. 7).

5.3 Can outflows in projection yield flat radio spectra?

We show in Fig. 13 maps of the projected synchrotron emissivity |$I_{1.4\, \mathrm{GHz}}$|⁠, resulting from primary (left-hand panels) and secondary electrons (middle panels) for our fiducial halo with |$M_{200}=10^{12}\, \mathrm{M_\odot }$| at |$t=2.3\, \mathrm{Gyr}$|⁠. Because the simulation shown here includes anisotropic diffusion, the radio luminosity from primary electrons dominates over the secondary contribution (as discussed in Section 4.2). From the top to bottom panels, we show projected views of the simulated galaxy, changing the galaxy inclination from a face-on to an edge-on view in steps of 30°. The right-hand panels show spectral index maps of the total synchrotron emission between 1 and 3 GHz.

Figure 13.

From the top to bottom panels, we show projected maps of the radio synchrotron intensity (equation 9) arising from primary (left-hand panels) and secondary electrons (middle panels) observed from different inclination angles (0º, 30º, 60º, and 90º) for our fiducial galaxy with |$M_{200}=10^{12}\, \mathrm{M_\odot }$| at |$t=2.3\, \mathrm{Gyr}$|⁠. The right-hand hand panels show the corresponding maps of the spectral index of the total radio synchrotron emission (i.e. the primary plus secondary contribution) between frequencies of 1 and 3 GHz.

As discussed above, the spectral index of the steady-state electron population remains unchanged, if advection is the dominant loss process, i.e. α_e = α_inj = 2.2. Hence, the radio spectral index is α_ν = (α_e − 1)/2 = 0.6 in the outflows, where CRs predominantly are advected with the gas into the halo. By contrast, the observed spectral index maps of the radio emission of starburst galaxies like M82 (Adebahr et al. 2013), show flat spectral indices mainly in their central regions. We conclude that CR advection in combination with galaxy inclination cannot be responsible for the observed flat radio spectral indices because regions with flat spectral indices (due to advection losses) are seen off-centre and only for inclinations larger than about 45°.

6 DISCUSSION AND CONCLUSIONS

In this work, we model for the first time the radio emission of galaxies in three-dimensional MHD simulations, which evolve the energy density of CR protons self-consistently. To this end, we determine the steady-state spectra of CR protons, primary, and secondary electrons in post-processing and assess the relative importance of primaries versus secondaries for models with and without anisotropic CR diffusion. The detailed modelling of the synchrotron and free–free radio emission of star-forming galaxies enables new insight into the underlying physics of the FRC over a broad range of SFRs. In particular, this novel approach sheds light on the long-standing puzzle whether the almost linear FRC implies electron calorimetry and how this compares to the observed flat radio spectra.

While in simulations that only account for CR advection (‘CR adv’) the secondary electron population is responsible for most of the radio synchrotron emission and independent of SFR, we find in our model that additionally includes anisotropic CR diffusion (‘CR diff’) that primary electrons dominate the total radio luminosity at 1.4 GHz (see Fig. 3). The main reason for this difference is the steepening of the steady-state CR proton spectra due to energy-dependent diffusion, which modifies the shape of the secondary electron source function. Thus, in steady state, the secondary CR electron spectra differ from the primary electron spectra (see Fig. 4). As a result, due to their steeper spectra, secondary electrons are more affected by the ν_c-effect, which means that the observed radio synchrotron emission at a given frequency probes lower electron energies for larger magnetic field strengths that are realised in highly star-forming galaxies. This leads to an increasing contribution of secondary emission to the total radio luminosity with SFR in the ‘CR diff’ model, in contrast to the ‘CR adv’ model, where this trend is absent.

To test electron calorimetry, we determine the fraction of available injected CR electron luminosity that is actually converted into synchrotron emission. This calorimetric synchrotron fraction |$\eta _\rm {syn}$| varies from 0.03 to 0.67 among all our simulations with an increasing trend towards larger SFRs (see Fig. 6). The galaxies in the ‘CR adv’ model show overall larger calorimetric fractions with a weaker dependence on SFR in comparison to the ‘CR diff’ model, where diffusive losses and a weaker saturated magnetic field strength in dwarfs (in comparison to the thermal energy density, Pfrommer et al. 2021) decrease |$\eta _\rm {syn}$| towards smaller SFRs. We anticipate that an improved ISM model (e.g. Rathjen et al. 2021) in combination with an improved two-moment CR transport scheme (e.g. Thomas & Pfrommer 2019) within a cosmological setting that exhibits several epochs of accretion-driven star-forming phases can further amplify the magnetic field strength in these systems via star formation and accretion-driven turbulence as well as modulate the influence of CR transport on the magnetic dynamo amplification via a more consistent self-generated CR diffusion coefficient. This is expected to weaken the dependence of |$\eta _\rm {syn}$| on SFR. Generally, we find that CR electrons with energies ≳ 10 GeV lose most of their energy through IC interactions with starlight and CMB photons in comparison to synchrotron losses and thus are approximately in the calorimetric limit. However, we also find that the contribution to IC cooling is on average largely independent of SFR at fixed observed radio frequency and hence, enables to use the synchrotron emission as a quasi-calorimetric measure of the SFR (see Fig. 7).

Furthermore, we show that the increasing secondary contribution to the total radio luminosity in starbursts is balanced by increasing bremsstrahlung and Coulomb losses, that diminish the radio synchrotron luminosity with increasing densities in those systems (see Fig. 7). This finding is in accordance with one-zone models by Lacki et al. (2010). However in our model, the effect is not as strong: Lacki et al. (2010) find an increase in the logarithm of the TIR-to-radio luminosity ratio q of 0.5 to 0.6 if secondaries are not included in starbursts, as well as a corresponding decrease if secondaries are included but non-synchrotron losses are disregarded. By contrast, our ‘CR diff’ model predicts a change of q of about 0.2 in our highest star-forming galaxy when we neglect the secondary radio emission. If non-synchrotron cooling is disregarded, q changes by 0.2 to 0.4. Note that we refrain from analysing this effect as a function of gas surface density, because this observable strongly depends on the viewing angle of the galaxy, whereas the SFR is a more robust global quantity.

Finally, we examine three possible solutions for the problem of the discrepancy between the observed flat radio spectra and the theoretically expected steep radio synchrotron spectra in the calorimetric limit, i.e. that CR electrons lose their energy primarily to radiative (synchrotron and IC) processes before they can escape from the galaxies.

• Bremsstrahlung and Coulomb cooling processes are only able to affect low-energy electron spectra, and hence, yield flat radio spectra at low frequencies. However, these processes only play a minor role in the observed spectral flattening because (diffusive) escape losses are more important, especially at lower SFRs (see Fig. 12).

• We also find that advection losses cannot be primarily responsible for the observed low spectral indices because they only generate hard CR electron spectra in outflows. This implies flat radio spectra off-centre, which is in direct conflict with the observed flat central radio spectra (see Fig. 13).

• Our preferred solution is thermal free–free emission that starts to dominate the total radio spectrum at frequencies of several to tens of GHz and flattens the observed radio spectrum (see Figs 8 and 9). Thermal free–free emission is not only able to flatten the radio spectra of our simulated galaxies similar to NGC 253, M82 and NGC 2146 at high frequencies, but the involved absorption process, i.e. free–free absorption, coincides with the spectral flattening of the central spectrum of NGC 253.

This interpretation of thermal free–free emission being the main driver for flattening the radio spectra at high frequencies is corroborated by the observed spectrum of the central ∼450 pc of M82: subtracting this central emission region from the total spectrum yields a spectrum that is as steep as the radio synchrotron spectrum resulting from CR electrons close to the fully cooled limit, with spectral indices up to 1.1. Consequently, the flattening of the spectra at high frequencies must mainly originate from the central region of starbursts and can be reconciled with steep radio synchrotron spectra, that dominate the emission outside the dense, central starburst region and the total spectrum at lower frequencies, i.e. below ∼3, 10 or 20 GHz in our modelling of NGC 253, NGC 2146 and M82, respectively.

These insights into the radio emission from galaxies and the physics of the FRC necessarily require a spectral modelling of the different CR populations in full three-dimensional MHD–CR simulations, which demonstrates the power of this approach. While we confirm several findings of one-zone models (e.g. Lacki et al. 2010), our approach exhibits less free parameters and thus can be considered to be more predictive. Most importantly, because these simulations evolve the magnetic field and the CR energy density, we can use this modelling to link the simulated radio emission to observational data to quantify the effect of CR feedback on galaxy formation.

Clearly, the presented radio analysis of steady-state CR spectra need to be complemented by full spectral-dynamical simulations of CR protons (Girichidis et al. 2020) and electrons (Winner et al. 2019, 2020; Ogrodnik, Hanasz & Wóltański 2020). Moreover, the employed one-moment approach for CR transport (Pfrommer et al. 2017a) will be improved with a two-moment CR hydrodynamics model that is coupled to the Alfvén wave energy density, which delivers a realistic (spatially and temporally varying) CR diffusion coefficient in the self-confinement picture (Thomas & Pfrommer 2019, 2021; Thomas, Pfrommer & Pakmor 2021). Finally, pursuing MHD–CR simulations of galaxy formation in a cosmological setting (Buck et al. 2020; Hopkins et al. 2021) will be required to obtain realistic (bursty) star formation histories. These have different epochs of turbulent driving, which can potentially amplify the magnetic field in dwarf galaxies furthermore to come into equipartition with the thermal energy density so that it does not saturate at a sub-equipartition level (Pfrommer et al. 2021) and reaches the FRC.

ACKNOWLEDGEMENTS

The authors acknowledge support by the European Research Council under ERC-CoG grant CRAGSMAN-646955.

DATA AVAILABILITY

The data underlying this paper will be shared on a reasonable request to the corresponding author. The arepo code is publicly available.

Footnotes

1

This relation can be derived by approximating the synchrotron kernel with a Dirac delta distribution, but see also Appendix A1.

2

Note that this implies that we only sum over cells with an SFR |$\dot{M}_\star \gt 0$|⁠.

3

The energy dependence of the diffusion coefficient does not significantly change the radio spectra of these galaxies because radiative losses of CR electrons beat their diffusive losses which is different for CR protons that do not suffer radiative losses. For the considered snapshots, the effect of changing δ = 0.5 → 0.3 only changes |$L_{1.4\, \mathrm{GHz}}$| by less than 8 per cent.

4

The influence of advection losses in outflows is discussed in Section 5.3.

5

For simplicity, we use in the estimation of the characteristic electron energy entering the calculation of the time-scales in Fig. 12 the total magnetic field of each cell here. Adapting B_⊥ and calculating the resulting time-scale ratios for edge-on and face-on configurations, respectively, changes the ratios shown in Fig. 12 by at most a factor of 2.

APPENDIX A: RADIATION PROCESSES

A1 Synchrotron emission

Each relativistic electron with charge e that is accelerated by an ambient magnetic field B emits synchrotron radiation. The resulting power emitted per unit frequency ν is given by (Rybicki & Lightman 1986)

$$\begin{eqnarray*} P(\nu , \gamma _{\mathrm{e}})=\frac{\mathrm{d}E}{\mathrm{d}\nu \mathrm{d}t}=\frac{\sqrt{3}e^{3}B\sin \alpha _{\mathrm{pitch}}}{m_{\mathrm{e}}c^{2}}F(x), \end{eqnarray*}$$

(A1)

where α_pitch is the pitch angle, x = ν/ν_c, the critical frequency is |$\nu _\mathrm{c}=3/(4\pi)\, \gamma _{\mathrm{e}}^{3}\omega _{\rm B}\sin \alpha _\mathrm{pitch}$|⁠, and the frequency of gyration is ω_B = eB/(γ_em_ec). Furthermore, the dimensionless synchrotron kernel F(x) is defined as

$$\begin{eqnarray*} F(x)=x\intop _{x}^{\infty }K_{5/3}(\xi)\mathrm{d}\xi , \end{eqnarray*}$$

(A2)

where K_5/3 is the modified Bessel function of the order of 5/3. For a population of electrons with a distribution f_e(p), the total emissivity, i.e. the emitted energy per unit time, volume and frequency, is obtained by integrating over the electron distribution:

$$\begin{eqnarray*} j_{\nu }(\nu)= E\frac{\mathrm{d}N_{\gamma }}{\mathrm{d}\nu \mathrm{d}V \mathrm{d}t}=\frac{\sqrt{3}e^{3}B_{\perp }}{m_{\rm e}c^{2}}\intop _{0}^{\infty }f_{\mathrm{e}}(p)F(x)\mathrm{d}p, \end{eqnarray*}$$

(A3)

where B_⊥ is the component of the magnetic field perpendicular to the line of sight and the argument of the synchrotron kernel depends on the electron momentum via the dependence of the critical frequency on the Lorentz factor, |$\gamma _\rm {e}=\sqrt{1+p^2}$|⁠. Since the integration over the modified Bessel function is numerically expensive, we use an analytical approximation by Aharonian et al. (2010), which reads

$$\begin{eqnarray*} \widetilde{F}(x) &\approx & 2.15x^{1/3}(1+3.06x)^{1/6} \nonumber \\ && \times \,\,\frac{1+0.884x^{2/3}+0.471x^{4/3}}{1+1.64x^{2/3}+0.974x^{4/3}}{\rm e}^{-x}. \end{eqnarray*}$$

(A4)

This function peaks at x = 0.2858, but its expectation value lies at x = 2.13.

Figure A1.

We compare the central radio spectrum of our simulated NGC 253-like galaxy (see Table 2) to observational data by Kapińska et al. (2017) of the central region. We compare the spectra within three different radial regions (as indicated in the legend) and show the resulting thermal free–free emission (magenta) and radio synchrotron spectrum (green).

A2 Thermal free–free emission and absorption

Thermal electrons are deflected in the Coulomb field of ions of charge Ze and emit thermal free–free emission. The resulting emissivity from a medium with electron density n_e, ion density n_i, and temperature T is given by (Rybicki & Lightman 1986)

$$\begin{eqnarray*} j_{\nu ,\mathrm{ff}} &=& 6.8 \times 10^{-38}~\mathrm{erg\, s^{-1}\, cm^{-3}\, Hz^{-1}}\nonumber \\ &&\times Z^2 \frac{n_{\mathrm{e}}\, n_{\mathrm{i}}}{\mathrm{cm^{-6}}} \, T_1^{-0.5}\, {\rm e}^{-h\nu /(k_{\mathrm{B}} T)}\, \overline{g}_{\mathrm{ff}}, \end{eqnarray*}$$

(A5)

where |$T_1 = T/(1\, \mathrm{K})$|⁠, h is Planck’s constant, and k_B denotes Boltzmann’s constant. The corresponding absorption coefficient of free–free absorption in the Rayleigh–Jeans regime, i.e. for hν ≪ k_BT, reads (Rybicki & Lightman 1986)

$$\begin{eqnarray*} \kappa _{\mathrm{ff}}(\nu) = 0.018\, T_1^{-1.5}\, Z^2 \frac{n_{\mathrm{e}}\, n_{\mathrm{i}}}{\mathrm{cm^{-6}}} \, \left(\frac{\nu }{\mathrm{Hz}} \right)^{-2}\, \overline{g}_{\mathrm{ff}} \, \mathrm{cm^{-1}}. \end{eqnarray*}$$

(A6)

Here, the mean Gaunt factor |$\overline{g}_\mathrm{ff}$| in the ‘small-angle, classical region’ is given by (Novikov & Thorne 1973)

$$\begin{eqnarray*} \overline{g}_{\mathrm{ff}} = \frac{\sqrt{3}}{\pi } \ln \left[ \frac{1}{4\xi ^{5/2}Z} \left(\frac{k_{\mathrm{B}}T}{h\nu } \right) \left(\frac{k_{\mathrm{B}}T}{13.6\, \mathrm{eV}} \right)^{0.5} \right], \end{eqnarray*}$$

(A7)

where |$\xi =\exp (\gamma _\rm {E})\approx 1.781$| and |$\gamma _\rm {E}$| denotes Euler’s constant. We assume a temperature of T = 8000 K for the warm-ionized medium and approximate the electron density n_e by taking a constant fraction ξ_e of the electron density provided in our pressurized ISM model (Springel & Hernquist 2003) to parametrize the effect of radiation feedback in the dense star-forming/bursting regions and further assume n_e ≈ n_i and Z = 1. The optical depth due to free–free absorption is calculated via the integral of the absorption coefficient along the line of sight, i.e.

$$\begin{eqnarray*} \tau _{\mathrm{ff}}(\nu ,s_0,s_1) =\intop _{s_0}^{s_1} \kappa _{\mathrm{ff}} (\nu)\mathrm{d}s \end{eqnarray*}$$

(A8)

$$\begin{eqnarray*} = 0.018\, T_1^{-1.5}\, Z^2 \left(\frac{\nu }{\mathrm{Hz}} \right)^{-2}\, \overline{g}_{\mathrm{ff}} \intop _{s_0}^{s_1} \left(\frac{ n_{\mathrm{e}} }{\mathrm{cm^{-3}}}\right)^2 \mathrm{d}s, \nonumber\\ \end{eqnarray*}$$

(A9)

where |$\mathrm{d}s=\sin \phi \, \mathrm{d}y + \cos \phi \, \mathrm{d}z$| and ϕ denotes the inclination angle. Because for a constant temperature, j_{ν, ff}/κ_ff is constant along the line of sight and the solution of the radiative transfer equation (e.g. Rybicki & Lightman 1986) is straight forward and yields an observed intensity of free–free emission

$$\begin{eqnarray*} 4 \pi \, I_{\nu ,\mathrm{ff}}(s_1) = \frac{j_{\nu ,\mathrm{ff}}}{\kappa _{\mathrm{ff}}} \left\lbrace 1-\exp [-\tau _{\mathrm{ff}}(\nu ,s_0,s_1)] \right\rbrace . \end{eqnarray*}$$

(A10)

Also, the emitted synchrotron spectrum is affected by free–free absorption, which is typically relevant at low frequencies, as well as synchrotron self-absorption (SSA). The absorption coefficient and the corresponding optical depth of SSA reads (Rybicki & Lightman 1986)

$$\begin{eqnarray*} \kappa _{\mathrm{SSA}} (\nu)& =& -\frac{c^2}{8\pi \nu ^2} \intop _0^{\infty } \mathrm{d}p_{\mathrm{e}} P(\nu , p_{\mathrm{e}}) \left(1+p_{\mathrm{e}}^2 \right) \frac{\partial }{\partial p_{\mathrm{e}} }\left[\displaystyle\frac{f_{\mathrm{e}} (p_{\mathrm{e}})}{\beta \left(1+p_{\mathrm{e}} ^2\right)}\right] \nonumber\\ \end{eqnarray*}$$

(A11)

and

$$\begin{eqnarray*} \tau _{\mathrm{SSA}}(\nu ,s_0,s_1) &=\intop _{s_0}^{s_1} \kappa _{\mathrm{SSA}} (\nu)\mathrm{d}s. \end{eqnarray*}$$

(A12)

In order to calculate the absorbed synchrotron spectrum, we use the formal solution of the radiative transfer equation for an emissivity j_ν and the optical depth τ = τ_SSA + τ_ff., i.e.

$$\begin{eqnarray*} 4 \pi \, I_{\nu }(s_1) = \intop _{s_0}^{s_1} j_{\nu }(s) \exp [-\tau (\nu , s,s_1)] \mathrm{d}s. \end{eqnarray*}$$

(A13)

In practice, we solve equation (A13) by computing j_ν on slices that are oriented perpendicular to the line of sight s and equidistantly spaced along s through the simulation cube and cumulatively add the optical depth. If we observe a simulated galaxy edge-on, for instance, ds = dy and we obtain two-dimensional slices of the optical depth and the synchrotron emissivity in the x−z plane, i.e. τ(x, z) and j_ν(x, z). This results in a two-dimensional projection of the intensity I_ν(x, z) after adding up the slices from the front to the back of the simulation. This procedure also allows us to construct the absorbed spectrum from a different viewing angle by initially rotating the cube around the x-axis, which yields |$I_\nu (\boldsymbol {r}_\perp)$|⁠, where |$\boldsymbol {r}_\perp$| is the vector in the plane perpendicular to the line of sight. The observed flux density from an object emitting at a luminosity distance d is then obtained by integrating over the area A and solid angle Ω

$$\begin{eqnarray*} F_\nu = \frac{1}{4 \pi d^2} \int _\Omega \int _A I_\nu (\boldsymbol {r}_\perp) \mathrm{d}^2r_\perp \mathrm{d}\Omega . \end{eqnarray*}$$

(A14)

A3 Radio spectra for different central radii

The radio spectra of the central regions of our analysed galaxies shown in the middle panels of Figs 8 and 9 are calculated by adapting the radii summarized in Table 2. Here, we analyse the effect of choosing different central radii on the spectral shape of free–free-absorbed synchrotron spectra and the spectra of free–free emission. As an example, we show the central spectrum of NGC 253 in Fig. A1. Decreasing the radius of the central region of the galaxy yields a lower synchrotron flux, which additionally suffers from more absorption losses due to the higher central gas densities, leading to a stronger turn-over at low frequencies.

APPENDIX B: PARAMETER VARIATION OF THE FRC

In Fig. A2, we show the FRC with variations of the parameters of the model discussed in Section 4, where we adapt ζ_SN = 0.05 and |$B_0=10^{-10}\, \mathrm{G}$| (see Fig. 3). Instead, the simulations shown in the top panels of Fig. A2 have an initial magnetic field of |$B_0=10^{-12}\, \mathrm{G}$|⁠. Still, the observed FRC is reproduced both in our ‘CR diff’ and ‘CR adv’ models, respectively. There is only one outlier in the |$10^{12}~\rm {M}_\odot$| halo mass simulation in the ‘CR diff’ model, that overshoots the FRC beyond the observed scatter, at an SFR of 4.6 M_⊙ yr⁻¹. This illustrates the different dynamo actions taking place in the different simulations. As shown in Pfrommer et al. (2021), after t ≈ 1 Gyr, the averaged magnetic energy density of a simulation with |$M_{200}=10^{12}\, \mathrm{M_\odot }$| initiated with a lower magnetic field of B₀ = 10⁻¹² G manages to overtake the simulations that started with a higher initial magnetic field of B₀ = 10⁻¹⁰ G. However, for smaller haloes, the magnetic field grows at a smaller rate in the simulations with B₀ = 10⁻¹² G (Pfrommer et al. 2021) and hence, they tend to fall short of the FRC.

Figure A2.

Parameter variation of the FRC. In contrast to Fig. 3, here we adopt a lower initial magnetic field strength |$B_0 = 10^{-12}\, \mathrm{G}$| (top panels). In the bottom panels, we additionally increase the CR injection efficiency by a factor of 2, ζ_SN = 0.1. Generally, our radio luminosities match the observed scatter.

The lower panels of Fig. A2 show the FRC of the simulation model with |$B_0=10^{-12}\, \mathrm{G}$|⁠, but where CRs are injected with a higher acceleration efficiency and obtain a fraction of ζ_SN = 0.10 of the SN explosion energy, that is a factor of two larger than in our fiducial model. The combination of a lower initial magnetic field, but on the other hand a higher acceleration efficiency of CRs at SNe yields very similar results in comparison to our fiducial model and also matches the FRC within the observed scatter. However, we found in Paper II that this acceleration efficiency of 10 per cent is inconsistent with the observed FIR–gamma-ray relation, because it overpredicts the gamma-ray luminosity L_{0.1−100 GeV} in both the ‘CR diff’ and ‘CR adv’ models. This highlights the importance of comparing theoretical models with observations across all accessible electromagnetic frequencies.