Exploring the physics behind the non-thermal emission from star-forming galaxies detected in gamma rays

We revisit the $L_{\gamma}$–SFR correlation studied by K20 and reconstruct it with the most recent $\gamma$-ray fluxes, taken from the Fermi Large Area Telescope Fourth Source Catalogue Data Release 2 (4FGLR2, Ballet et al., 2020). This catalogue presents a new detection (Arp 299) compared to the previous release used by K20 (4FGL, Abdollahi et al., 2020).

We compute the $\gamma$-ray luminosities using the observed fluxes in the $0.1$–$100\,\mathrm{GeV}$ band and the luminosity distances $D_{L}$ from the self-consistent set reported in K20. For two objects, M33 and NGC 2403, data in the $0.1$–$100\,\mathrm{GeV}$ band are not reported in 4FGLR2; we thus compute their fluxes from the $0.1$–$800\,\mathrm{GeV}$ energy fluxes and the best-fitting spectral energy distributions (SEDs) provided by Ajello et al. (2020). We also include an estimate of the Milky Way (MW) $\gamma$-ray luminosity, computed from a CR propagation model (Strong et al., 2010; Ackermann et al., 2012). A power-law fit to these data as a function of the SFR, $L_{\mathrm{\gamma}}=A\dot{M}_{*}^{m}$, yields a value of $m=1.42\pm 0.17$ and $\log{A}=39.21\pm 0.15$, consistent with those reported in K20.

To obtain a radio-continuum luminosity sample centred at 1.4 GHz ($L_{\mathrm{1.4GHz}}$) as homogeneous as possible, we measured the radio fluxes directly from available images whenever possible. For eight galaxies (NGC 253, NGC 1068, NGC~2146, NGC 2403, M82, Arp220, Arp299, NGC 3424) continuum images at 1.4 GHz of Stokes-I are available in the NRAO/VLA Sky continuum Survey catalogue (NVSS, Condon et al., 1998). For all the radio images, with the exception of NGC 253, we obtained an average rms ($1\sigma$) background fluctuation value of $\sim$ 0.5 mJy $\rm beam^{-1}$, consistent with what was reported by Condon et al. (1998). We obtained the fluxes at 1.4 GHz ($F_{\mathrm{1.4\,GHz}}$) using the radio interferometry data reduction package MIRIAD (Sault et al., 1995), integrating the emission above $3\sigma$ to include all the flux related to the star formation and not only their central core and nucleus. We estimated the flux errors by integrating the emission above $2\sigma$ ($F_{\mathrm{2\sigma}}$) and subtracting it from the total flux. Whenever possible, we compared our fluxes with those reported by Yun et al. (2001), although these authors do not report flux errors. Our values are consistent with their nominal values.

For the particular case of NGC 253, we excluded two background sources not associated with its intrinsic emission (Kapińska et al., 2017). Due to the extension of NGC 2403, we calculated the radio flux within the optical isophote of $25\ \mathrm{mag\ arcsec}^{-2}$. We excluded patchy emission clearly disconnected from the main source and not associated with the optical emission of the galaxy, assuming that it does not come from the stellar populations. For NGC 4945 and Circinus, we calculate the 1.4 GHz flux using the same method but with the new Stokes-I continuum image provided by The Rapid Square Kilometre Array Pathfinder (ASKAP) Continuum Survey²²2https://research.csiro.au/casda/ (RACS, McConnell et al., 2020). For both galaxies we obtain an rms background of $\sim 0.3$ mJy/beam, consistent with the mean value $\sim 0.25$ mJy/beam reported by McConnell et al. (2020). The nearby galaxy Circinus shows extended radio lobes (Wilson et al., 2011); we report an average flux between the one including the lobes (1.94 Jy) and the flux excluding them (1.51 Jy) in order to be consistent with the calculated fluxes for unresolved galaxies. We note that, in these cases, the fluxes we calculate may be slightly contaminated by the emission produced by AGN jets or other outflows that may arise in SBGs.

For the closest galaxies, M31, M33, and the Magellanic Clouds (SMC and LMC), the flux calculation is more complex due to the large spatial extent of these galaxies and the presence of resolved structures. Therefore, we adopt the fluxes provided by detailed studies of their emission available in the literature (Dennison et al., 1975; Loiseau et al., 1987; Hughes et al., 2007).

