{'fraction_non_alphanumeric': 0.03629002136464823, 'fraction_numerical': 0.026623638563937072, 'mean_word_length': 4.351723035100344, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exploring the physics behind the non-thermal emission from star-forming galaxies detected in gamma rays

Star-forming galaxies (SFGs) emit non-thermal radiation from radio to gamma-rays. We aim to investigate the main mechanisms of global CR transport and cooling in SFGs. The way they contribute in shaping the relations between non-thermal luminosities and SFR could shed light onto their nature. We develop a model to compute the CR populations of SFGs, taking into account their production, transport, and cooling. The model is parameterised only through global galaxy properties, and describes the non-thermal emission in both radio and gamma-rays. We focus on the role of diffusive and advective transport by galactic winds, either driven by turbulent or thermal instabilities. We compare model predictions to observations, for which we compile a homogeneous set of luminosities in these radio bands, and update those available in gamma-rays. Our model reproduces reasonably well the observed relations between the gamma-ray or 1.4 GHz radio luminosities and the SFR, assuming a single power-law scaling of the magnetic field with the latter with index beta=0.3, and winds blowing either at Alfvenic speeds or typical starburst wind velocities. Escape of CR is negligible for>30 Mo/yr. A constant ionisation fraction of the interstellar medium fails to reproduce the 150 MHz radio luminosity throughout the whole SFR range. Our results reinforce the idea that galaxies with high SFR are CR calorimeters, and that the main mechanism driving proton escape is diffusion, whereas electron escape also proceeds via wind advection. They also suggest that these winds should be CR or thermally-driven at low and intermediate SFR, respectively. Our results globally support that magnetohydrodynamic turbulence is responsible for the dependence of the magnetic field strength on the SFR and that the ionisation fraction is strongly disfavoured to be constant throughout the whole SFR range.


Introduction
Star forming galaxies (SFGs) comprise one of the extragalactic γ-ray sources classes detected at GeV energies by the Fermi-LAT telescope (Abdo et al. 2010a;Lenain et al. 2010;Ackermann et al. 2012;Abdo et al. 2010b). Their emission is believed to be associated with young stellar populations, which explains its correlation with the 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). Beside SNRs, other sources related to star forming regions such as young massive star clusters or microquasars could also play an important role as CR acceleration sites, possibly up to 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 SFR. 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 the aforementioned interpretation is the observed correlation between the integrated γ-ray luminosity (L γ ) of SFGs and their SFR (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;Lacki et al. 2011;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 SFR, so that the CR transport in these environments is inferred to approach the calorimetric limit, namely all Article number, page 1 of 21 arXiv:2107.00823v2 [astro-ph.HE] 23 Sep 2021 A&A proofs: manuscript no. aa 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 −1 , the SFR range where proton calorimetry might be achieved is still under debate due to our current poor knowledge of the escape mechanisms.
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, P. 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.4GHz ) 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 IR emission. Bell (2003) argued that these two effects conspire to achieve a linear relation between IR luminosity and L 1.4GHz . At high SFR, 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, also thermal processes such as free-free emission and absorption characterise the radio emission of SFGs. In fact, numerous young stars create several Hii regions, whose hot ionised gas radiates thermally, contributing to the total radio emission. A spectral identification of such freefree 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 arises above a few tens of GHz. Many works have discussed the fraction of the thermal contribution in such band with respect to the nonthermal 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, D. J. B. et al. 2021), where the free-free absorption might become severe. Typically, the free-free absorption significantly reduces the emission up to ∼ 10 − 10 2 MHz in galaxies with moderate SFR (Schober et al. 2017;Chyży et al. 2018), and up to a few GHz in ultraluminous IR galaxies (ULIRG; Condon et al. 1991;Clemens et al. 2010).
Despite most of the radio emission is expected from CRs populating galactic disks, some additional radio emission pro-duced in SFGs might have different origin. In particular, a nonnegligible 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 disk via diffusion or advection in winds (Heesen et al. 2009;Romero et al. 2018). It is also possible a smaller contribution of emission coming from pulsar wind nebulae (Ohm & Hinton 2013) and millisecond pulsars (Sudoh et al. 2021); in particular, the latter can become important at low frequencies and low SFR. 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 the SFR. This multiwavelength approach allowed to better constrain 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. Indeed, the latter are inferred to be responsible for most of the radio emission at high SFR, 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. For instance: How does the CR escape from galaxies evolve with the SFR and which 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 in exploring possible dependencies of magnetic and radiation fields in SFGs? What is the actual role/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 nonthermal 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.4HGz and L 150MHz -SFR correlations for all SFGs detected in γ-rays so far. In particular, we update the L γ -L 1.4HGz 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 taking into account also 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 SFR. 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. Besides, 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 mentioned different wavelength.
The paper is organised as follows: In Section 2 we present the data set used and the luminosity-SFR correlations. In Sections 3 and 4 we present the model and the scale relationships used, taking as a starting point the ones proposed in K20. In Section 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 Section 6 we present our conclusions.

