© ESO 2022

1. Introduction

Star-forming galaxies (SFGs) are one of the extragalactic γ-ray source classes detected at gigaelectronvolt (GeV) energies by the Fermi-LAT telescope (Abdo et al. 2010a,b; Lenain et al. 2010; Ackermann et al. 2012). Their emission is believed to be associated with young stellar populations, which explains its correlation with star formation rate (SFR). In particular, γ rays are thought to be mostly produced in interactions between hadronic cosmic rays (CRs) with particles of the interstellar medium (ISM), while CRs are accelerated at supernova remnant (SNR) shocks (Jokipii & Morfill 1985; Blasi 2013; Bustard et al. 2017). In addition to SNRs, other sources related to star-forming regions, such as young massive star clusters or microquasars, could also play an important role in CR acceleration, possibly up to petaelectronvolt (PeV) energies (see e.g. Bykov & Toptygin 2001; Aharonian et al. 2019; Bykov et al. 2020; Escobar et al. 2021; Morlino et al. 2021). These acceleration sites are abundant in SFGs and are expected to be more numerous in sources with higher SFRs. Similarly, higher gas densities are found in higher SFR environments since the gas is the seed of star formation. As a result, an increase in the γ-ray emission is expected for higher star formation rates.

The main evidence supporting this interpretation is the observed correlation between the integrated γ-ray luminosity (L_γ) of SFGs and their SFRs (Ackermann et al. 2012; Ajello et al. 2020). The latter, in particular, is typically estimated from tracers such as the total infrared (IR) luminosity (L_IR, 8 − 1000 μm). This relation was predicted by Thompson et al. (2007) and Lacki et al. (2010) in terms of the γ-ray flux and the gas surface density of SFGs. Several works have endeavoured to explain the high-energy emission of these objects individually (e.g. Völk et al. 1996; Blom et al. 1999; Romero & Torres 2003; Domingo-Santamaría & Torres 2005; de Cea del Pozo et al. 2009; Persic & Rephaeli 2010; Lacki et al. 2011; Rephaeli & Persic 2014; Krumholz et al. 2020) and others as a result of the common properties of the population of SFGs (e.g. Lacki et al. 2010; Pfrommer et al. 2017; Zhang et al. 2019; Peretti et al. 2020; Kornecki et al. 2020). There is a general consensus on the increasing dynamical relevance of radiative energy losses for higher SFRs, so that the CR transport in these environments is inferred to approach the calorimetric limit; specifically, all injected particles cool before being able to escape. In particular, while electron calorimetry is believed to be reached already in mild starbursts, SFR ≳ 5 M_⊙ yr⁻¹, the SFR range where proton calorimetry might be achieved is still under debate due to our current poor knowledge of the escape mechanisms.

Star forming galaxies have been observed and modelled throughout the whole electromagnetic spectrum (e.g. Völk et al. 1996; Persic & Rephaeli 2002, 2014; Persic et al. 2008; de Cea del Pozo et al. 2009; Rephaeli et al. 2010; Yoast-Hull et al. 2013; Peretti et al. 2019). Diffuse thermal and non-thermal emission from these objects has been observed in radio frequencies (e.g. Condon 1992; Cox et al. 1988; Loiseau et al. 1987; Strickland et al. 2002; Kapińska et al. 2017), X-rays (Bauer et al. 2008; Sarkar et al. 2016), and even very high-energy γ rays (Acciari et al. 2009; Acero et al. 2009).

Non-thermal radio and X-rays are produced by relativistic primary and secondary CR electrons, where primaries are injected together with protons at acceleration sites (see e.g. Cristofari et al. 2021; Aharonian et al. 2019 in the context of SNRs and star clusters, respectively) and secondaries are a by-product of the CR interaction with the ISM.

CR electrons in galactic magnetic fields (B) emit synchrotron radiation detectable mostly at radio frequencies. A connection between the synchrotron radiation and supernova rate was one of the first attempts to explain the linear correlation between the observed luminosity at 1.4 GHz (L_1.4 GHz) and far-IR (FIR) luminosity (van der Kruit 1971; Condon 1992; Yun et al. 2001). However, the nature of this link is not trivial in low-SFR galaxies, due to their low efficiency in converting starlight into IR radiation, and the predominance of escape as the main driver of electron transport (Condon 1992; Yun et al. 2001; Bell 2003). This missing energy causes both radio and IR luminosities to underestimate the SFR for galaxies with low levels of IR emission. Bell (2003) argued that these two effects conspire to achieve a linear relation between IR luminosity and L_1.4 GHz. At high SFRs the efficient cooling of protons typical of starburst galaxies (SBGs) leads to a copious production of secondary leptons whose synchrotron emission is inferred to overcome that of primaries (e.g. Torres 2004; Domingo-Santamaría & Torres 2005; Lacki et al. 2010). This could possibly impact the global shape of the correlation.

In addition to synchrotron processes, thermal processes such as free–free emission and absorption also characterise the radio emission of SFGs. Numerous young stars create several H II regions, whose hot ionised gas radiates thermally and contributes to the total radio emission. A spectral identification of such free–free emission is in general very complex due to the presence of other components such as synchrotron, dominating at GHz, and the FIR thermal peak produced by the dust, which occurs above a few tens of GHz. Many works have discussed the fraction of the thermal contribution in this band with respect to the non-thermal emission (e.g. Bell 2003; Klein et al. 2018). However, no definitive conclusion has been reached yet due to the complexity of this task.

Focusing on the synchrotron emission, many studies explored the SFR–luminosity correlation at lower radio frequencies, in particular at 150 MHz (e.g. Gürkan et al. 2018; Wang et al. 2019; Smith et al. 2021), where the free–free absorption might become severe. Typically, the free-free absorption significantly reduces the emission up to ∼10 − 10² MHz in galaxies with moderate SFRs (Schober et al. 2017; Chyży et al. 2018), and up to a few GHz in ultraluminous IR galaxies (ULIRGs; Condon et al. 1991; Clemens et al. 2010).

Although most of the radio emission is expected from CRs populating galactic discs, some additional radio emission produced in SFGs might have a different origin. In particular, a non-negligible contribution might come from individual sources such as SNRs (Lisenfeld & Völk 2000), activity of a coexisting active galactic nucleus (AGN; Mancuso et al. 2017), or particles populating the galactic halo as a consequence of the escape from the galactic disc via diffusion or advection in winds (Heesen et al. 2009; Romero et al. 2018). A smaller contribution of emission coming from pulsar wind nebulae (Ohm & Hinton 2013) and millisecond pulsars (Sudoh et al. 2021) is also possible; in particular, the latter can become important at low frequencies and low SFRs. Nevertheless, the contribution of these objects to the diffuse emission is still unclear.

Lacki et al. (2010) and Martin (2014) have explored the transport of CRs throughout the SFR range, studying at the same time both the radio and γ-ray correlations with SFR. This multiwavelength approach allowed better constraints on the role of CRs in shaping these correlations highlighting, in particular, the key role played by protons in the production of both γ rays and secondary electrons. The latter are inferred to be responsible for most of the radio emission at high SFRs, where protons are expected to be efficiently confined. Despite these recent improvements, our knowledge of the CR transport properties in SFGs remains poor and several questions are still open: How does the CR escape from galaxies evolve with the SFR and what is its impact on calorimetry? What is the interplay between primary and secondary electrons in shaping the radio luminosity of SFGs? Can the luminosity–SFR correlation help us explore possible dependences of magnetic and radiation fields in SFGs? What is the actual role and/or contamination of free–free emission and absorption in SFGs? Can it help us in exploring alternative cosmic acceleration sites?

In this work we address these issues by exploring the non-thermal and thermal emissions in galaxies detected in γ rays, and their relation with the SFR at different bands (GeV γ rays and 1.4 GHz and 150 MHz in radio). We provide for the first time a nearly homogeneous set of observations, and compute the general trend for the L_γ, L_1.4 GHz, and L_150 MHz–SFR correlations for all SFGs detected in γ rays so far. In particular, we update the L_γ − L_1.4 GHz correlation presented by Ackermann et al. (2012). In order to understand the nature of these correlations, and improving on the work of Kornecki et al. (2020) (hereafter K20), we develop a detailed model of the non-thermal emission also taking into account the thermal component produced by the ionised gas and the associated absorption. Different from Lacki et al. (2010), we allow for advection escape in the entire range of SFRs. In addition, we attempt to capture the main properties of CR diffusion proposing a model for CR escape adequate to the source activity level, comprising SFR–dependent diffusion coefficient. In addition, we explore different scaling laws of B with the SFR. Unlike Martin (2014), we focus on the singular contribution of primary and secondary electrons to the radio correlations. This allows us to explore the impact of calorimetric transport on such correlations. Additionally, we include and discuss the role of free-free emission and absorption. We analyse carefully the effects of model parameters (as the total ionisation fraction of the ISM) and scaling relations needed to describe the emission at the three different wavelengths mentioned.