The total IR luminosity in the 8–1000 $\mu\mathrm{m}$ band ($L_{\mathrm{IR}}$) is usually used as a proxy for the SFR. A power-law fit to $L_{\mathrm{1.4GHz}}$ as a function of $L_{\mathrm{IR}}$, $L_{\mathrm{1.4GHz}}=AL_{\mathrm{IR}}^{m}$, with $L_{\mathrm{IR}}$ from K20, yields $m=1.00\pm 0.04$, $\log{A}=11.41\pm 0.37$ and a dispersion of 0.2 dex, indicating a very tight relation. These values are in excellent agreement with the the ones reported by Bell (2003) for their sample of 162 SFGs, and are also consistent with the best-fit values for a much bigger sample of 1809 galaxies reported by Yun et al. (2001). A power-law fit to SFR data, $L_{\mathrm{1.4GHz}}=A\dot{M}_{*}^{m}$, yields a value of $m=1.13\pm 0.09$ and $\log{A}=21.14\pm 0.10$, with a larger dispersion of 0.34 dex concentrated close to $\log{\dot{M_{*}}}\sim 0$.

In the upper left panel of Fig. 1 we show the $L_{\mathrm{1.4GHz}}$–SFR correlation derived from the galaxies in our sample. The galaxies follow a tight correlation except for NGC 4945 which deviates from the relation marginally ($\Delta_{\dot{M_{*}}1.4}$ = 2.8). The latter may be due to an additional contribution linked to the hosted Seyfert 2 nucleus (Lenc & Tingay, 2009), or a misguided estimate of the SFR (see discussion in K20). For comparison we also show the estimated correlation by Yun et al. (2001). These authors derived their correlation estimating the SFR using the linear $L_{\mathrm{IR}}$–SFR relationship from Kennicutt (1998), which is not a good proxy for galaxies with low SFRs (Bell, 2003). This makes the $L_{\mathrm{1.4GHz}}$–SFR relationship presented by Yun et al. (2001) inconsistent with our fit in the low SFR range. We show the latter on upper right panel of Fig. 1.

The previous analysis shows that the luminosities we calculate for the restricted sample of galaxies detected in $\gamma$-rays are consistent with those obtained by other authors, and that they follow the global trend found for more general samples of SFGs observed in radio. This supports that the global radio emission of these galaxies does not present any particular characteristic with respect to that shown by more general samples of galaxies.

As a by product, a fit of a power law to the relation between $\gamma$-ray and radio luminosities ($L_{\gamma}=A(L_{\mathrm{1.4GHz}}/{10}^{21}\,\mathrm{W\,Hz^{-1}})^{m}$) gives $m=1.26\pm 0.10$, $\log{A}=39.04\pm 0.10$ and a dispersion of 0.5 dex. This slope is marginally steeper than the value previously reported by Ackermann et al. (2012), who finds a slope of 1.10 $\pm$ 0.05. These differences may be due to the size of our sample, which is almost twice as large as that of Ackermann et al. (2012), and also to the self-consistent and updated values of distances and fluxes we use. We show the resulting correlation on the lower left panel of Fig. 1, together with the Ackermann et al. (2012) relation.

At low radio frequencies ($\lesssim$ 300 MHz) some SFGs show a flattening of their spectra. This may be caused by partial absorption of the synchrotron emission by ionised gas, or by energy losses and propagation effects of CRs (Marvil et al., 2015; Chyży et al., 2018). Therefore, low-frequency data can be useful to infer the absorbing mechanisms and properties of the medium.

Eleven galaxies in our sample (M31, NGC 253, SMC, M33, LMC, NGC 2146, NGC 2403, M82, NGC 3424, Arp 299, Arp 220) have fluxes at 150 MHz available in the literature. For NGC 1068 we take the flux at 145 MHz and the spectral index given by Kuehr et al. (1981), and extrapolate the flux to 150 MHz. In all cases, we derive the corresponding luminosities at this frequency ($L_{\mathrm{150\,MHz}}$) using the K20 distances. We report these values in Table 1, together with the respective references. In the cases of M31 and Arp 299, we assumed a typical flux error of 10%.