The Luminosity-SFR relation
We aim to 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, that 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. We note that, 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 from K20 the 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.

The L γ -SFR relation
We revisit the L γ -SFR correlation studied by K20 and reconstruct 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 compute 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, M33 and NGC 2403, data in the 0.1-100 GeV band are not reported in 4FGLR2; we thus compute 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 include 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 γ = AṀ m * , yields a value of m = 1.42 ± 0.17 and log A = 39.21 ± 0.15, consistent with those reported in K20.

The L 1.4 GHz -SFR relation
To obtain a radio-continuum luminosity sample centred at 1.4 GHz (L 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σ) background fluctuation value of ∼ 0.5 mJy beam −1 , consistent with what was 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 −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 1 (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, 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 µm band (L IR ) is usually used as a proxy for the SFR. A power-law fit to L 1.4GHz as a function of L IR , L 1.4GHz = AL m IR , 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 the ones reported by Bell (2003) for their sample of 162 SFGs, and are also consistent with the bestfit values for a much bigger sample of 1809 galaxies reported by Yun et al. (2001). A power-law fit to SFR data, L 1.4GHz = AṀ 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.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 (∆Ṁ * 1.4 = 2.8). The latter may be A&A proofs: manuscript no. aa   Israel et al. (1992) Kuehr et al. (1981) Hales et al. (1991) 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.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 γ-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.4GHz /10 21 W Hz −1 ) 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 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

The L 150MHz -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 (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 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 150MHz = AṀ m * , yields a value of m = 0.98 ± 0.13 and log A = 21.7 ± 0.15, with a dispersion of 0.5 dex. On 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, D. J. B. et al. 2021), reporting values are between m = 1.07 ) and 1.37 ). 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 ∼ 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 −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 ≈ 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.

The model
We develop a non-thermal emission model to understand the observed L 1.4 GHz -SFR correlation (Sec. 2.2), while also contemplating the L γ -SFR (Sec. 2.1) and L 150 MHz -SFR (Sec. 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 valid over a broad range of SFRs.

CR injection and transport
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 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 = 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 15 eV), and ζ 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 R SN = (83 M ) −1Ṁ * , consistent with the Chabrier (2003) initial mass function (IMF), and that a fraction ξ CR = 10% of the supernova energy (E SN = 10 51 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; 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) upscattering 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: Here τ esc is the characteristic CRs escape time computed as τ −1 where each τ 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) τ loss (E), with the benefit that it yields a more precise normalisation for injection indices α > 2 when cooling dominates over escape.