The paper is organised as follows. In Sect. 2 we present the data set used and the luminosity–SFR correlations. In Sects. 3 and 4 we present the model and the scale relationships used, taking as a starting point those proposed in K20. In Sect. 5 we present our main results, test the robustness of the model, and compare its predictions with observations to assess its reliability and delimit the most probable region of the parameter space. We also discuss our results and compare them with the literature in this section. Finally, in Sect. 6 we present our conclusions.

2. The luminosity–SFR relation

We address the physics behind the observed luminosity–SFR relations in a panchromatic fashion, specifically in the radio and γ-ray bands that trace non-thermal emission from relativistic particles. The existence of a tight correlation between the SFR of a galaxy and its luminosities in these bands (L_γ and L_radio, respectively) spanning more than four decades in SFR, indicates that non-thermal emission processes are common to all SFGs. The most compelling interpretation for these relations is that they arise from the emission of the population of cosmic rays within SFGs. Other scenarios such as the emission in jets or super-winds have been associated with non-thermal emission, but there is no evidence supporting that they are ubiquitous in SFGs.

Among these two bands, the sample of galaxies detected in γ rays is the most restrictive. We therefore focus on this sample, that comprises 14 objects. These SFGs have a steady γ-ray emission (Ackermann et al. 2012; Ajello et al. 2020), with the exception of the Circinus Galaxy, which shows marginal evidence of variability (Guo et al. 2019). This favours the interpretation that their high-energy emission is associated with their stellar populations instead of jets. Due to their different sizes and distances, and also because of the different resolving power of γ-ray and radio telescopes, some of these galaxies are detected as point sources, whereas others are resolved. We do not take into account the structure of the resolved sources, because our main goal is to study the link between the diffuse emission of the galaxies and their global SFR.

In the following sections, we construct our sample by taking from the literature the observed fluxes of each of the 14 galaxies currently detected in γ rays. We use the homogeneous set of distances and SFRs determined by K20 to calculate their radio and γ-ray luminosities in a self-consistent way. We also adopt the K20 SFR values for the galaxies in our sample. These SFRs were computed with L_IR-independent methods, when possible, in order to obtain a reliable SFR estimate for unobscured systems. In Table 1 we present the data on distances, and fluxes and luminosities in the relevant bands, together with their respective errors and references.

2.1. The L_γ–SFR relation

We revisited the L_γ–SFR correlation studied by K20, and reconstructed it with the most recent γ-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 computed the γ-ray luminosities using the observed fluxes in the 0.1 − 100 GeV band and the luminosity distances D_L from the self-consistent set reported in K20. For two objects, M 33 and NGC 2403, data in the 0.1 − 100 GeV band are not reported in 4FGLR2; we thus computed their fluxes from the 0.1 − 800 GeV energy fluxes and the best-fitting spectral energy distributions (SEDs) provided by Ajello et al. (2020). We also included an estimate of the Milky Way (MW) γ-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_{\gamma }=A{\dot{M}}_{\ast }^m$, yields a value of m = 1.42 ± 0.17 and log A = 39.21 ± 0.15, consistent with those reported in K20.

2.2. The L_1.4 GHz–SFR relation

To obtain a radio-continuum luminosity sample centred at 1.4 GHz (L_1.4 GHz) that is 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, M 82, Arp 220, Arp 299, 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σ) background fluctuation value of ∼0.5 mJy beam⁻¹, consistent with the value reported by Condon et al. (1998). We obtained the fluxes at 1.4 GHz (F_1.4 GHz) using the radio interferometry data reduction package MIRIAD (Sault et al. 1995), integrating the emission above 3σ 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σ (F_2σ) 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 mag arcsec⁻². 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 calculated 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¹ (RACS, McConnell et al. 2020). For both galaxies we obtain an rms background of ∼0.3 mJy beam⁻¹, consistent with the mean value ∼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, M 31, M 33, and the Small and Large Magellanic Clouds (SMC and LMC, respectively), 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 μm band (L_IR) is usually used as a proxy for the SFR. A power-law fit to L_1.4 GHz as a function of L_IR, $L_{1.4\mathrm{G}\mathrm{H}\mathrm{z}}=A{L}_{\mathrm{I}\mathrm{R}}^m$, with L_IR from K20, yields m = 1.00 ± 0.04, log A = 11.41 ± 0.37, and a dispersion of 0.2 dex, indicating a very tight relation. These values are in excellent agreement with 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_{1.4\mathrm{GHz}}=A{\dot{M}}_{\ast }^m$, yields a value of m = 1.13 ± 0.09 and log A = 21.14 ± 0.10, with a larger dispersion of 0.34 dex concentrated close to log Ṁ_∗ ∼ 0.

In the upper left panel of Fig. 1 we show the L_1.4 GHz–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 (Δ_Ṁ∗1.4 = 2.8), which 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_IR–SFR relationship from Kennicutt (1998), which is not a good proxy for galaxies with low SFRs (Bell 2003). This makes the L_1.4 GHz–SFR relationship presented by Yun et al. (2001) inconsistent with our fit in the low SFR range (we show both in the upper right panel of Fig. 1).

Fig. 1. Observed correlations for our galaxy sample. Upper left panel: L1.4 GHz as a function of the SFR. The blue solid line is the best fit to the data, and the shaded region its 68% confidence level. For comparison, the yellow dashed and green dotted lines show the Yun et al. (2001) and Bell (2003) relations, respectively. Upper right panel: L1.4 GHz as a function of LIR. The grey solid line is the best fit to the data, and the shaded region its 68% confidence level. For comparison, the green dotted line shows the Bell (2003) relation. Lower left panel: Lγ as a function of L1.4 GHz. The cyan solid line shows the best fit to the data, and the shaded region its 68% confidence level. For comparison, the brown dashed line shows the Ackermann et al. (2012) relation. Lower right panel: L150 MHz as a function of the SFR. The green solid line is the best fit to the data, and the shaded region its 68% confidence level. For comparison, the red dotted and violet dashed lines show the Gürkan et al. (2018) and Wang et al. (2019) relations, respectively. In all cases, the corresponding residuals weighted by their standard deviation are shown in the bottom panel; the black dotted lines correspond to 3σ.

This analysis shows that the luminosities we calculate for the restricted sample of galaxies detected in γ 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 γ-ray and radio luminosities (L_γ = A(L_1.4 GHz/10²¹ W Hz⁻¹)^m) gives m = 1.26 ± 0.10, log A = 39.04 ± 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 ± 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 due to the self-consistent and updated values of distances and fluxes we use. We show the resulting correlation in the lower left panel of Fig. 1, together with the Ackermann et al. (2012) relation.

2.3. The L_150 MHz–SFR

At low radio frequencies (≲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 (M 31, NGC 253, SMC, M 33, LMC, NGC 2146, NGC 2403, M 82, 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_150 MHz) using the K20 distances. We report these values in Table 1, together with the respective references. In the cases of M 31 and Arp 299 we assumed a typical flux error of 10%.

A power-law fit to the data, $L_{150\mathrm{MHz}}=A{\dot{M}}_{\ast }^m$, yields a value of m = 0.98 ± 0.13 and log A = 21.7 ± 0.15, with a dispersion of 0.5 dex. In the lower right panel of Fig. 1 we show the derived L_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 et al. 2021), reporting values 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). These authors state that the SFRs provided by their method differ by ∼0.2 dex from those obtained from Hα 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 (Ṁ_∗ ≈ 3 M_⊙ yr⁻¹). 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 ≈ 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 on the properties of non-thermal emission of SFGs. We keep this caveat in mind when interpreting results based on these data.

3. The model

We developed a non-thermal emission model to understand the observed L_1.4 GHz–SFR correlation (Sect. 2.2), while also contemplating the L_γ–SFR (Sect. 2.1) and L_150 MHz–SFR (Sect. 2.3) relations. The fact that two integrated properties of galaxies are involved in these correlations (the integrated non-thermal luminosity and the SFR) makes it reasonable to describe their emission based on the global CR energy balance. It is sensible to address the comprehension of the correlations between the non-thermal luminosity and the SFR in terms of a population of non-thermal particles that is injected at a rate that depends on the SFR of the galaxy. We further assume that the galactic environment in which these particles are injected is homogeneous and stationary. We self-consistently calculate the energy losses, diffusion, and advection of these relativistic particles and their non-thermal radiative output, corrected by the relevant absorption effects. In what follows we detail how we incorporate all this physics in a model that is valid over a broad range of SFRs.

3.1. CR injection and transport

We model the galaxy emission region as a disc of radius R and height 2H with constant gas density, and in which acceleration sites are uniformly distributed. We adopt R = 1 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

$$\begin{array}{c}\hfill Q_i(E)=Q_{0,i}{E}^{-\alpha }exp[-{(E/E_{\max ,i})}^{{\zeta }_i}],\end{array}$$(1)