A power-law fit to the data, $L_{\mathrm{150MHz}}=A\dot{M}_{*}^{m}$, yields a value of $m=0.98\pm 0.13$ and $\log{A}=21.7\pm 0.15$, with a dispersion of 0.5 dex. On the lower right panel of Fig. 1 we show the derived $L_{\mathrm{150MHz}}$–SFR correlation, together with those estimated by Wang et al. (2019) and Gürkan et al. (2018). Different works studied this correlation using much larger samples (Calistro Rivera et al., 2017; Read et al., 2018; Smith, D. J. B. et al., 2021), reporting values are between $m=1.07$ (Gürkan et al., 2018) and 1.37 (Wang et al., 2019). Our slope is compatible with that of Gürkan et al. (2018) but our data points fall systematically below the fit of these authors. Given the precision of radio data, this points to a systematic difference in the SFRs whose origin may be the different estimators used (UV–IR proxies in our case, multiwavelength photometric fits in theirs). Indeed these authors state that the SFRs provided by their method differ in $\sim 0.2$ dex from those obtained from H$\alpha$ luminosities, which are in general agreement with ours (see K20). A correction of 0.2 dex in SFR to their fit would bring it into agreement with our data, strongly suggesting that the difference in SFR estimation methods is responsible for this systematic departure. Regarding Wang et al. (2019), the value that they report is derived from a sample using only galaxies above their completeness limit ($\dot{M}_{*}\approx 3\,\mathrm{M}_{\odot}\,\mathrm{yr}^{-1}$). Their best fit using all the data, which has better statistics due to the larger number of galaxies at low SFRs, has a lower slope of $m\approx 1.1$ in better agreement with ours. Although the 150 MHz data are not homogeneous because the fluxes were obtained from different surveys, and with different criteria for the flux density integration, they deserve to be studied because they provide unique constraints to the properties of non-thermal emission of SFGs. We will keep in mind this caveat when interpreting results based on these data.

We model the galaxy emission region as a disk of radius $R$ and height $2H$ with constant gas density, and in which acceleration sites are uniformly distributed. We adopt $R=1\,\mathrm{kpc}$ and $H=0.2\,R$ suggested by K20. Following the standard prescriptions of diffusive shock acceleration at SNRs (e.g. Axford et al., 1977; Bell, 1978; Blandford & Ostriker, 1978), we assume that CRs are steadily injected in the system according to:

where $i=\mathrm{e,p}$ (electrons and protons), $\alpha=2.2$ is the spectral index (Lacki & Thompson, 2013), $E_{\mathrm{max}}$ is the maximum particle energy (assumed to be $10^{15}$ eV), and $\zeta_{\mathrm{p,e}}=(1,2)$ (Zirakashvili & Aharonian, 2007; Blasi, 2010). The normalisation constant $Q_{0,i}$ is computed as

assuming that the supernova rate of the galaxy is $\mathcal{R}_{\mathrm{SN}}=(83\,\mathrm{M}_{\odot})^{-1}\dot{M}_{*}$, consistent with the Chabrier (2003) initial mass function (IMF), and that a fraction $\xi_{\rm CR}=10\%$ of the supernova energy ($E_{\rm SN}=10^{51}\,\rm erg$) is converted in CRs. We adopt $E_{\mathrm{min,p}}$ = 1.2 GeV as the minimum proton energy. Moreover, this energy is distributed as $Q_{0,p}/Q_{0,e}=50$; such a proton-to-electron ratio is supported by multiwavelength modelling and observations (Torres, 2004; Merten et al., 2017; Peretti et al., 2019; Yoast-Hull et al., 2013). We note that a slightly different value compared to our benchmark assumption (within a factor two), only affects the emission from primary electrons (also within a factor two), which does not impact significantly our global results.

Young massive star clusters, as discussed in Morlino et al. (2021) can be additional acceleration sites for CRs where diffusive shock acceleration is inferred to be efficient possibly up to PeV. Even though these sources can play an important role as acceleration sites, it is not clear yet to which extent they can contribute to the bulk of galactic CRs where, so far, SNRs are believed to be dominant. Therefore, we work under the assumption that SNRs are dominating the injection of the bulk of CRs. Nevertheless, a second population of galactic accelerators such as young massive star clusters is unlikely to change our main results as far as the acceleration efficiency of both components is adjusted to produce a total CR power similar to that used here.

Electrons and protons experience different losses, which can be divided into radiative, non-radiative, and escape. Electrons suffer radiative losses via synchrotron, inverse Compton (IC) up-scattering of IR photons, and relativistic Bremsstrahlung (BS). Protons lose energy via inelastic $p$-$p$ collisions. For both particle species the non-radiative losses are due to ionisation. We refer to K20 for the exact expressions used to calculate the corresponding timescales. In particular, for the IC process we adopt the formalism developed by Khangulyan et al. (2014) which accounts for the complete differential cross-section in both Thomson and Klein – Nishina regimes. Finally, diffusion and advection regulate the particle escape from the emission region. We discuss the details of escape timescales in the following section.