Secondary products
Here we present the formulae used to calculate γ rays and secondary electrons 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 For E i > 100 GeV the emissivity formula is: with σ pp the cross section for inelastic p-p scattering, where E th = 1.22 GeV is the threshold energy for π 0 production, L = ln (E/TeV) and κ = 0.5 is the process inelasticity.F γ and F e are given in equations (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): Here n is the number density of the ISM, m p and m π are the masses of the proton and the pion, respectively, E π is the pion energy, E i,min = E i + m 2 π c 4 /4E i , K π ≈ 0.17 is the fraction of kinetic energy transferred from the parent proton to the single pion,f γ (E γ /E π ) = 1 andf e (E e /E π ) are defined in equations (36-39) from Kelner et al. (2006) (see also Appendix C).

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 between ionised and total number densities. Thus, we compute the absorbed free-free thermal luminosity at a frequency ν, L ff (ν), following the Olnon (1975) formalism for an homogeneous cylinder seen faceon, valid in the Rayleigh-Jeans regime (hν kT e , where h and k are the Planck and Boltzmann constants, respectively), 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), where T e the ionised gas temperature andḡ ff is the mean Gaunt factor (e.g. Leitherer et al. 1995), 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 forming activity and supernova rate, most of the ionised component 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 −4 (Krumholz et al. 2020). We take this into account exploring scenarios where ρ is assumed to such a low value. 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 e = 10 4 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: where L 0 is the unabsorbed luminosity produced by non-thermal radiative processes in our model.

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 Sec. 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.

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: 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 2 σ SB T 4 ), 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 4.2. Magnetic field 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 β we explore three possible scenarios, β = 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 ∝ Σ 0.5 g , 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 (Σ SFR ∝ Σ 1.41 g ) 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 (πGΣ 2 g ) and the magnetic field pressure (B 2 /8π) 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 β = 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 0 = 100 µ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 ∼ µG in low SFR galaxies (Ferrière 2010;Jurusik et al. 2014) to ∼mG for ULIRGs (Robishaw et al. 2008).

Advection
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 M yr −1 with velocities up to 10 3 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 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,Ṁ * 2 M yr −1 , are characterised by a starburst-driven wind, v sw = 400 km s −1 (see also K20, Yoast-Hull et al. 2013;Peretti et al. 2019). For objects withṀ * < 2 M 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 A = B/ 4π m p n i . Putting all together, we define a global advection velocity over the entire SFR range as 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.

Diffusion
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 quasilinear 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 δ characterises the nature of the turbulence. In particular, we have δ = 1/3 for the Kolmogorov turbulence, and δ = 1, 0.5 for Bohm and Kraichnan. In this work we assume the Kolmogorov-like 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 c = (2HπR 2 R −1 SN τ −1 SNR ) 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 4 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 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 .
We give the explicit dependency of τ D on the SFR in Appendix B.
Article number, page 7 of 21  Table 2. Values of the free parameters of our model (index of B(Ṁ * ) relation, its normalisation, ionised gas fraction, temperature of the ionised gas) for all the scenarios considered.

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 −2 -10 3 M yr −1 ), according to the scaling relations discussed in Sec. 4. We compare the timescales of the different processes involved, and solve the transport equation (Eq. 3) following the prescriptions of Secs. 3 and 4, to determine the main characteristics of CR transport in different ranges of SFR. We present the results for this scenario in Sec 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. Such a 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 max γ , L max 1.4 GHz and L max 150MHz 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 of each luminosity-SFR correlation in Secs. 5.2, 5.3 and 5.5, respectively.

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  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 −1 , 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 2 M yr −1 . Electron calorimetry is achieved forṀ * 5 M yr −1 , at which Bremsstrahlung dominates the cooling of GeV electrons; higher energy electrons are coolingdominated 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 SFR,Ṁ * 1 M yr −1 , most of the protons escape via diffusion for E p 10 GeV and via advection at lower energies. At higher SFR, up to a few tens of M yr −1 , 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 as good proton calorimeters for SFR 30 M yr −1 .

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 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 adv = v A or v adv = v sw . Since v sw > v A (v A 50 km s −1 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.
Variations of β, B 0 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), (6)) modifies the CR escape through v 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 γ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 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 γ -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 −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 L 1.4GHz -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.4GHz 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Ṁ * > 100 M yr −1 (see Fig B.1).
The emission of primary electrons dominates the L 1.4GHz -SFR correlation in the low SFR range ( 1 M 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 30 M 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 10 2 M yr −1 , L 1.4GHz is dominated by the synchrotron emission, hereafter L 1.4 syn . In order to understand the shape of the correlation, we note that L 1.4 syn 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 Sec. 3.2). This leads to a complex dependence of L 1.4 syn on the SFR, which we parametrise as where η is the slope of the synchrotron luminosity in the logarithmic L 1.4GHz -SFR correlation. Notice that such component dominates the L 1.4GHz -SFR correlation up to ∼ 10 2 M yr −1 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 SFR, also the transport of protons becomes calorimetric and the slope of the correlation softens (see Appendix B for a detailed discussion). Notice 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 SFR 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.4GHz 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 M82. At high SFRs ( 100 M yr −1 ) 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 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.4GHz -SFR correlation. On the right panel of Fig. 4 we present observational data in the L γ -L 1.4GHz plane and compare them with the outcomes of our model; we obtain a good agreement throughout the whole SFR range.