where i = e, p (electrons and protons), α = 2.2 is the spectral index (Lacki & Thompson 2013), E_max is the maximum particle energy (assumed to be 10¹⁵ eV), and ζ_p, e = (1, 2) (Zirakashvili & Aharonian 2007; Blasi 2010). The normalisation constant Q_0, i is computed as

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{SN}}{E}_{\mathrm{SN}}{\xi }_{\mathrm{SN}}={\displaystyle {\int }_{E_{\min ,p}}^{E_{\max ,p}}}E{Q}_{0,\mathrm{p}}(E)\mathrm{d}E,\end{array}$$(2)

assuming that the supernova rate of the galaxy is ℛ_SN = (83 M_⊙)⁻¹ Ṁ_∗, consistent with the Chabrier (2003) initial mass function (IMF), and that a fraction ξ_CR = 10% of the supernova energy (E_SN = 10⁵¹ erg) is converted in CRs. We adopt E_min, p = 1.2 GeV as the minimum proton energy. Moreover, this energy is distributed as Q_0, p/Q_0, e = 50; this 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 of two) only affects the emission from primary electrons (also within a factor of two), which does not significantly impact 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 levels. Even though these sources can play an important role as acceleration sites, it is not clear yet to what extent they can contribute to the bulk of galactic CRs where, to date, 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 (e.g. young massive star clusters) is unlikely to change our main results because the acceleration efficiency values of the two components are 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 radiation, 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 the Thomson and the 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

$$\begin{array}{c}\hfill N_i(E)=\frac{1}{\mid \dot{E}\mid }{\displaystyle {\int }_E^{E_{\max ,i}}}{Q}_i(E)\mathrm{d}E,\end{array}$$(3)

where $\dot{E}=E{\tau }_{\mathrm{loss}}^{-1}(E)$ and ${\tau }_{\mathrm{loss}}=1/[{\tau }_{\mathrm{esc}}^{-1}+{\tau }_{\mathrm{cool}}^{-1}]$. Here τ_esc is the characteristic CR escape time computed as ${\tau }_{\mathrm{esc}}^{-1}={\tau }_{\mathrm{adv}}^{-1}+{\tau }_{\mathrm{diff}}^{-1}$ and ${\tau }_{\mathrm{cool}}^{-1}={\displaystyle \sum _j}{\tau }_j^{-1}$, where each τ_j represents the CR cooling time by each different cooling mechanism.

We note that this equation is almost equivalent to the standard leaky box solution, N_i(E)=Q_i(E) τ_loss(E), with the benefit that it yields a more precise normalisation for injection indices α >  2 when cooling dominates over escape.

3.2. Secondary products

Here we present the formulae used to calculate γ rays and secondary electron spectra produced in p − p interactions. Following Kelner et al. (2006), the emissivity of these secondary products i = γ, 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_i >  100 GeV the emissivity formula is

$$\begin{array}{c}\hfill Q_i(E_i)=cn{\displaystyle {\int }_{x_{i,\min }}^1}{\sigma }_{\mathrm{pp}}\left(\frac{E_i}{x_i}\right)N_{\mathrm{p}}\left(\frac{E_i}{x_i}\right){\stackrel{\sim }{F}}_i(x_i,\frac{E_i}{x_i})\frac{\mathrm{d}{x}_i}{x_i},\end{array}$$(4)

with σ_pp the cross section for inelastic p − p scattering

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{pp}}(E)=(34.3+1.88L+0.25L^2)\times {[1-{\left(\frac{E_{\mathrm{th}}}{E}\right)}^4]}^2\mathrm{mb},\end{array}$$(5)

where E_th = 1.22 GeV is the threshold energy for π₀ production, L = ln(E/TeV), and κ = 0.5 is the process inelasticity; ${\stackrel{\sim }{F}}_{\gamma }$ and ${\stackrel{\sim }{F}}_{\mathrm{e}}$ are given in Eqs. (58)–(61) and (62)–(65) of Kelner et al. (2006), respectively.

For energies E_i ≤ 100 GeV we describe the emissivity by the δ-functional approximation (Aharonian & Atoyan 2000):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Q_i(E_i)=& cn{\displaystyle {\int }_{E_{i,\min }}^{\infty }}\frac{2}{K_{\pi }}{\sigma }_{\mathrm{pp}}(\frac{E_{\pi }}{K_{\pi }}+m_{\mathrm{p}}{c}^2)N_{\mathrm{p}}(\frac{E_{\pi }}{K_{\pi }}+m_{\mathrm{p}}{c}^2)\hfill \\ \hfill & \times \frac{\stackrel{\sim }{f_i}(E_i/E_{\pi })}{\sqrt{E_{\pi }^2-m_{\pi }{c}^2}}\mathrm{d}{E}_{\pi }.\hfill \end{array}\end{array}$$(6)

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

3.3. Free-free absorption and emission

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 ρ, defined here as the ratio of ionised to total number densities. Thus, we compute the absorbed free-free thermal luminosity at a frequency ν, L_ff(ν), following the Olnon (1975) formalism for a homogeneous cylinder seen face-on, valid in the Rayleigh-Jeans regime (hν ≪ kT_e, where h and k are the Planck and Boltzmann constants, respectively)

$$\begin{array}{c}\hfill L_{\mathrm{ff}}(\nu )=\frac{8{\pi }^{2}k{T}_{\mathrm{e}}{\nu }^2}{c^2}{R}^2(1-e^{-{\tau }_{\mathrm{ff}}(\nu )}),\end{array}$$(7)

where τ_ff = α_ff(ν) l is the optical depth, l = 2H is the length of the ionised region along the line of sight, and α_ff is the absorption coefficient for photons of frequency ν. This coefficient is given by Rybicki & Lightman (1986)

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{ff}}=3.7\times {10}^8{\rho }^2{n}^2{T}_{\mathrm{e}}^{-1/2}{Z}^2{\nu }^{-3}[1-exp(-\frac{h\nu }{k{T}_{\mathrm{e}}})]{\overline{g}}_{\mathrm{ff}},\end{array}$$(8)

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

$$\begin{array}{c}\hfill {\overline{g}}_{\mathrm{ff}}=9.77(1+0.13\log \frac{T_{\mathrm{e}}^{3/2}}{Z\nu })·\end{array}$$(9)

Following Tielens (2005), we assume a constant ionisation fraction ρ = 0.05 throughout the whole SFR range, and we explore the impact of higher values, up to 0.1, on the luminosity–SFR correlations. In SBGs and ULIRGs the ISM conditions are less well-known. Despite the high star formation activity and supernova rate, most of the ionised components may be removed from the central galactic region and load galactic winds. This would lead to extremely low values of the ionisation ratio down to 10⁻⁴ (Krumholz et al. 2020). We take this into account by exploring scenarios where ρ is assumed to have such a low value.

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

In a similar manner, we obtain the outgoing non-thermal flux of the whole emission region by solving the corresponding radiative transport equation (see Rybicki & Lightman 1986). Modelling the emission region as a homogeneous face-on cylinder, the absorption-corrected synchrotron luminosity is

$$\begin{array}{c}\hfill L_{\mathrm{abs}}(\nu )=L_0(\nu )\left(\frac{1-e^{-{\tau }_{\mathrm{ff}}(\nu )}}{{\tau }_{\mathrm{ff}}(\nu )}\right),\end{array}$$(10)

where L₀ is the unabsorbed luminosity produced by non-thermal radiative processes in our model.

4. Modelling and scale relations

The SFR is the only independent variable in our model. We assume that the energy injected into protons and electrons depends linearly on the SFR (Eq. (2)). The observed correlations presented in Sect. 2 differ from linearity, which suggests that there is also a dependence of the main galactic properties on the SFR. In this section we are interested in finding simple parametrisations for the density n, the IR radiation field, and the magnetic and velocity fields with the SFR, motivated by physical and observational arguments.

4.1. Densities and L_IR

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

$$\begin{array}{c}\hfill {\dot{M}}_{\ast }[M_{\odot }{\mathrm{yr}}^{-1}]=1.72\times {10}^{-10}\frac{L_{\mathrm{IR}}}{L_{\odot }}(1+\sqrt{\frac{{10}^9}{L_{\mathrm{IR}}/L_{\odot }}})ϵ,\end{array}$$(11)

with ϵ = 0.79 (Crain et al. 2010) and L_IR the total infrared luminosity between 8 and 1000 μ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 Σ = L_IR/(8πR²σ_SBT⁴), where σ_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 (Σ_SFR) and cold gas mass (Σ_gas) of galaxies. Assuming that the emitter is a cylinder of radius R and thickness 2H (K20), we get

$$\begin{array}{c}\hfill n=\frac{{\mathrm{\Sigma }}_{\mathrm{gas}}}{2H{m}_{\mathrm{H}}}\propto {\dot{M}}_{\ast }^{0.71}{R}^{-1.42}{H}^{-1}.\end{array}$$(12)

4.2. Magnetic field