The balance between injection and loss mechanisms leads to the following approximate steady-state solution for the particles distribution in the galaxy:

where $\dot{E}=E\,\tau_{\mathrm{loss}}^{-1}(E)$ and $\tau_{\mathrm{loss}}\leavevmode\nobreak\ =\leavevmode\nobreak\ 1/[\tau_{\mathrm{esc}}^{-1}+\tau_{\mathrm{cool}}^{-1}]$. Here $\tau_{\mathrm{esc}}$ is the characteristic CRs escape time computed as $\tau_{\mathrm{esc}}^{-1}=\tau_{\mathrm{adv}}^{-1}+\tau_{\mathrm{diff}}^{-1}$ and $\tau_{\mathrm{cool}}^{-1}=\sum_{j}\tau_{j}^{-1}$, where each $\tau_{j}$ represents the CRs cooling time by each different cooling mechanisms.

We note that this equation is almost equivalent to the standard leaky box solution, $N_{i}(E)=Q_{i}(E)\,\tau_{\rm loss}(E)$, with the benefit that it yields a more precise normalisation for injection indices $\alpha>2$ when cooling dominates over escape.

Here we present the formulae used to calculate $\gamma$ rays and secondary electrons spectra produced in $p$-$p$ interactions. Following Kelner et al. (2006), the emissivity of these secondary products $i=\gamma$,e can be summarised in two different formalisms depending on the energy $E_{i}$ of the resulting particle. In the following equations we denote $x_{i}=E_{i}/E_{p}$.

For $E_{\mathrm{i}}>100$ GeV the emissivity formula is:

with $\sigma_{\mathrm{pp}}$ the cross section for inelastic $p$-$p$ scattering,

where $E_{\mathrm{th}}=1.22$ GeV is the threshold energy for $\pi_{0}$ production, $L=\ln{(E/\mathrm{TeV})}$ and $\kappa=0.5$ is the process inelasticity. $\tilde{F}_{\gamma}$ and $\tilde{F}_{e}$ are given in equations (58–61) and (62–65) of Kelner et al. (2006), respectively.

For energies $E_{\mathrm{i}}\leq 100$ GeV we describe the emissivity by the $\delta$-functional approximation (Aharonian & Atoyan, 2000):

Here $n$ is the number density of the ISM, $m_{p}$ and $m_{\pi}$ are the masses of the proton and the pion, respectively, $E_{\pi}$ is the pion energy, $E_{i,\mathrm{min}}=E_{i}+m_{\pi}^{2}c^{4}/4E_{i}$, $K_{\pi}\approx 0.17$ is the fraction of kinetic energy transferred from the parent proton to the single pion, $\tilde{f_{\gamma}}(E_{\gamma}/E_{\pi})=1$ and $\tilde{f_{e}}(E_{e}/E_{\pi})$ are defined in equations (36–39) from Kelner et al. (2006) (see also Appendix C).

Part of the ionised gas in SFGs fills uniformly the ISM, whereas another part is concentrated in Hii regions randomly distributed over the disc. This gas emits free-free radiation in the radio band. To estimate this emission, we consider that the ionised gas is homogeneously distributed in the same volume as the CRs. In addition, this gas can be optically thick to the low-frequency (thermal and non-thermal) radio emission. Our model has two free parameters that determine the impact of the free-free absorption and emission processes: the electron temperature $T_{e}$ and the ionisation fraction $\rho$, defined here as the ratio between ionised and total number densities. Thus, we compute the absorbed free-free thermal luminosity at a frequency $\nu$, $L_{\mathrm{ff}}(\nu)$, following the Olnon (1975) formalism for an homogeneous cylinder seen face-on, valid in the Rayleigh-Jeans regime ($h\nu\ll kT_{e}$, where $h$ and $k$ are the Planck and Boltzmann constants, respectively),

where $\tau_{\mathrm{ff}}=\alpha_{\mathrm{ff}}(\nu)\,l$ is the optical depth, $l=2H$ is the length of the ionised region along the line of sight and $\alpha_{\mathrm{ff}}$ is the absorption coefficient for photons of frequency $\nu$. This coefficient is given by Rybicki & Lightman (1986),

where $T_{\mathrm{e}}$ the ionised gas temperature and $\bar{g}_{\mathrm{ff}}$ is the mean Gaunt factor (e.g. Leitherer et al., 1995),