Parameter dependence of the L 1.4GHz -SFR relation
In this section, we discuss the impact of parameter variations on the L 1.4GHz -SFR correlation.
Magnetic field index. As the thermal emission is unaltered by changes in β, hereafter we focus only on the synchrotron component. Top panels in Fig. 5 illustrate how the variation of β (scenarios (1) and (2)) affects the slope of the L 1.4GHz -SFR and L 1.4GHz -L IR correlations.
In particular, the larger the value of β, the steeper the slope η of the L 1.4GHz -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 SFR 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), looks somehow 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 0 (scenarios (3), (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.4GHz on the magnetic field strength is not trivial, and a change in the normalisation constant B 0 introduces an overall shift which slightly depends on the SFR. Despite this com-plexity, we show that the parameter B 0 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 0 = 100 is the one that best fits the data. The resulting correlations for B 0 = 50 and 500, limits 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.4GHz -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 −4 (violet lines) and ρ = 10 −1 (orange lines).
Although a value of ρ=0.1 produce 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 M82 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 multi-wavelength 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 −4 underestimates the data at low SFR 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 4 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 4 K might be consistent with the average temperature in these environments (Westmoquette et al. 2009). The variation of T e (scenarios (7), (8)) affects the thermal emission and absorption. Such behaviour is shown in Fig. 5 (bottom left panel) from which it can be concluded that values between T e = 10 4 and 10 5 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, τ ff ∝ T −3/2 eṀ 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.4GHz across the whole SFR range.

The L 150MHz -SFR relation
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 τ 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 150MHz -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 SFR 10 2 M yr −1 . Although the effects of the variation of parameters in the L 150MHz -SFR are similar to the L 1.4GHz -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 β = 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 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 ρ and/or T 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.

Conclusions
We investigated how the transport of CRs, produced by sources related to the star formation, shapes the observed L γ -SFR, L 1.4GHz -SFR and L 150MHz -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, M82, Arp 220, Arp 299, NGC 3424. We also compute 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 sub-sample 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.4GHz correlation and found a steeper slope m = 1.26 ± 0.10 (with L γ ∝ L m 1.4GHz ) with respect to the one reported by Ackermann et al. (2012).
-We obtained a slope of m = 0.98 ± 0.13 for the L 150MHz -SFR correlation (with L 150MHz ∝Ṁ 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 multi-wavelength 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.4GHz -SFR and L 1.4GHz -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. Indeed, 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 −4 since the hot gas would be advected in the wind. -A scaling law B = B 0 (Ṁ * /10 M yr −1 ) β 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 −1 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.4GHz -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 4 -10 5 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 −1 . Such a calorimetry can be assessed via observations in γ-rays. On the contrary, relativistic electrons, despite reaching also the calorimetric limit for SFR 5 M yr −1 , 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 might be 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 have used the K20 value (0.027 M yr −1 ) 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 of the SFR of the SMC of 0.5-1 dex results almost always in a better agreement between observations and our L 1.4GHz -SFR and L 150MHz -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 AGNi 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.4GHz correlation, which is important for two reasons. The first one is that it could provide an independent method for the classification of γ-ray galaxies, in addition to the standard one based on the spectral index determination. The second one 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 150MHz -SFR correlation is important to assess the extent to which L 150MHz traces star formation. Furthermore, a more refined L 150MHz -SFR correlation would allow 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/acceleration region (such as a co-existing AGN) not directly related to the SFR .
The next generation of telescopes with better angular resolution and sensitivity will allow unravelling 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 SFRluminosity correlations in SFGs is a key to assess 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 et al. 2019). The latter can provide a new powerful tool for investigating the CR transport properties.  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.