Robishaw et al. (2008) showed that strong magnetic fields (B of the order of milligauss, mG, hundreds of 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

$$\begin{array}{c}\hfill B=B_0{\left(\frac{{\dot{M}}_{\ast }}{10M_{\odot }{\mathrm{yr}}^{-1}}\right)}^{\beta },\end{array}$$(13)

where the normalisation factor B₀ is a free parameter and for the index β we explore three possible scenarios, β = 0.3, 0.5, 0.7, based on the equipartition arguments described below.

The weakest scaling is obtained relying on a stationary condition of magnetohydrodynamic (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 {\mathrm{\Sigma }}_{\mathrm{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 (${\mathrm{\Sigma }}_{\mathrm{SFR}}\propto {\mathrm{\Sigma }}_{\mathrm{g}}^{1.41}$) one can retrieve the index β ≈ 0.3. The second scenario is based on an equipartition condition between magnetic field and starlight, namely U_B ≈ U_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 β ≈ 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{\mathrm{\Sigma }}_{\mathrm{g}}^2$) and the magnetic field pressure (B²/8π) provides a natural limit for the magnetic field strength; higher values would make the galactic disc unstable. This condition, once the Kennicutt scaling is applied, leads to β = 0.7.

We assumed β = 0.3 as the fiducial value for the magnetic field index. In addition, we explored other prescriptions with a stronger B scaling (β = 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₀ = 100 μG as a fiducial normalisation, which yields B-values in agreement with the expected values for normal starbursts galaxies (see also Peretti et al. 2019) and with that inferred for the Galactic central region (see Ferrière 2009; Lacki & Thompson 2013). We also explored different values of B₀ across an order of magnitude to cover the different inferred values for the average magnetic field in SFGs, ranging from ∼μG in low SFR galaxies (Ferrière 2010; Jurusik et al. 2014) to ∼mG for ULIRGs (Robishaw et al. 2008).

4.3. Advection

Winds are ubiquitous in SFGs and are possibly one of the most spectacular consequences of their star formation activity. They are capable of driving mass out of galaxies at a rate of several M_⊙ yr⁻¹ with velocities up to 10³ km s⁻¹ (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 mainly responsible for wind launching in galaxies with mild SFRs; 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

$$\begin{array}{c}\hfill {\tau }_{\mathrm{adv}}({\dot{M}}_{\ast })=\frac{H}{v_{\mathrm{adv}}({\dot{M}}_{\ast })},\end{array}$$(14)

where v_adv is the SFR-dependent advection speed.

We divide the whole galaxy population into two classes: high SFR objects and mild SFR objects. We assume that objects with a star formation higher than the MW, Ṁ_∗ ≳ 2 M_⊙ yr⁻¹, are characterised by a starburst-driven wind, v_sw = 400 km s⁻¹ (see also K20, Yoast-Hull et al. 2013; Peretti et al. 2019). For objects with Ṁ_∗ < 2 M_⊙ yr⁻¹, 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, ${\mathit{v}}_{\mathrm{A}}=B/\sqrt{4\pi m_{\mathrm{p}}{n}_i}$. Putting all this together, we define the combined advection velocity over the entire SFR range as

$$\begin{array}{c}\hfill v_{\mathrm{b}}(\dot{M})=\{\begin{array}{cc}{v}_{\mathrm{A}},\hfill & {\dot{M}}_{\ast }<2M_{\odot }{\mathrm{yr}}^{-1}\hfill \\ v_{\mathrm{sw}},\hfill & {\dot{M}}_{\ast }>2M_{\odot }{\mathrm{yr}}^{-1}\hfill \end{array},\end{array}$$(15)

and we assume v_adv = v_b(Ṁ) to compute the Eq. (14). In addition, we explore alternative scenarios in which the nature of the advection speed is unique (either v_A or v_sw) for the entire SFR range considered.

4.4. Diffusion

Diffusion is a fundamental aspect of CR physics. At the microscopic level this phenomenon consists of 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 the 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

$$\begin{array}{c}\hfill D(E,{\dot{M}}_{\ast })\approx \frac{1}{3}{r}_{\mathrm{L}}v(E){\left(\frac{l_{\mathrm{c}}}{r_{\mathrm{L}}}\right)}^{1-\delta },\end{array}$$(16)

where v(E) is the particle speed, r_L is the Larmor radius, l_c represents the coherence length of the magnetic field, and δ characterises the nature of the turbulence. In particular, we have δ = 1/3 for the Kolmogorov turbulence, and δ = 1, 0.5 for Bohm and Kraichnan turbulence. In this work we assume Kolmogorov turbulence, δ = 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_{\mathrm{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 τ_SNR = 3 × 10⁴ yr (Truelove & McKee 1999). We also impose that such a parameter cannot be smaller than the typical dimension of a middle-aged SNR, namely l_c ≳ 10 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, τ_D, in which particles diffuse away from the galaxy

$$\begin{array}{c}\hfill {\tau }_{\mathrm{D}}(E,{\dot{M}}_{\ast })=\frac{H^2}{D(E,{\dot{M}}_{\ast })}·\end{array}$$(17)

We give the explicit dependency of τ_D on the SFR in Appendix B.

5. Results and discussion

To investigate how star formation influences the transport of CRs in SFGs, we set up a benchmark parameter configuration (scenario (0) in Table 2), to compute the global properties of galaxies throughout the whole SFR range (10⁻² − 10³ M_⊙ yr⁻¹), according to the scaling relations discussed in Sect. 4. We compare the timescales of the different processes involved, and solve the transport equation (Eq. (3)) following the prescriptions of Sects. 3 and 4, to determine the main characteristics of CR transport in different ranges of SFR. We present the results for this scenario in Sect. 5.1.

The transport of CRs shapes the non-thermal SED, and in turn each luminosity–SFR relation. We focus especially on the impact of different CR transport conditions on these relations. We therefore explore the outcome of our model in terms of the photon SED, from which we extract the integrated luminosity in the Fermi-LAT energy band (0.1−100 GeV) and the monochromatic luminosities at 1.4 GHz and 150 MHz.

We define an upper bound to the luminosity as a function of the SFR (calorimetric limit) in our model, to assess the ability of SFGs to cool CRs, converting their energy into radiation. This limit is built under the following conditions: (i) the escape is artificially turned off so that all injected particles cool inside the source; (ii) the non-thermal emission is not absorbed. This leads to maximum values $L_{\gamma }^{\mathrm{m}\mathrm{a}\mathrm{x}}$, $L_{1.4\mathrm{GHz}}^{\max }$, and $L_{150\mathrm{MHz}}^{\max }$ for each SFR.

We finally analyse the robustness of the model by exploring its parameter space, in order to assess the impact of different assumptions on the luminosity–SFR correlations. The scenarios considered are listed in Table 2. Each scenario from (1) to (8) represents a single parameter variation with respect to the benchmark scenario. We present, discuss, and compare our results for the three luminosity–SFR correlations in Sects. 5.2, 5.3, and 5.5.

5.1. Transport properties of the benchmark scenario

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 M_⊙ yr⁻¹ and 10 M_⊙ yr⁻¹.

Fig. 2. Cooling and escape times for primary electrons (upper panels) and protons (bottom panels) as a function of energy in scenario (0). Left panels: typical case of a low SFR galaxy (Ṁ∗ = 0.1 M⊙ yr−1). Right panels: typical case of a high SFR system (Ṁ∗ = 10 M⊙ yr−1). Solid lines represent different cooling processes: synchrotron (pink), ionisation (yellow), Bremsstrahlung (red), IC (green), and p − p (blue). The black dashed lines are for diffusion, whereas the grey double-dot-dashed lines are for advection.

Electrons. Ionisation losses dominate the transport of electrons with E_e <  0.1 GeV throughout the whole SFR range, while at higher electron energies the dominant loss mechanism depends on the SFR.

In the low SFR range, Ṁ_∗ ≲ 1 M_⊙ yr⁻¹, electrons with energies between 0.5 GeV and 10 GeV escape efficiently from galaxies, whereas synchrotron confines electrons at higher energies. As the SFR increases, Bremsstrahlung losses compete with the escape at E_e ∼ 1 GeV, while synchrotron keeps dominating the high energy domain together with IC scattering. In particular, the latter increases its relative importance for increasing SFR and eventually becomes the dominant cooling mechanism for SFR ≳ 10² M_⊙ yr⁻¹. Electron calorimetry is achieved for Ṁ_∗ ≳ 5 M_⊙ yr⁻¹, at which Bremsstrahlung dominates the cooling of GeV electrons; higher energy electrons are cooling-dominated for all SFR ranges.

Protons. The transport of protons is regulated by the competition between escape (diffusion and advection) and energy losses (p − p collisions). At low SFRs, Ṁ_∗ ≲ 1 M_⊙ yr⁻¹, most of the protons escape via diffusion for E_p ≳ 10 GeV and via advection at lower energies. At higher SFRs, up to a few tens of M_⊙ yr⁻¹, p − p losses become the dominant loss mechanism. However calorimetry is not completely achieved because the advection and p − p timescales remain comparable. A transition to a diffusive-dominated transport also occurs in the multi-TeV range. Finally, we note that galaxies behave like good proton calorimeters for SFRs ≳ 30 M_⊙ yr⁻¹.

5.2. The L_γ–SFR relation

We define L_γ as the integrated luminosity in the 0.1−100 GeV energy band, and we show in Fig. 3 the L_γ–SFR relation obtained for the benchmark scenario. The main contribution to L_γ comes from the p − p interaction, while γ rays of leptonic origin are subdominant (see Fig. B.1). In agreement with K20, the shape of the L_γ–SFR relation at low SFRs is influenced by the escape of most of the injected protons, whereas at high SFRs 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_adv = v_A or v_adv = v_sw. Since v_sw >  v_A (v_A ≳ 50 km s⁻¹ for the fiducial scenario), the CR escape is more efficient when v_adv = v_sw and consequently the γ-ray luminosity is reduced. Our benchmark case is in good agreement with the best fit of the observed sample.

Fig. 3. Lγ–SFR relation for scenario (0) using three different prescriptions for vadv: va (solid blue line), vsw (solid pink line), and vb from Eq. (15) (violet dashed line). The grey line and shaded band indicate the best fit to the data (purple circles) and its 1σ confidence region, respectively. The yellow region is forbidden by the calorimetric limit.

Variations of β, B₀, and T_e do not affect the emission of protons and, consequently, the L_γ–SFR correlation is not altered. On the contrary, the variation of the ionisation fraction ρ (scenarios (5) and (6)) modifies the CR escape through v_A at low SFRs, affecting the total source luminosity.

As described in Sect. 3.3, we explore values for the ionisation fraction down to 10⁻⁴. While these values can characterise powerful starbursts such as ULIRGs, they are unlikely to be physical in standard SFGs with low SFRs. Such a small ionisation fraction leads to an Alfvén speed larger than the average velocity of starburst superwinds resulting in an advection dominated transport at odds with observations of CRs (Blasi 2013) and γ rays (Ajello et al. 2020). The scenario ρ = 0.1 does not show any relevant modifications with respect to the benchmark scenario since the dominant escape mechanism is diffusion.

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

The L_γ–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 Ṁ_∗ ≥ 30 M_⊙ yr⁻¹, 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.

5.3. The L_1.4 GHz–SFR relation

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 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. The thermal contribution eventually becomes comparable to that of synchrotron at Ṁ_∗ > 100 M_⊙ yr⁻¹ (see Fig. B.1).

Fig. 4. Results for scenario (0) using vadv = vb from Eq. (15) (violet dashed line). Left panel: L1.4 GHz–SFR relation. The green dashed and dot-dashed lines represent the primary and secondary synchrotron electron contribution, respectively, and the dotted violet line the free-free contribution. The total emission at L1.4 GHz is also shown for two different prescriptions for vadv, vA (solid blue line) and vsw (solid pink line). Right panel: Lγ–L1.4 GHz relation. In both panels the grey shaded band indicates the 1σ confidence region of each fit (solid grey line) to the data sample (purple points).

The emission of primary electrons dominates the L_1.4 GHz–SFR correlation in the low SFR range (≲1 M_⊙ yr⁻¹), whereas synchrotron from secondaries gradually increases with the SFR and overcomes that of primaries at high SFRs (see respectively green dot-dashed and dashed lines in left panel of Fig. 4). This 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 thus 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 ≳30 M_⊙ yr⁻¹, 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 SFRs ≲ 10² M_⊙ yr⁻¹, L_1.4 GHz is dominated by the synchrotron emission, hereafter $L_{\mathrm{syn}}^{1.4}$. In order to understand the shape of the correlation, we note that $L_{\mathrm{syn}}^{1.4}$ is dominated by the emission of ∼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 α − 1 when ionisation dominates to α + 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 Sect. 3.2). This leads to a complex dependence of $L_{\mathrm{syn}}^{1.4}$ on the SFR, which we parametrise as

$$\begin{array}{c}\hfill L_{\mathrm{syn}}^{1.4}({\dot{M}}_{\ast })\propto {\dot{M}}_{\ast }^{\eta ({\dot{M}}_{\ast })},\end{array}$$(18)

where η is the slope of the synchrotron luminosity in the logarithmic L_1.4 GHz–SFR correlation. We note that this component dominates the L_1.4 GHz–SFR correlation up to ∼10² M_⊙ yr⁻¹ for the benchmark scenario, where we constrain η to the range 0.8−1.8. In particular, the slope of the correlation is steeper in the low SFR range, where the relative contribution of secondaries increases with the SFR. At high SFRs, the transport of protons also becomes calorimetric and the slope of the correlation softens (see Appendix B for a detailed discussion). We note that in calorimetric conditions η <  1. This results from the higher efficiency of ionisation and Bremsstrahlung in cooling GeV electrons.

The relevance of thermal processes increases for higher SFRs as a consequence of the higher density of the ISM and our assumption of constant ionisation fraction ρ throughout the whole SFR range. In particular, the thermal contribution to L_1.4 GHz is at the level of ∼5% for low SFR galaxies; this is consistent with values reported by Condon (1992) for the thermal emission in the MW and M 82. At high SFRs (≳100 M_⊙ yr⁻¹) the thermal contribution is ∼50%, and the free-free absorption also becomes significant. Nonetheless, assumptions in the ionisation fraction ρ impact on the relative importance of the free-free processes, as we 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 GHz–SFR correlation. In the right panel of Fig. 4 we present observational data in the L_γ–L_1.4 GHz plane and compare them with the outcomes of our model; we obtain a good agreement throughout the whole SFR range.

5.4. Parameter dependence of the L_1.4 GHz–SFR relation

In this section we discuss the impact of parameter variations on the L_1.4 GHz–SFR correlation.

Magnetic field index. As the thermal emission is unaltered by changes in β, here we focus only on the synchrotron component. The top panels in Fig. 5 illustrate how the variation of β (scenarios (1) and (2)) affects the slope of the L_1.4 GHz–SFR and L_1.4 GHz–L_IR correlations.

Fig. 5. Model relations for the different scenarios considered (Table 2). The grey shaded band indicates the 1σ confidence region of each fit (grey solid line) to the data sample (purple points). In all panels the dotted lines represent the thermal emission, the dashed lines the non-thermal emission, and the solid line the total emission. Upper panel: for scenarios (0), (1), and (2) the L1.4 GHz–SFR relation (left) and the L1.4 GHz–LIR relation (right) are shown. Middle panel: for scenarios (0), (3), and (4) (left) and (0), (5), and (6) (right) the L1.4 GHz–SFR relation is shown. Bottom panel: for scenarios (0), (7), and (8) the L1.4 GHz–SFR relation (left) is shown and for scenarios (0), (5), and (6) the L150 MHz–SFR relation (right) is shown.

In particular, the larger the value of β, the steeper the slope η of the L_1.4 GHz–SFR correlation. This is a natural consequence of the dependence of the synchrotron luminosity on the magnetic field which, for larger β, presents a larger variation throughout the SFR range (see Appendix B).

At the same time, lower values of B resulting from β >  0.3 at low SFRs reduce the escape time, which in turn reduces the luminosity and makes the correlation steeper. The benchmark scenario with β = 0.3 provides the best fit to the general trend of the observed data. However, the assumption β = 0.5, in agreement with Lacki et al. (2010), is slightly more appropriate while fitting data at high SFR. The scenario assuming β = 0.7 is disfavoured when compared with data in the entire SFR range, whereas, as in the case of β = 0.5, it can be adopted to describe the high SFR range.

Magnetic field normalisation. A variation of B₀ (scenarios (3) and (4)) results in a different normalisation of the correlation curve of about a factor ∼1.5 − 4. We show these results in the middle left panel of Fig. 5. This behaviour is in agreement with the dependence of the synchrotron luminosity of a non-thermal electron population with the magnetic field (see e.g. Ghisellini 2013). However, the synchrotron radiation at 1.4 GHz is mostly produced by electrons in the ∼GeV range. At these energies several physical processes compete in shaping the electron energy distribution. This leads to an energy distribution different from a simple power law. Thus, in the context of SFGs, the dependence of L_1.4 GHz on the magnetic field strength is not trivial, and a change in the normalisation constant B₀ introduces an overall shift which slightly depends on the SFR. Despite this complexity, we show that the parameter B₀ does not affect the global slope of the correlation, but only influences its normalisation (see Appendix B for a more detailed and quantitative discussion). The value B₀ = 100 is the one that best fits the data. The resulting correlations for B₀ = 50 and 500 limit the data above and below, respectively, excluding only NGC 4945.

Ionisation fraction. The impact of different assumptions for the ionisation fraction ρ on the L_1.4 GHz–SFR correlation is threefold since it affects: (1) the free-free emission, (2) the free-free absorption, and (3) the advection timescale in low SFGs. The middle right panel of Fig. 5 illustrates the results obtained in the benchmark case (blue lines) compared with ρ = 10⁻⁴ (violet lines) and ρ = 10⁻¹ (orange lines).

Although a value of ρ = 0.1 produces results consistent with the observed trend, at high SFR the thermal emission becomes significant and possibly dominates. This is the opposite of what is expected for galaxies like M 82 where the estimated thermal fraction is even lower (∼2%, Basu et al. 2012) than for normal galaxies. The scenario where ρ = 0.1 is therefore unlikely to be realised in nature and at odds with multiwavelength modelling of ULIRGs (see e.g. Yoast-Hull et al. 2017). On the other hand, lower values of ρ reduce the absorption and the relative importance of thermal emission at high SFR. This yields a thermal contribution that is too low with respect to what is generally accepted to hold for normal galaxies (Condon 1992). Finally, ρ = 10⁻⁴ underestimates the data at low SFRs, and the model behaviour is similar to that of γ rays.

We conclude that the scenario where ρ = 0.05 is the one that better accommodates the observed data and remains consistent with synchrotron dominating the emission at 1.4 GHz. Although it is tempting to speculate on a possible SFR dependence of ρ, this would introduce further complexity in the model, and we therefore leave it for upcoming investigations.

Electron temperature. T_e ∼ 10⁴ K is the standard temperature assumed for the warm ionised ISM. In starburst regions, this gas coexists with a hot and pressurised component resulting from the collective effect of supernovae and stellar winds. This implies that T_e >  10⁴ K might be consistent with the average temperature in these environments (Westmoquette et al. 2009). The variation of T_e (scenarios (7) and (8)) affects the thermal emission and absorption. This behaviour is shown in Fig. 5 (bottom left panel) from which it can be concluded that values between T_e = 10⁴ and 10⁵ K reproduce well the observations, whereas lower values of T_e lead to an excessive absorption at high SFR compared to the data trend.

These results can be understood by looking at the physical dependence of the opacity on the temperature. In particular, at a given photon frequency, ${\tau }_{\mathrm{ff}}\propto T_{\mathrm{e}}^{-3/2}{\dot{M}}_{\ast }^{1.7}$, and so the optically thick regime is reached at higher SFRs for higher values of T_e. Thus, the absorption of synchrotron radiation decreases with T_e.

Finally, we observe that the different assumptions about the temperature do not strongly impact the relative contribution of thermal emission to the L_1.4 GHz across the whole SFR range.

5.5. The L_150 MHz–SFR relation

The physical processes governing the emission at 150 MHz are the same as at 1.4 GHz (Sect. 5.3). However, at this frequency free-free absorption plays a crucial role. The turnover frequency where τ_ff >  1 increases up to 150 MHz for Ṁ_∗ ∼ 10 M_⊙ yr⁻¹. Therefore, the 150 MHz synchrotron emission is completely absorbed at high SFRs in our benchmark scenario. This is reflected in the L_150 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, ρ ∼ 5%. In addition, the calorimetric limit determined for the benchmark scenario also shows that it can hardly match the observed fluxes for SFRs ≳ 10² M_⊙ yr⁻¹.

Fig. 6. L150 MHz–SFR relation obtained for scenario (0) using three different prescriptions for vadv: vA (solid cyan line), vsw (solid pink line), and vb from Eq. (15) (violet dashed line). The green dashed and green dot-dashed lines represent the primary and secondary electron contributions, respectively, and the dotted violet line represents the free-free contribution. The grey shaded band indicates the 1σ confidence region of the fit (solid grey line) to the data sample (purple points).

Although the effects of the variation of parameters in the L_150 MHz–SFR correlation are similar to those in the L_1.4 GHz–SFR correlation, some details are required to determine if it is possible to match the trend of data with a model prediction under different parametric assumptions.

The assumption β = 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 ρ and T_e can also mitigate or enhance the absorption as evidenced in the bottom right 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 ρ and/or T_e would provide a better fit at a 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.

6. Conclusions

We investigated how the transport of CRs, produced by sources related to the star formation, shapes the observed L_γ–SFR, L_1.4 GHz–SFR, and L_150 MHz–SFR correlations for SFGs detected in γ rays. We presented a flexible emission model for SFGs over a broad range of SFRs and compared the theoretical estimates with the most updated observational data available. The main results from our work are summarised as follows:

• We provide homogeneous radio continuum fluxes and their errors computed directly from images at 1.4 GHz of Stokes-I available in the NRAO/VLASky continuum Survey catalogue for NGC 253, NGC 1068, NGC 2146, NGC 2403, M 82, Arp 220, Arp 299, and NGC 3424. We also computed for the first time the fluxes and the errors for NGC 4945 and Circinus from the new 1.4 GHz Stokes-I continuum image provided by The ASKAP Continuum Survey.

• We provide the first self-consistent data set at 150 MHz and 1.4 GHz for all galaxies detected in γ rays. We present updated correlations of L_γ–SFR, L_1.4 GHz–SFR, and L_150 MHz–SFR for these galaxies.

• The trend followed by the L_1.4 GHz–SFR correlation observed in the subsample of SFGs detected at γ rays agrees very well with that obtained from larger surveys of SFGs (Yun et al. 2001; Bell 2003).

• We revisited the L_γ–L_1.4 GHz correlation and found a steeper slope m = 1.26 ± 0.10 (with $L_{\gamma }\propto L_{1.4\mathrm{G}\mathrm{H}\mathrm{z}}^m$) with respect to the value reported by Ackermann et al. (2012).

• We obtained a slope of m = 0.98 ± 0.13 for the L_150 MHz–SFR correlation (with $L_{150\mathrm{MHz}}\propto {\dot{M}}_{\ast }^m$), which is marginally shallower than results reported by other authors (ranging from 1 to 1.3).

• We developed a self-consistent model for CR transport and multiwavelength emission suitable to interpret the luminosity–SFR correlation at different energy ranges over a very broad range of SFRs. This model can reproduce reasonably well the observed L_γ–SFR, L_1.4 GHz–SFR, and L_1.4 GHz–L_γ correlations.

• A constant ionisation fraction ρ fails to reproduce the L_150 MHz–SFR correlation throughout the whole SFR range. In particular, an ionisation fraction of ∼5% at low SFRs is compatible with observations, whereas at high SFRs a value < 1% is favoured to achieve the emission observed in ULIRGs. Small values of ρ would be in agreement with the results of Valdés et al. (2005) who suggest that a significant fraction of the ionising photons may be absorbed by dust. Moreover, Krumholz et al. (2020) propose that typical ionisation fractions in starbursts should be ∼10⁻⁴ since the hot gas would be advected in the wind.

• A scaling law B = B₀(Ṁ_∗/10 M_⊙ yr⁻¹)^β is compatible with the observed data for β = 0.3, which suggests that the magnetic field strength is regulated by stationary MHD turbulence. Higher beta values are also compatible with data limited to the range of high SFR. The normalisation value that best describes the magnetic field at 10 M_⊙ yr⁻¹ is of the order of 100 μG.

• The coexistence of synchrotron components from both primary and secondary electrons is mandatory to reproduce the slope of the L_1.4 GHz–SFR correlation. Such a slope reflects the transport conditions of primary protons and could provide indications of a transition to a calorimetric regime.

• High temperatures of 10⁴ − 10⁵ K for the ionised gas are more suitable to explain the emission of the ULIRGs.

• Galaxies behave as good proton calorimeters for SFR ≳ 30 M_⊙ yr⁻¹. Such a calorimetry can be assessed via observations in γ rays. On the contrary, relativistic electrons, although they reach the calorimetric limit for SFR ≳ 5 M_⊙ yr⁻¹, are hard to assess given that their emission is not completely dominant for all SFRs in either the radio or the γ-ray band.

• Different wind prescriptions depending on the SFR of the galaxy yield a good description of the L_γ–SFR correlation. In particular, an Alfvenic wind is proposed to dominate in low-SFR galaxies, whereas a thermally driven wind would better describe the outflow at high SFRs.

The SMC deserves a more thorough discussion. It is the galaxy at the lower end of the SFR range, separated by far from the next one, and therefore may play an important role in model fitting. Its non-thermal luminosity is almost always underpredicted by our models. It might be argued that this is the result of additional contributions from either unresolved point-like sources or from the halo. However, this seems very improbable, due to the short distance and extended nature of this galaxy. The discrepancy can be mitigated taking into account the uncertainty in the current SFR estimates for the SMC. We used the K20 value (0.027 M_⊙ yr⁻¹) obtained from an Hα/IR proxy, sensitive to the SFR in a time interval of ∼10 Myr, whereas it has been found that the SMC has a complex SFR history with recent bursts of up to ten times this value (Harris & Zaritsky 2004). An increase in the SFR of the SMC of 0.5−1 dex almost always results in a better agreement between observations and our L_1.4 GHz–SFR and L_150 MHz–SFR modelled correlations at low SFRs.

Enlarging the current sample of γ-ray emitting SFGs, and developing reliable techniques capable of distinguishing the emission associated with AGN hosted by some galaxies of the sample is necessary to obtain a more accurate L_γ–SFR correlation. A larger sample of galaxies detected in γ rays would also yield a more accurate L_γ − L_1.4 GHz correlation, which is important for two reasons. The first is that it could provide an independent method for the classification of γ-ray galaxies, in addition to the standard method based on the spectral index determination. The second reason is that it could be used to estimate the contribution of SFGs to the γ-ray background by extrapolating the radio luminosity function to high redshift.

A careful study of the L_150 MHz–SFR correlation is important to assess the extent to which L_150 MHz traces star formation. Furthermore, a more refined L_150 MHz–SFR correlation would allow us to discern if the SFGs detected in γ rays follow a distinguishable trend with respect to more general samples. Due to the difficulties in reproducing the emission at these low frequencies, in upcoming studies we will address the effect produced on the modelled correlation by an ionisation fraction diminishing with the SFR, or a highly clumped gas distribution resulting in a porous medium. An alternative scenario can be a non-negligible contribution from particles radiating from a region larger than our galactic box, for instance electrons diffusing or advected in the halo (Heesen et al. 2009), or the presence of an additional source or acceleration region (such as a co-existing AGN) not directly related to the SFR (Gürkan et al. 2018).

The next generation of telescopes with better angular resolution and sensitivity will allow us to unravel the morphology of the diffuse non-thermal emitting region in SFGs, possibly resolving the emission coming from the disc from that originated in the halo or wind (see e.g. Romero et al. 2018; Müller et al. 2020; Peretti et al. 2021). The panchromatic view of the SFR–luminosity correlations in SFGs is a key to assessing the physical processes that govern the non-thermal emission of these sources. In particular, we plan to extend the study of luminosity correlations to the VHE and hard X-ray energy bands. The former can provide a unique identification tool for possible new detections by the upcoming Cherenkov Telescope Array (CTA; Cherenkov Telescope Array Consortium 2019). The latter can provide a new powerful tool for investigating the CR transport properties.

Acknowledgments

P.K. and P.B. acknowledge B. Koribalski for help in providing the ASKAP data. The research activity of E.P. was supported by Villum Fonden (project n. 18994) and by the European Unions Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 847523 “INTERACTIONS”. P.B. acknowledges support from ANPCyT PICT 0773–2017. L.J.P, P.K. and P.B. acknowledge also support from CONICET grant PIP 2014-0265.

References

Appendix A: Radio images

Here we present the Stokes-I continuum images at 1.4 GHz for the ten galaxies from which we extract the flux. In Fig. A.1 we show the emission above 3σ and the selected contours for the flux integration of each galaxy as explained in Sect. 2.2.

Fig. A.1. Stokes-I continuum image at 1.4 GHz for each galaxy in our sample. The flux scale is in Jy beam−1. The magenta contours indicate the 3σ level and the green solid contour the region selected for the flux integration. The synthesised beam is shown with a blue circle at the bottom left corner of each image.

Appendix B: Model details

In this appendix we specify the dependencies of the escape and cooling times with the free parameters of the model (β, B₀, ρ) and the SFR. We also show in Fig. B.1 the modelled SEDs for the benchmark scenario at the two fixed SFRs presented in Sect. 5.1 (Ṁ_∗ = 0.1 and 10 M_⊙ yr⁻¹). Finally, we make a detailed derivation of the slope η presented in Sect. 5.1, using approximations in different ranges of SFRs.

B.1. Cooling times

Here we show the dependencies of the cooling times on the SFR and the free parameters of the model,

$$\begin{array}{c}\hfill {\tau }_{\mathrm{adv}}(\dot{M})\propto \frac{\sqrt{\rho n}}{B}\propto \frac{\sqrt{\rho }}{B_0{\dot{M}}_{\ast }^{\beta -0.35}},{\dot{M}}_{\ast }<2M_{\odot }{\mathrm{yr}}^{-1},\end{array}$$(B.1)

which is constant for Ṁ_∗ ≥ 2 M_⊙ yr⁻¹,

$$\begin{array}{cc}& {\tau }_D({\dot{M}}_{\ast })\propto B^{\delta }{\dot{M}}_{\ast }^{\frac{1-\delta }{3}}\propto B_0^{1/3}{\dot{M}}_{\ast }^{\beta /3+2/9},\hfill \end{array}$$(B.2)

$$\begin{array}{cc}& {\tau }_{\mathrm{sync}}({\dot{M}}_{\ast })\propto \frac{1}{B^2}\propto \frac{1}{B_0^2{\dot{M}}_{\ast }^{2\beta }},\hfill \end{array}$$(B.3)

$$\begin{array}{cc}& {\tau }_{\mathrm{IC}}({\dot{M}}_{\ast })\propto \frac{1}{L_{\mathrm{IR}}}\propto \frac{1}{{\dot{M}}_{\ast }},\hfill \end{array}$$(B.4)

while τ_brem, τ_ion, and τ_pp are all $\propto n^{-1}\propto {\dot{M}}_{\ast }^{-0.71}$.

B.2. Model SEDs

We show the modelled SEDs for Ṁ_∗ = 0.1 M_⊙ yr⁻¹ in the left panel of Fig. B.1. The dominant contribution at 150 MHz and 1.4 GHz comes from the synchrotron emission of primary electrons. In the case of Ṁ_∗ = 10 M_⊙ yr⁻¹, which we show in right panel of Fig. B.1, the dominant contribution comes instead from synchrotron emission of secondary electrons. In both cases the free-free thermal component becomes important at higher frequencies. In the Fermi energy range the emission is dominated by hadronic emission for both SFRs.

B.3. The slope of the L_1.4 GHz–SFR relation

In order to understand the behaviour of the slope η of Eq. (18) in the benchmark scenario, we approximate the electron distribution as a power law (N(E)=N₀ E^−p) in the electron-energy range producing the dominant synchrotron contribution to radiation at 1.4 GHz (ℰ_1.4). This ranges between ∼1 GeV and ∼5 GeV for the benchmark scenario. Under this assumption, the synchrotron luminosity at a given frequency can be estimated roughly as $L_{\mathrm{syn}}^{1.4}\propto N_0{B}^{(p+1)/2}$ (Ghisellini 2013). For simplicity we focus on two typical values at low and high SFRs, where we can assume that electrons lose their energy through a single dominant process near ℰ_1.4.

At low SFR $L_{\mathrm{syn}}^{1.4}$ is shaped by primary synchrotron emission, and near Ṁ_∗ = 0.1 M_⊙ yr⁻¹ the electron distribution is manly dominated by diffusion at ℰ_1.4 ∼ 4 GeV (see upper left panel in Fig. 2). Then p = α + δ and N_0,pri ∝ Ṁ_∗τ_D, which results in $L_{\mathrm{syn}}^{1.4}\propto {\dot{M}}_{\ast }^{1.8}$. At high SFRs the slope of the primary electrons subdominant contribution can be estimated in a similar way taking into account that Bremsstrahlung dominates cooling. Then p = α, N_0,pri ∝ Ṁ_∗τ_brem, and $L_{\mathrm{syn}}^{1.4}\propto {\dot{M}}_{\ast }^{0.8}$.

The emission of secondary electrons is dominant in $L_{\mathrm{syn}}^{1.4}$ at high SFR. The normalisation of their particle distribution depends on the proton distribution roughly as $N_{0,\sec }\propto n{Q}_{0,p}{\tau }_p{\tau }_e\propto {\dot{M}}_{\ast }{\dot{M}}_{\ast }^{0.7}{\tau }_p{\tau }_e$ (see Sect. 3.2), where τ_p and τ_e are the dominant loss times for protons and electrons, respectively.

Near Ṁ_∗ ∼ 10 M_⊙ yr⁻¹ proton-proton interactions and advection dominate the proton losses, and Bremsstrahlung controls the electron cooling near ℰ_1.4 (see lower right panel on Fig. 2). On the one hand, if we assume that proton-proton interactions dominate alone, then p = α, $N_{0,\sec }\propto {\dot{M}}_{\ast }^{1.7}{\tau }_{\mathrm{pp}}{\tau }_{\mathrm{brem}}$, and $L_{\mathrm{syn}}^{1.4}\propto {\dot{M}}_{\ast }^{0.8}$. On the other hand, if we assume that advection domains alone, then $N_{0,\sec }\propto {\dot{M}}_{\ast }^{1.7}{\tau }_{\mathrm{adv}}{\tau }_{\mathrm{brem}}$, where $L_{\mathrm{syn}}\propto {\dot{M}}_{\ast }^{1.5}$. The combined effect of these two losses returns a dependence in SFR between these last two values.

Fig. B.1. SEDs computed from scenario 0, for typical cases of low (Ṁ∗ = 0.1 M⊙/yr, left) and high (Ṁ∗ = 10 M⊙/yr, right) SFRs. The coloured lines (solid for the primary CRs and dashed for the secondary CRs) are the individual contributions of different non-thermal radiative mechanisms (the colour-coding is the same as in Fig. 2). The dashed brown line is the thermal free-free emission; the grey shaded region represents the Fermi-LAT energy range; and the vertical dotted black lines show the studied radio frequencies (150 MHz and 1.4 GHz).

Appendix C: Hadronic processes

For the sake of completeness we report here the functions we adopted in the computation of all hadronic byproducts,

$$\begin{array}{c}\hfill F_{\gamma }(x,E_p)=B_{\gamma }\frac{d}{\mathit{dx}}[\ln (x)(\frac{1-x^{{\beta }_{\gamma }}}{1+k_{\gamma }{x}^{{\beta }_{\gamma }}(1-x^{{\beta }_{\gamma }})}{)}^4]=\\ \hfill B_{\gamma }\frac{\ln (x)}{x}(\frac{1-x^{{\beta }_{\gamma }}}{1+k_{\gamma }{x}^{{\beta }_{\gamma }}(1-x^{{\beta }_{\gamma }})}{)}^4\times \\ \hfill [\frac{1}{\ln (x)}-\frac{4{\beta }_{\gamma }{x}^{{\beta }_{\gamma }}}{1-x^{{\beta }_{\gamma }}}-\frac{4k_{\gamma }{\beta }_{\gamma }{x}^{{\beta }_{\gamma }}(1-2x^{{\beta }_{\gamma }})}{1+k_{\gamma }{x}^{{\beta }_{\gamma }}(1-x^{{\beta }_{\gamma }})}],\end{array}$$(C.1)

where x = E_γ/E_p and

$$\begin{array}{c}\hfill B_{\gamma }=1.30+0.14L+0.011L^2,\end{array}$$(C.2)

$$\begin{array}{c}\hfill {\beta }_{\gamma }=\frac{1}{1.79+0.11L+0.008L^2},\end{array}$$(C.3)

$$\begin{array}{c}\hfill k_{\gamma }=\frac{1}{0.801+0.049L+0.014L^2},\end{array}$$(C.4)

and where L = ln(E_p/1 TeV);

$$\begin{array}{c}\hfill F_e(x,E_p)=B_e\frac{{[1+k_e{(\ln x)}^2]}^3}{x(1+0.3/x^{{\beta }_e})}{[-\ln (x)]}^5,\end{array}$$(C.5)

where

$$\begin{array}{c}\hfill B_e=\frac{1}{69.5+2.65L+0.3L^2},\end{array}$$(C.6)

$$\begin{array}{c}\hfill {\beta }_e=\frac{1}{{(0.201+0.062L+0.00042L^2)}^{1/4}},\end{array}$$(C.7)

$$\begin{array}{c}\hfill k_e=\frac{0.279+0.141L+0.0172L^2}{0.3+{(2.3+L^2)}^2}.\end{array}$$(C.8)

We also report the function ${\stackrel{\sim }{f}}_e$:

$$\begin{array}{c}\hfill {\stackrel{\sim }{f}}_e(x)=f_{{\nu }_{\mu }^{(2)}}(x)=g_{{\nu }_{\mu }}(x)\mathrm{\Theta }(x-r)+[h_{{\nu }_{\mu }}^{(1)}(x)+h_{{\nu }_{\mu }}^{(2)}(x)]\mathrm{\Theta }(r-x),\end{array}$$(C.9)

where r = 1 − λ = (m_μ/m_π)² = 0.573,

$$\begin{array}{c}\hfill g_{{\nu }_{\mu }}(x)=\frac{3-2r}{9{(1-r)}^2}(9x^2-6\ln x-4x^3-5),\end{array}$$(C.10)

$$\begin{array}{c}\hfill h_{{\nu }_{\mu }}^{(1)}(x)=\frac{3-2r}{9{(1-r)}^2}(9r^2-6\ln r-4r^3-5),\end{array}$$(C.11)

$$\begin{array}{c}\hfill h_{{\nu }_{\mu }}^{(2)}(x)=\frac{(1+2r)(r-x)}{9r^2}\times [9(r+x)-4(r^2+rx+x^2)],\end{array}$$(C.12)

and Θ is the Heaviside step function.

All Tables

All Figures

Fig. 1. Observed correlations for our galaxy sample. Upper left panel: L1.4 GHz as a function of the SFR. The blue solid line is the best fit to the data, and the shaded region its 68% confidence level. For comparison, the yellow dashed and green dotted lines show the Yun et al. (2001) and Bell (2003) relations, respectively. Upper right panel: L1.4 GHz as a function of LIR. The grey solid line is the best fit to the data, and the shaded region its 68% confidence level. For comparison, the green dotted line shows the Bell (2003) relation. Lower left panel: Lγ as a function of L1.4 GHz. The cyan solid line shows the best fit to the data, and the shaded region its 68% confidence level. For comparison, the brown dashed line shows the Ackermann et al. (2012) relation. Lower right panel: L150 MHz as a function of the SFR. The green solid line is the best fit to the data, and the shaded region its 68% confidence level. For comparison, the red dotted and violet dashed lines show the Gürkan et al. (2018) and Wang et al. (2019) relations, respectively. In all cases, the corresponding residuals weighted by their standard deviation are shown in the bottom panel; the black dotted lines correspond to 3σ.

In the text

Fig. 2. Cooling and escape times for primary electrons (upper panels) and protons (bottom panels) as a function of energy in scenario (0). Left panels: typical case of a low SFR galaxy (Ṁ∗ = 0.1 M⊙ yr−1). Right panels: typical case of a high SFR system (Ṁ∗ = 10 M⊙ yr−1). Solid lines represent different cooling processes: synchrotron (pink), ionisation (yellow), Bremsstrahlung (red), IC (green), and p − p (blue). The black dashed lines are for diffusion, whereas the grey double-dot-dashed lines are for advection.

In the text

	Fig. 3. Lγ–SFR relation for scenario (0) using three different prescriptions for vadv: va (solid blue line), vsw (solid pink line), and vb from Eq. (15) (violet dashed line). The grey line and shaded band indicate the best fit to the data (purple circles) and its 1σ confidence region, respectively. The yellow region is forbidden by the calorimetric limit.
In the text

Fig. 4. Results for scenario (0) using vadv = vb from Eq. (15) (violet dashed line). Left panel: L1.4 GHz–SFR relation. The green dashed and dot-dashed lines represent the primary and secondary synchrotron electron contribution, respectively, and the dotted violet line the free-free contribution. The total emission at L1.4 GHz is also shown for two different prescriptions for vadv, vA (solid blue line) and vsw (solid pink line). Right panel: Lγ–L1.4 GHz relation. In both panels the grey shaded band indicates the 1σ confidence region of each fit (solid grey line) to the data sample (purple points).

In the text

Fig. 5. Model relations for the different scenarios considered (Table 2). The grey shaded band indicates the 1σ confidence region of each fit (grey solid line) to the data sample (purple points). In all panels the dotted lines represent the thermal emission, the dashed lines the non-thermal emission, and the solid line the total emission. Upper panel: for scenarios (0), (1), and (2) the L1.4 GHz–SFR relation (left) and the L1.4 GHz–LIR relation (right) are shown. Middle panel: for scenarios (0), (3), and (4) (left) and (0), (5), and (6) (right) the L1.4 GHz–SFR relation is shown. Bottom panel: for scenarios (0), (7), and (8) the L1.4 GHz–SFR relation (left) is shown and for scenarios (0), (5), and (6) the L150 MHz–SFR relation (right) is shown.

In the text

Fig. 6. L150 MHz–SFR relation obtained for scenario (0) using three different prescriptions for vadv: vA (solid cyan line), vsw (solid pink line), and vb from Eq. (15) (violet dashed line). The green dashed and green dot-dashed lines represent the primary and secondary electron contributions, respectively, and the dotted violet line represents the free-free contribution. The grey shaded band indicates the 1σ confidence region of the fit (solid grey line) to the data sample (purple points).

In the text

	Fig. A.1. Stokes-I continuum image at 1.4 GHz for each galaxy in our sample. The flux scale is in Jy beam−1. The magenta contours indicate the 3σ level and the green solid contour the region selected for the flux integration. The synthesised beam is shown with a blue circle at the bottom left corner of each image.
In the text

Fig. B.1. SEDs computed from scenario 0, for typical cases of low (Ṁ∗ = 0.1 M⊙/yr, left) and high (Ṁ∗ = 10 M⊙/yr, right) SFRs. The coloured lines (solid for the primary CRs and dashed for the secondary CRs) are the individual contributions of different non-thermal radiative mechanisms (the colour-coding is the same as in Fig. 2). The dashed brown line is the thermal free-free emission; the grey shaded region represents the Fermi-LAT energy range; and the vertical dotted black lines show the studied radio frequencies (150 MHz and 1.4 GHz).

In the text