The warm neutral and ionised medium phases of the ISM are inferred to have an average temperature of $10^{4}$ K. However, depending on the fraction of cooler neutral or hot ionised component, the effective ISM temperature could possibly differ by about one order of magnitude (Quireza et al., 2006). Therefore, we assume $T_{\mathrm{e}}=10^{4}$ K as a reference value and explore the impact of an order of magnitude change in its value.

We link the IR luminosity to the SFR by using the relation of Bell (2003), adapted to a Chabrier (2003) IMF:

with $\epsilon=0.79$ (Crain et al., 2010) and $L_{\mathrm{IR}}$ the total infrared luminosity between 8 and 1000 $\mu$m. For the calculation of the IC interactions we introduce a dilution factor of the energy density of the IR radiation field with respect to the one corresponding to a black body. For a disc geometry this factor is $\Sigma=L_{\rm IR}/(8\pi R^{2}\sigma_{\mathrm{SB}}T^{4})$, where $\sigma_{\mathrm{SB}}$ is the Stephan-Boltzmann constant.

Finally, the scaling of the density $n$ with the SFR can be set by using the Kennicutt-Schmidt (K-S) law (Kennicutt, 1998; Kennicutt & Evans, 2012; de los Reyes & Kennicutt, 2019). This law gives the correlation between the surface densities of SFR ($\Sigma_{\mathrm{SFR}}$) and cold gas mass ($\Sigma_{\mathrm{gas}}$) of galaxies. Assuming that the emitter is a cylinder of radius $R$ and thickness $2H$ (K20), we get

Robishaw et al. (2008) showed that strong magnetic fields ($B$ of the order of mG – hundreds times stronger than in our Galaxy) pervade the ISM of high SFR objects, such as ULIRGs, suggesting a strong dependence of $B$ on the level of source activity. Therefore, a detailed modelling of the transport of CR-electrons and their associated non-thermal radio emission requires a careful treatment of the magnetic field and its dependence on the SFR. For this reason, we improved the model developed in K20 proposing the following scaling,

where the normalisation factor $B_{0}$ is a free parameter and for the index $\beta$ we explore three possible scenarios, $\beta=0.3,0.5,0.7$, based on equipartition arguments described below.

The weakest scaling is obtained relying on a stationary condition of the MHD turbulence in the ISM plasma. The MHD turbulence reaches its stationary state when the magnetic field growth produced by turbulence-induced field line stretching is balanced by the back-reaction of the magnetic field itself (Groves et al., 2003). Such a condition results in equipartition between the energy density of the magnetic field and the kinetic energy density of the gas. This yields a scaling $B\propto\Sigma_{\rm g}^{0.5}$, which is in good agreement with numerical simulations that report exponent values in the range 0.4–0.6 (Groves et al., 2003). By applying the K-S scaling ($\Sigma_{\mathrm{SFR}}\propto\Sigma_{\mathrm{g}}^{1.41}$) one can retrieve the index $\beta\approx 0.3$. The second scenario is based on an equipartition condition between magnetic field and starlight, namely $U_{\rm B}\approx U_{\rm rad}$. This condition, often quoted in the literature (e.g. Murphy, 2009; Lacki et al., 2010; Yoast-Hull et al., 2016), is motivated by the assumption that the radiation pressure drives the turbulence in the ISM, which in turn determines the strength of the magnetic field. The linear relation between the SFR and the starlight results in the scaling $\beta\approx 0.5$. The third scenario considered relies on the stability of the galactic disc. As discussed by Lacki et al. (2010) (see also Parker, 1966), the equipartition between the gas pressure at the mid-plane of the disc ($\pi G\Sigma_{g}^{2}$) and the magnetic field pressure ($B^{2}/8\pi$) provides a natural limit for the magnetic field strength; higher values would make the galactic disk unstable. This condition, once the Kennicutt scaling is applied, leads to $\beta=0.7$.

We assumed $\beta=0.3$ as the fiducial value for the magnetic field index. In addition, we explored other prescriptions with a stronger $B$ scaling ($\beta=0.5,0.7$), which could be consistent with a more relevant molecular phase of the ISM in environments with high SFR. We adopted $B_{0}=100\leavevmode\nobreak\ \mu{\rm G}$ as a fiducial normalisation, which yields $B$-values in agreement with both the expected ones for normal starbursts galaxies (see also Peretti et al., 2019) and with the one inferred for the Galactic central region (see Ferrière, 2009; Lacki & Thompson, 2013). We also explored different values of $B_{0}$ across an order of magnitude to cover the different inferred values for the average magnetic field in SFGs, ranging from $\sim\mu$G in low SFR galaxies (Ferrière, 2010; Jurusik et al., 2014) to $\sim$mG for ULIRGs (Robishaw et al., 2008).

Winds are ubiquitous in SFGs and are possibly one of the most spectacular consequences of their star forming activity. They are capable of driving mass out of galaxies at a rate of several $\rm M_{\odot}\,yr^{-1}$ with velocities up to $10^{3}\,\rm km\,s^{-1}$ (see Veilleux et al., 2005). The mechanism(s) powering galactic winds are likely to depend on the SFR. In SBGs they arise from the combined effect of mechanical energy and heating produced by supernovae, and radiation pressure from young stars (Zhang, 2018). On the other hand, these effects are unlikely to be the main responsible for wind launching in galaxies with mild SFR; instead, the CR pressure seems to play the key role. In particular, the gradient of escaping CRs can be strong enough to overcome the gravitational pull and lift the ISM into the halo (Breitschwerdt et al., 1991; Recchia et al., 2016). In addition, in the case of SBGs, CRs cool more efficiently via $p$-$p$ collisions and their contribution to powering galactic winds could be minor, but this is still under debate given the importance of the CR pressure exerted on the cold interstellar gas (Bustard & Zweibel, 2021).

The effect of galactic winds on the transport of CRs needs to be taken into account. We thus consider that particles are advected away from the galaxy with an energy-independent characteristic time:

where $v_{\rm adv}$ is the SFR-dependent advection speed.

We divide the whole galaxy population in two classes: 1) high SFR and 2) mild SFR objects. We assume that objects with a star formation higher than the MW, $\dot{M}_{*}\gtrsim 2\,\rm M_{\odot}\,yr^{-1}$, are characterised by a starburst-driven wind, $v_{\mathrm{sw}}=400\,\mathrm{km\,s^{-1}}$ (see also K20, Yoast-Hull et al., 2013; Peretti et al., 2019). For objects with $\dot{M}_{*}<2\,\rm M_{\odot}\,yr^{-1}$, we assume that the galactic wind is CR-driven. In this scenario, Recchia et al. (2016) have shown that, at the base of the wind, the dominant advection velocity is the Alfvén speed. Therefore, we consider that particles are advected away from the galactic box at the local Alfvén speed, $v_{\mathrm{A}}=B/\sqrt{4\pi\,m_{p}\,n_{i}}$. Putting all together, we define a global advection velocity over the entire SFR range as

Diffusion is a fundamental aspect of CR physics. At a microscopic level, this phenomenon consists in resonant pitch-angle scattering of CRs with the local average magnetic field, taking place when the Larmor radius approaches the wavelength of local magnetic field disturbances (e.g. Alfvén waves, as discussed in Blasi, 2013). At a macroscopic level, this small scale interaction results in a diffusive motion of CRs, where the confinement properties of the system are parametrised by the diffusion coefficient, $D$. This coefficient holds all information that effectively characterises diffusion, such as the turbulence properties of the medium.

We compute the diffusion coefficient adopting the quasi-linear theory formalism

where $v(E)$ is the particle speed, $r_{L}$ is the Larmor radius, $l_{c}$ represents the coherence length of the magnetic field, and $\delta$ characterises the nature of the turbulence. In particular, we have $\delta=1/3$ for the Kolmogorov turbulence, and $\delta=1,\,0.5$ for Bohm and Kraichnan. In this work we assume the Kolmogorov-like turbulence, $\delta=1/3$, as our benchmark case, in agreement with current observations of the boron-to-carbon ratio in the Galaxy (Aguilar et al., 2016). We assume that supernovae are responsible for the injection of the turbulence in the system. Therefore, we compute $l_{c}=(2H\pi R^{2}\mathcal{R}_{\mathrm{SN}}^{-1}\tau_{\mathrm{SNR}}^{-1})^{1/3}$ as the average separation between SNRs in the galactic box, where we adopt a typical value for the SNR lifetime $\tau_{\mathrm{SNR}}=3\times 10^{4}\,\mathrm{yr}$ (Truelove & McKee, 1999). We also impose that such a parameter cannot be smaller than the typical dimension of a middle age SNR, namely $l_{c}\gtrsim 10\rm\,pc$.

The diffusion coefficient defined in Eq. (16) is assumed to be homogeneous in the whole galactic box. This allows us to define the characteristic timescale, $\tau_{\rm D}$, in which particles diffuse away from the galaxy

We give the explicit dependency of $\tau_{D}$ on the SFR in Appendix B.

The transport of electrons and protons in SFGs is strongly influenced by the average properties of the ISM, and in turn by their SFR. In this section we present the main transport properties of the benchmark scenario throughout the SFR range, analysing separately electrons and protons. In order to highlight the evolution of the transport condition with the SFR, we show in the upper (lower) panels of Fig. 2 the timescales of the relevant processes for electrons (protons) at two different values of the SFR: 0.1 $\rm M_{\odot}\,yr^{-1}$ and 10 $\rm M_{\odot}\,yr^{-1}$.

We define $L_{\gamma}$ as the integrated luminosity in the 0.1–100 GeV energy band, and we show in Fig. 3 the $L_{\gamma}$–SFR relation obtained for the benchmark scenario. The main contribution to $L_{\gamma}$ comes from the $p$-$p$ interaction, while $\gamma$-rays of leptonic origin are subdominant (see Fig. 8). In agreement with K20, the shape of the $L_{\gamma}$–SFR relation at low SFR is influenced by the escape of most of the injected protons, whereas at high SFR sources approach the calorimetric limit, where the correlation becomes a linear scaling. We also show the results obtained assuming a unique nature of the advection speed throughout the whole SFR range, namely $v_{\rm adv}=v_{\mathrm{A}}$ or $v_{\rm adv}=v_{\mathrm{sw}}$. Since $v_{\mathrm{sw}}>v_{\mathrm{A}}$ ($v_{\mathrm{A}}\gtrsim 50$ km s^-1 for the fiducial scenario), the CR escape is more efficient when $v_{\rm adv}=v_{\mathrm{sw}}$ and consequently the $\gamma$-ray luminosity is reduced. Our benchmark case is in good agreement with the best fit of the observed sample.

Variations of $\beta$, $B_{0}$ and $T_{\mathrm{e}}$ do not affect the emission of protons and, consequently, the $L_{\gamma}$–SFR correlation is not altered. On the contrary, the variation of the ionisation fraction $\rho$ (scenarios (5), (6)) modifies the CR escape through $v_{\mathrm{A}}$ at low SFRs, affecting the total source luminosity.

As described in Sec. 3.3, we explore values for the ionisation fraction down to $10^{-4}$. While these values can characterise powerful starbursts such as ULIRGs, they are unlikely to be physical in standard SFGs with low SFR. Indeed, such low ionisation fraction leads to an Alfvén speed larger than the average velocity of starburst superwinds resulting in an advection dominated transport at odd with observations of CRs (Blasi, 2013) and $\gamma$-rays (Ajello et al., 2020). The scenario $\rho=0.1$ does not show any relevant modifications with respect to the benchmark scenario since the dominant escape mechanism is diffusion.

As mentioned in Sec. 4, we improved the model of K20 introducing a SFR-dependent diffusion coefficient and a two-fold nature of the wind, while both diffusion coefficient and the wind velocity do not depend on the SFR in K20. Despite of these differences, the results presented here are compatible with the benchmark scenario discussed in K20.

The $L_{\gamma}$–SFR relation obtained here is in good agreement with the trend of data, and consistent with the results presented in K20. The model approaches the calorimetric limit (limit of the yellow shaded region) for $\dot{M}_{*}\geq 30\,M_{\odot}\rm yr^{-1}$, above which galaxies show a calorimetric behaviour. NGC 3424 and 4945 luminosities are incompatible with this limit by 4.1 and 6.4 times their errors. As discussed in K20, the former probably harbours a low-luminosity AGN, whereas the SFR of the latter may have been underestimated due to its obscuration and inclination to the line of sight.

The emission at 1.4 GHz is set by three processes: synchrotron radiation from both primary and secondary electrons, free-free (thermal) emission of the ionised gas, and free-free absorption produced by this same gas. Synchrotron radiation dominates the $L_{1.4\mathrm{GHz}}$ throughout the entire SFR range, as shown in Fig. 4 (left panel), while free-free emission and absorption become relevant only at high SFR. Eventually the thermal contribution becomes comparable to that of synchrotron at $\dot{M}_{*}>100\,\rm M_{\odot}\,yr^{-1}$ (see Fig 8).

The emission of primary electrons dominates the $L_{1.4\rm GHz}$–SFR correlation in the low SFR range ($\lesssim 1\,\rm M_{\odot}\,yr^{-1}$), whereas synchrotron from secondaries gradually increases with the SFR and overcomes that of primaries at high SFR (see respectively green dot–dashed and dashed lines in left panel of Fig. 4). Such a behaviour is a direct consequence of the transport condition of the parent protons. In particular, the higher the SFR, the higher the gas density and so the rate of interaction of protons.

The secondary–to–primary ratio increases with the SFR as a consequence of the increasing level of calorimetry up to $\gtrsim 30\,\rm M_{\odot}\,yr^{-1}$, where full calorimetry is achieved and the secondary–to–primary ratio becomes constant. This result is found to be in agreement with those of Peretti et al. (2019).

In SFGs with SFR $\lesssim 10^{2}\,\rm M_{\odot}\,yr^{-1}$, $L_{1.4\rm GHz}$ is dominated by the synchrotron emission, hereafter $L_{\rm syn}^{1.4}$. In order to understand the shape of the correlation, we note that $L_{\rm syn}^{1.4}$ is dominated by the emission of $\sim$ GeV electrons, independently of the SFR. This electron-energy band is characterised by the competition of several cooling and escape mechanisms, so that the slope of the particle energy distribution can range from $\alpha-1$ when ionisation dominates to $\alpha+1$ when synchrotron dominates. In addition, secondary electron emission depends on the SFR through both the density of the medium and the parent proton transport (see Sec. 3.2). This leads to a complex dependence of $L_{\rm syn}^{1.4}$ on the SFR, which we parametrise as

The relevance of thermal processes increases for higher SFR as a consequence of the higher density of the ISM and our assumption of constant ionisation fraction $\rho$ throughout the whole SFR range. In particular, the thermal contribution to $L_{1.4\mathrm{GHz}}$ is at the level of $\sim 5\%$ for low SFR galaxies; this is consistent with values reported by Condon (1992) for the thermal emission in the MW and M82. At high SFRs ($\gtrsim 100\,\rm M_{\odot}\,yr^{-1}$) the thermal contribution is $\sim 50\%$, and the free-free absorption also becomes significant. Nonetheless, assumptions in the ionisation fraction $\rho$ impact on the relative importance of the free-free processes, as we will discuss below.

Finally, we compare our benchmark model with the data presented in Table 1. We highlight that the coexistence of both primary and secondary synchrotron electron components is necessary to reproduce the slope of the $L_{1.4\rm GHz}$–SFR correlation. On the right panel of Fig. 4 we present observational data in the $L_{\gamma}$–$L_{1.4\rm GHz}$ plane and compare them with the outcomes of our model; we obtain a good agreement throughout the whole SFR range.

In this section, we discuss the impact of parameter variations on the $L_{1.4\rm GHz}$–SFR correlation.

The physical processes governing the emission at 150 MHz are the same as at 1.4 GHz (Sec. 5.3). However, at this frequency free-free absorption plays a crucial role. In fact, the turnover frequency where $\tau_{\mathrm{ff}}>1$ increases up to 150 MHz for $\dot{M}_{*}\sim 10\leavevmode\nobreak\ M_{\odot}/\mathrm{yr}$. Therefore, the 150-MHz synchrotron emission is completely absorbed at high SFRs in our benchmark scenario. This is reflected in the $L_{150\mathrm{MHz}}$–SFR correlation shown in Fig 6. This scenario departs considerably from the observed trend of data at high SFRs and cannot reproduce the emission of the ULIRGs. Such a strong absorption is the result of the ionisation fraction assumed, $\rho\sim 5\%$. In addition, the calorimetric limit determined for the benchmark scenario also shows that it can hardly match the observed fluxes for SFR $\gtrsim 10^{2}\,\rm M_{\odot}\,yr^{-1}$.

Although the effects of the variation of parameters in the $L_{150\mathrm{MHz}}$–SFR are similar to the $L_{1.4\rm GHz}$–SFR correlation, some details are required to figure out if it is possible to match the trend of data with some model prediction under different parametric assumptions.

The assumption $\beta=0.5$, as in the case at 1.4 GHz, increases the slope of both the correlation and the calorimetric limit. Moreover, ionisation losses are dynamically important for electrons emitting synchrotron at 150 MHz. In addition, different assumptions for $\rho$ and $T_{\mathrm{e}}$ can also mitigate or enhance the absorption as evidenced in the right bottom panel of Fig 5. As described above for 1.4 GHz, it is hard to find a unique parameter configuration that can accommodate the entire SFR range without being at odds with observations. In particular, we observe that the introduction of a dependence on the SFR for $\rho$ and/or $T_{\mathrm{e}}$ would provide a better fit at high SFR. Such a lack of emission at low radio frequencies was also reported by Schober et al. (2017) who found a similar result studying the correlation at 200 MHz.