Constraining the hadronic properties of star-forming galaxies above 1 GeV with 15-years Fermi-LAT data

,

These results are typically interpreted as an evidence for the existence of common properties shared by the entire population of SFGs and SBGs [6].For instance, the fact that these sources present hard power-law spectra E −γ , with γ ∼ 2.2 − 2.3, might indicate that the physics of CRs is dominated by energy-independent mechanisms, such as the pp inelastic timescale or advection [16].However, some works [17][18][19] have recently proposed that CR transport inside these source might be dominated by diffusion for E CR ≥ 100 GeV, leading to a suppression of the γ-ray and neutrino production rates at higher energies.Physically, this means that the calorimetric fraction F cal -i.e. the fraction of high-energy CRs which actually lose energy inside SFGs and SBGs producing γ-rays and neutrinos -might be smaller than previously predicted and also energy-dependent.In order to discriminate these two scenarios, however, new and more precise measurements (especially in the TeV energy range) are required [20].Several authors have attempted in modelling the calorimetric fraction along a large star formation rate (SFR) range (10 −2 − 103 M ⊙ yr −1 ) [5, 6, 8-12, 19, 21, 22].All the theoretical studies point towards the conclusion that F cal increases with SFR and so CRs lose most of their energy in the environment of SBGs.However, the actual degree of calorimetry is still under debate due to a lack of knowledge of the CR escape mechanisms from these astrophysical environments [6].In this paper, we provide new data-driven constraints on the calorimetric fraction of SFGs and SBGs by analysing a catalogue of 70 sources introduced by [12], using ∼ 15.3 years of Fermi-LAT data 1 and the publicly-available fermitools. 2 In particular, we search for γ-ray emission between 1 GeV and 1 TeV, dividing the catalogue in two samples: the 56 sources not yet detected 3 and the 14 sources which have been previously detected.We find strong hints of γ-ray emission in coincidence of M 83 and NGC 1365 at level of ∼ 4σ.Furthermore, for M 33 which was previously reported as a discovered source by [14], we find that its gamma-ray emission still stands right below the 4σ detection level (as already emphasised by [10,11]).
Then, we test a physically-motivated relation between F cal and the rate of supernovae explosion R SN , in contrast with the simplistic power-law function previously tested [12][13][14][15], finding a good agreement with the data.We emphasise that the correct estimate of the systematic uncertainty on R SN is crucial in order to extract the correct information on this correlation.Moreover, undiscovered sources place strong constraints to F cal , thus slightly modifying the F cal -R SN correlation.Therefore, in order to interpret these emissions as shared properties of all SFGs and SBGs is also important to take into account sources which present no evidence for γ-ray emission.
Finally, we employ such a correlation to evaluate the diffuse γ-rays and neutrinos flux from the whole source population.In order to do this, we make use of the recently-updated cosmic star-formation-rate distribution obtained through the James Webb Space Telescope reported by [23].We find that SFGs and SBGs might contribute (12±3)% to the extragalactic gamma-ray background (EGB) [24] above 50 GeV, while their contribution to the diffuse neutrino flux measured by IceCube with 6-year cascade events [25] might vary from 4% to 18% crucially depending on the assumed distribution of the spectral indexes along the source class.
The paper is structured as follows.In Sec. 2, we describe the sample of galaxies analysed.In Secs. 3 and 4, we describe the statistical analysis of the Fermi-LAT data and report the corresponding results, respectively.In Sec. 5, we discuss the theoretical model we adopt to evaluate the γ-ray and neutrino fluxes from each source.In Sec. 6, we describe the F cal -R SN correlation and discuss our findings.In Sec. 7, we extrapolate our results to the diffuse γ-ray and neutrino fluxes.Finally, in Sec. 8, we draw our conclusions.The paper has five appendices: in appendix A, we report all the new spectral energy distributions (SEDs) for the sources above the 5σ discovery threshold; in appendix B, we discuss the properties of the diffuse spectrum; in appendix C we comment on the impact of the systematic uncertainty affecting R SN on our results; in appendix D, we discuss the fit of data with a power-law function and the comparison with the function used in the main text and finally in appendix E, we discuss on the impact of the sources with potential AGN contamination in the data fit presented in the main text.

Sample of galaxies
We investigate the gamma-ray emission of 70 sources which we divide into two samples: • Sample A (see Tab. 1): it contains the galaxies introduced by [12] (see also [13]) for which no γ-ray detection has been reported yet.These galaxies exhibit a galactic latitude coordinate |b| ≥ 10 • and, therefore, the contamination from the diffuse galactic γ-ray emission is negligible.For these sources, we take the distances and the total infrared luminosity between 8 − 1000 µm (L IR ) from [13], consistently rescaled for the different hubble parameter H 0 used. 4 Sample B (see Tab. 2): it refers to the 14 sources discovered in γ-rays reported by [14], including also the Circinus Galaxy reported by [5,6].For this sample, we use the distances and the infrared luminosity reported by [5].For NGC 3424, ARP 220 and ARP 299, we use the updated values reported by [11].
Some of these sources are not only classified as SFGs but also AGNs, with Seyfert activity.
For this reason, we focus on γ-ray emission above E ≥ 1 GeV, where the photons from seyfert activity are expected to be negligible [26][27][28].

Data analysis
We analyse the latest Fermi-LAT data which have been collected in sky-survey mode from August 2008 and November 2023, from a Mission Elapsed Time 239557417 s to 720724699 s, with a total lifetime of ∼15.3 years.We select photons in the energy range [1 − 1000] GeV, which strongly reduces the possibility of mis-identification of sources due to a limited PSF dimension of ∼10 • at lower energies.We consider events belonging to P8R3_v3 version of the Pass 8 photon dataset and the corresponding P8R3_SOURCE_V3 instrument response functions.Data are analysed using the publicly-available fermitools provided by the Fermi-LAT collaboration and their analysis threads. 5 We consider a 15 • × 15 • Region of Interest (RoI) centred at the equatorial coordinates of each source, selecting only the data passing the filter for being considered of good-quality (DATA_QUAL>0)&&(LAT_CONFIG==1).
In order to reduce the contamination from the Earth's limb, following the default suggestions in the fermitools, the events with zenith angle >90 • are excluded.We emphasise that the Fermi-LAT collaboration has recently updated the selection for events above 1 GeV, selecting events for zenith angle <105 • [29].However, we have verified that the results do not change either for sample A or for sample B, even with this new selection.Therefore, we prefer to leave the event selection suggested in the fermitools in order to work with a photon sample with higher purity.
These data are analysed following the binned maximum likelihood ratio method, which is officially released by the Fermi-LAT collaboration.The likelihood function is defined as [30] where P(E i , X i | M i (Ω)) is the Poisson probability distribution function for observing a photon of a given energy E i and direction X i , given the expected number of photons M i provided by the model which depends on the Ω parameters.The index i runs over the bins for the events in the RoI.We determine the test statistic TS for each source as where L 0 is the maximised likelihood in the background-only hypothesis, namely in the hypothesis the source does not emit photons, and L is the maximised likelihood including the source under study.The conversion from the TS to the significance level can be performed using a chi-squared χ 2 distribution with degrees of freedom equal to the number of the free parameters for the source model [30].For instance, for power-law spectra, considering both normalisation and spectral index as free parameters, TS = 25 (also defined as discovery threshold for the TS) corresponds to ∼ 4.6σ significance.
In order to maximise the likelihood in Eq. (3.2), the data count maps are binned in angular coordinates, with 0.1 • bin per pixel, and in energy with 37 logarithmically spaced bins. 6he background hypothesis comprises all the sources in the 4FGL catalogue gll_psc_v32.xml[31,32], the standard isotropic extragalactic emission iso_P8R3_SOURCE_V3_v1 and the galactic diffuse emission gll_iem_v07.In order to account for the finite dimension of the PSF, we also consider sources outside the RoI with a further radius of 5 • .As suggested by [30], the fit is performed in an iterative way and at each step sources with very low TS, such as spurious solutions with TS < 0, are eliminated from the likelihood.
In this work, for the signal hypothesis, we test power-law spectra ϕ γ = ϕ 0 E −γ added at the nominal position of the source7 .In the likelihood maximisation, we fit all the sources leaving free the source parameters (normalisation ϕ 0 and spectral index γ) within 5 • of the RoI centre.Furthermore, we leave free the normalisation of extremely variable sources up to 15 • of the RoI centre8 as well as the normalisation for the isotropic extragalactic and the galactic diffuse templates.The other parameters are fixed to their best-fit values of the 4FGL catalogue.Finally, we also account for the energy dispersion using edisp_bins = -2 as advised in the fermitools threads.All the sources except for the Small Magellanic Cloud (SMC) and the Large Magellanic Cloud (LMC) are considered as point-like sources.For SMC and LMC, we instead utilise the official templates provided in the 4FGL catalogue.For these sources, we leave the source parameters to be free within 8 • and 6 • from the RoI centres, respectively.

Results of the statistical analysis
We report the obtained results in Tabs. 1 and 2 for the sample A and B, respectively.For each source of the sample A for which the TS is smaller than the discovery threshold (TS < 25), we report the luminosity distance, the infrared luminosity, and the 95% CL upper limit on the flux in the range 1 − 1000 GeV assuming a spectral index γ = 2.3 as typical value for known SBGs (see results for the sample B).We do not find any excess, except for M83 and NGC 1365 which shows TS ∼ 15.For these cases, we also report the best-fit values and the 68.3% CL limits in brackets.Differently from [33], we do not find any hint for NGC 3079: this is probably due to the fact that they look for photons with E ≥ 50 MeV where the limited Fermi-LAT PSF might cause mis-identification of sources.This problem has already been studied by [13] who pointed out that increasing the energy threshold leads to a better probe of the emission from single sources (and potentially reducing previous evidence of emission).Moreover, we find no evidence for γ-ray emission from the sources NGC 6946 and IC 342 which correlates with the most energetic CRs observed [34].
For each source of the sample B, we also report the best-fit interval of the flux normalisation and spectral index at 68.3% CL, and the corresponding value for the test statistics TS.Our results are in fair agreement with previous ones [14,31,35].For SMC, we find a slightly softer spectrum than [14] being in agreement with [31].For M31, along with the other sources of sample B, we have used the official point-like model present in the 4FGL catalogue, despite some other works have reported it as an extended source of 0.4 • [14,35].We obtain convergence anyway (with a TS ∼ 75), although with a very soft power-law spectrum ∼ E −3.0±0. 3.Finally, for M 33, there is not any match with sources present in the 4FGL catalogue.So, as for the sources in the sample A, we have added a point-like source in its position.Differently from [14], we find only an excess with (TS ∼ 16) which is below the discovery threshold.In appendix A, we report the SEDs for each source above the discovery threshold according to our analysis.We stress that some of the sources of sample B (LMC, SMC, M31, M82, NGC253) are not reported as simple power-laws in the 4FGL catalogue.Therefore, we show that there is not a statistical difference in using those signal models as opposed to simple power-laws.To this purpose, we define where L 4FGL is the likelihood maximised using the signal model in the 4FGL catalogue, while L PL is the likelihood maximised in the power-law model.We find that all the TS SM are much below the discovery threshold and so our signal assumption is justified.This result is given by the fact that the spectrum curvature is helpful to better describe the pion bump which is below 1 GeV.Furthermore, the Fermi-LAT sensitivity degrades above 10 GeV, leading to signal models being degenerate.This was also emphasised by [20] which highlighted the importance of the upcoming CTA to discriminate between different spectral assumptions for local SFGs and SBGs.For all the sources, we compute the γ-ray luminosity L γ between 1 − 1000 GeV, using where is the integration of the differential flux measured weighted by the energy, and z is the redshift of the source, directly related to the luminosity distance D L (z).Fig. 1 shows the L γ in the energy range [1,1000] GeV versus the L IR for the samples A and B. We report the best-fit values and the corresponding 1σ uncertainty for all the discovered sources as well as for the three sources which give us a 4σ hint of emission.On the other hand, for the undiscovered sources, we report the 95% CL upper limit assuming a E −2.3 spectrum.In the plot, we also take into account a 10% uncertainty in each distance and 5% in L IR as reported by [14].Table 1.Results for the sample A. From left to right: the source name, the luminosity distance, the infrared luminosity, the upper limit on the integrated γ-ray flux, the spectral index assumed to evaluate the upper limit.For hints of γ-ray emissions, we report the best-fit values in brackets.
Source Table 2. Results for the sample B. From left to right: the source name, the luminosity distance, the infrared luminosity, the integrated γ-ray flux, the flux normalisation, the spectral index, the value of the test statistics.† Note that since M 33 is below the discovery threshold, we also compute the 95% CL upper limit fixing γ = 2.3, obtaining F 1−1000 GeV = 1.65 Cyan points represent discovered sources with a significance level of σ ≥ 5, while orange points denote hint sources with a significance level of approximately σ ≃ 4. In both cases, the best fit scenarios along with their 1σ uncertainty are considered.For undiscovered sources (σ ≤ 4), represented by black triangles, we provide 95% confidence level upper limits, fixing γ = 2.3.

On the Non-thermal emission from SFGs and SBGs
The results presented in the previous section have important repercussion on the CR transport mechanisms occurring inside these sources.Indeed, since photons produced by hadronic interactions usually carry 10% of the parent energy of CRs, the γ-ray spectra are expected to inherit the properties of the CR distribution inside the sources.In order to assess such implications, we use a model describing the non-thermal emission of the sources.In general, since we expect the emission of SBGs to be dominated by their nuclei, we can neglect the spatial dependence of the CR diffusion.Hence, we can study the CR transport under the leaky-box model equation where the CR transport is modelled by a balance among different competing processes: the injection term of the sources such as SNRs, the escape phenomena (advection and diffusion) and the energy-loss mechanism such as hadronic collisions: where dE/dt = −E/τ loss , with τ loss being the energy-loss timescale, τ esc is the escape timescale, and Q(E) is the injection spectrum of SNRs.We assume the injected spectrum to be a powerlaw with a E max = 10 PeV exponential cut-off consistent with our previous results [4,20] and we neglect any other spectral feature of the injected spectrum (see below for further remarks about the chosen cut-off).The normalisation is set as Hence, the total energy injected into CRs is η = 10% of the total E SN = 10 51 erg emitted by SNRs.The quantity R SN is the SNRs rate which is expected to be tightly connected to the infrared luminosity according to the empirical relation [6,36] which takes advantage of the Chabrier Initial mass function (IMF), consistent with 83 M ⊙ , converted in new stars for each supernova explosion.In other words, the SFR is connected to ).We emphasise that Eq. 5.3 is not linear because the infrared luminosity itself is not a perfect tracer of the SFR.
In general, the solution to Eq. 5.1 can be approximated as [6] N where τ tot = (τ −1 loss + τ −1 esc ) −1 and the last passage holds for Q(E) ∝ E −γ .For SBGs, pp interactions should be the dominant CR energy-loss mechanism.Therefore, τ loss = τ pp .In turn, the escape timescale is given by the competition between CR advection and diffusion phenomena.While it is expected that their relative contribution to change across the whole SFR range (10 −2 − 10 3 M ⊙ yr −1 ) [8-10, 21], these timescales are strongly model-dependent as well as dependent on the assumption for their scaling relation with the SFR.Indeed, although Refs.[1][2][3][4]20] have shown advection to be important as escape phenomenon for SFGs and SBGs, Refs.[17][18][19] have argued that advection should be suppressed in interstellar medium (ISM) ambient in SBGs, proposing a major role played by diffusion phenomena.Furthermore, whereas Refs.[1][2][3][4]20] have modelled the diffusion coefficient using quasi linear theory assuming a pre-existent magnetic field turbulence, Refs.[17][18][19] have used self-generated diffusion from streaming instability.Given that it is not possible to distinguish between these scenarios with Fermi-LAT data (see previous section), here we introduce an overall parameter-F caldefined as in order to test if the γ-ray measurements of sample A and sample B might be interpreted in terms of star-forming activity.F cal is defined between 0 and 1 and it can be interpreted as an average fraction of CRs between 10 ≤ E CR /GeV ≤ 10 4 actually losing their energy onto pp collisions producing γ-rays and neutrinos.A very small F cal value would correspond to a very strict constraint on the ability to confine high-energy protons by the source.F cal can be expressed as For the following analysis, we assume F cal to be constant, which allows us to estimate it directly from the γ-ray data without any assumption on the magnetic field, gas density, wind velocity and energy dependence of the diffusion coefficients of the sources because we treat it as an effective number for each of the source in the sample.However, this restricts our analysis to assume that τ esc is only mildly energy-dependent in the whole SFR range analysed, leading to negligible diffusion phenomena.This might slightly overestimate F cal for low-SFR sources, where the role of diffusion might be more relevant [6].However, we stress that all the SFGs discovered, from SMC to ARP 220, show the same spectral behaviour (∼ E −2.3−2.4 ) totally consistent with the injected spectrum inferred for the Milky-way [37,38].Therefore, an energy-independent escape timescale cannot be, at the moment, completely ruled out.From Eq. 5.5, we can quantify the photon production rate following the analytical procedure of [39] (see also [6]).For E γ > 100 GeV, we have where Fγ (x, E γ /x) is defined in [39] (see Eqs. (58-61)) and x min = 10 −3 .For lower energies, we can assume that the pions produced by pp collisions take K π = 17% of the kinetic energy of the parent high-energy proton (delta-function approximation), having At E γ = 100 GeV, Eq. (5.8) is scaled in order to match Eq. (5.7).The final γ-ray flux at Earth is given by where z is the redshift of the source, D L (z) is the luminosity distance, and τ (E, z) is the optical depth for photons travelling through EBL and CMB.For the computation of the opacity, we employ the model of [40].From pp interactions, we expect production of highenergy neutrinos as well and we estimate their flux using the same procedures as for γ-rays.
In particular, for E ν > 100 GeV, we have that where Fν 1 µ (x, E ν /x), Fν 2 µ (x, E ν /x) and Fνe (x, E ν /x) take into account all the neutrinos produced in the interactions and are defined by [39].The factor 1/3 is due to the fact that we expect an equal flavour ratio at Earth.The final neutrino flux is given by Before concluding this section, we emphasise that the γ-ray spectra of SFGs and SBGs might be contaminated also by leptonic contributions such as Inverse compton and Bremsstrahlung as well as by the AGN related activity hosted by some of the sources in sample A and B [41][42][43][44] (see app.E).Therefore, our results might be relatively interpreted as upper limits for F cal which corresponds to conservative constraints on the star-forming activity of the sources.
Regarding the leptonic contributions, in our approximation where diffusion is negligible, the contribution from leptonic photons is expected very limited [4] above 1 GeV.However, this intrinsically assume a proton to primary electron ratio of K ep = 50 consistently with the Milky-way.Indeed, lower values would lead to a major role for primary electrons since they are usually trapped in the SFG environments cooling down much faster than protons [1].
6 On the correlation between gamma-rays and star-forming activity In this section, we discuss our constraints on the calorimetric fraction F cal from γ-ray observations and its correlation with R SN .Previous studies [12][13][14][15] have tested the relation . However, this relation cannot be valid for a wide SFR range [10 −3 − 10 3 ] M ⊙ yr −1 , since the calorimetric limit cannot be exceeded.In order to test a physically motivated relation between F cal and R SN , we exploit the fact that τ pp = (k • n gas σ pp c) −1  where k = 0.5 is the mean inelasticity of the process and τ esc = H/v wind with H being the height of the nucleus and v wind is the velocity of the galactic winds.Both n gas and v wind are expected to scale with R SN .Indeed, according to the kennicutt relation [45,46], we have a strict connection between n gas and R SN , namely SN [3,20].By contrast, the wind velocities have been found to correlate with the SFR as v wind ∝ R 0.15−0.30SN [47].All of this leads to τ esc /τ pp ∼ AR 0.30−0.50SN .Therefore, in the present paper we probe the following relation between F cal and R SN with A and β free parameters to be deduced from data.We notice that, for small value of R SN , Eq. (6.1) becomes consistent with a pure powerlaw relation as tested by previous study (in appendix D, we discuss also the power-law fit).
In order to test Eq.(6.1), for each source we estimate R SN from the infrared luminosity according to Eq. ( 5.3) and we calculate the calorimetric fraction as described in Eq. (5.5) by matching the measured integrated spectrum with the theoretical one using the model described in the previous section.
For the discovered sources (sample B), we evaluate the best-fit scenario and the ±1 σ values.For the undiscovered sources (sample A), we utilise the best-fit scenario for the fixed γ = 2.3 and for the uncertainty, we consider the difference between F 95%CL and F best .
In addition to the statistical uncertainties inferred by Fermi-LAT data, we also take into account the systematic uncertainties affecting F cal .On this regard, uncertainties on the source distance and rate of supernovae explosions as well as the detector systematics play a crucial role.As we mentioned above, the distance and the infrared luminosity provide an uncertainty of the order of 10% and 5%, respectively.By contrast, the uncertainty on R SN might also come from the IMF and the amount of mass converted in new star from each supernova explosions.The total uncertainty on R SN is difficult to reliably assess and it may vary within 20 − 40% (see [12] for further details).For the following discussion, we consider a systematic uncertainty of 20% on our estimates of R SN and in the appendix C we discuss the impact of a higher uncertainty.Regarding the detector systematic uncertainty, we consider a conservative uncertainty of 10%.9 Summing all the systematic uncertainties in quadrature, we obtain an overall 30% uncertainty on each value of F cal .
Fig. 2 shows the obtained F cal both for undetected and detected sources as a function of R SN , as well as the 1σ bands from the fit of Eq. (6.1) according to two different samples of galaxies: • Discovered sources, for which we find A = 2.2 ± 0.8 and β = 0.55 ± 0.08; • Combined sources, namely discovered + undiscovered sources, for which we find A = 0.7 +0.3 −0.2 and β = 0.39 ± 0.07.
Interestingly, even though the undiscovered sources are characterised by higher uncertainties, they are anyway able to constrain the F cal fit especially in the range [0.01 − 1] yr −1 .In the lowest range for R SN , the fit is totally dominated by the galaxies of local group (SMC, LMC, M 31 and M 33).On this regard, we have verified, in app.E, that removing M31 from the fit (because of its soft spectrum) does not impact our results.We a posteriori notice that our results are completely in agreement with the expected scaling values for CR escape phenomena dominated by advection.We have also verified that our results for low SFRs are consistent with the γ-ray measurements of the central part of our galaxy, the central molecular zone (CMZ).In fact, the Fermi-LAT and the HESS collaborations have reported ∼ E −2.4 spectra for the CMZ [48,49].In particular, assuming the observed R SN = 2 • 10 −4 yr −1 [50] for the CMZ and the corresponding F cal from our combined fit, we obtain that our predicted γ-ray flux is consistent with the diffuse measurements from the galactic ridge of the CMZ.We mention that a greater component coming from leptonic processes would lead to a smaller calorimetric fraction potentially constraining even more the properties of the population.We can extract that F cal ≳ 50% for SFR ≳ 190 +1230 −130 M ⊙ yr −1 when considering the whole sample.This might reduce the degree of calorimetry of Ultra Luminous Infrared Galaxies (ULIRGs) (sources with SFR ≳ 100 M ⊙ yr −1 ), although, this information at the moment is mainly driven by galaxies with lower IR luminosity, since Fermi-LAT is not yet sensitive enough to directly probe the Cyan points denote discovered sources, whereas black points denote undiscovered sources.Specifically, for sources exhibiting a flux compatible with zero, we present 95% CL upper limits indicated by black triangles.For all the sources, F cal and R SN values are reported respectively with 30% and 20% systematic uncertainties.We also report the best-fit and the corresponding 1σ uncertainty band of the fit performed over the whole sample (orange) and over discovered sources (blue).
calorimetric scenario within ULIRGs, due to their large distances.We highlight that MHD simulations (e.g.[8]) have theoretically predicted that calorimetric limit (F cal ≃ 1) cannot be reached by SBGs, although this conclusion is driven by an assumed diffusion coefficient of 3 • 10 28 cm 2 s −1 at 3 GeV, which is higher than expected in extreme environments such as ULIRGs [1].On the contrary, our results are entirely driven by the latest data, making them the most current and robust constraints.

Extrapolation to the diffuse emissions
We can use the calorimetric fraction of local SFGs and SBGs evaluated in the previous section to constrain the diffuse non-thermal emission of the entire source population.The diffuse emission, per solid angle, is given by where z is the redshift, E(z) = Ω M (1 + z) 3 + Ω Λ , S SFR (L IR , z) is the density of the sources as a function of the infrared luminosity, Q γ,ν are the γ and neutrino production rate for each source, and τ ν = 0 and τ γ (E, z, L IR ) accounts for the CMB+EBL absorption of photons as well as for internal absorption phenomena [3].We highlight that in Eq. 7.1 we use L IR = 10 6 L ⊙ as a lower limit for the infrared luminosity corresponding at SFR ∼ 4 • 10 −3 M ⊙ yr −1 .Increasing such a value to 10 10 L ⊙ (∼ 1 M ⊙ yr −1 ) results in a reduction of the flux by ∼5% only, since the bulk of the emission comes from sources with higher star formation rates.For the density of the sources, we use the approach described by [23], who have recently updated the distribution of the cosmic SFR using also JWST data.The distribution is given in terms of a Schechter function which behaves as a power-law for L IR ≪ L * (z) and as a Gaussian in log 10 (L IR ) for L IR ≫ L * (z).The redshift parameter evolutions are not simply set by power-laws, but rather follow skew Gaussian distributions [51] where erf(x) is the error function, k is called the shape parameter, ω the scale factor, A L and A Φ are the normalisation for the evolution of L * (z) and Φ * (z), respectively.Eqs.(7.3) and (7.4) provide physical representations of the evolution, allowing for different peaking redshifts as well as asymmetric increasing/decreasing rates for several populations [51].In fact, one of the main advantages of such a parameterisation is that it can be divided for distinct source classes.Here, we consider SFGs and SBGs, taking the values reported in Tab. 3.They provide excellent agreement with the ones reported by [23] (see their Fig.10).Some parameters are also in agreement with the ones reported by [51].Finally, for F cal (L IR ) ).The first four parameters (from left to right) are the ones reported by [23], while the remaining parameters are obtained to match the SFG and SBG distributions reported in their Fig. 10.
we use Eq.(6.1) with parameters inferred by the data of both the discovered sources and the total sample (discovered + undiscovered) sources.Fig. 3 shows the final γ-ray (in dark red colour) and neutrino (in cyan colour) fluxes for the combined fit including discovered and undiscovered sources.
On the left panel, we have fixed spectral index to γ = 2.3, while on the right we have used a spectral index distribution (blending scenario) provided by a superposition of Gaussian distributions with mean values equal to the best-fit spectral index for discovered sources and with standard deviation equal to their corresponding uncertainty (for the spectral index blending flux calculation, we employ the same technique as in [4]).In this approach, the injected spectral index follows a continous distribution which allows also for spectral indexes lower than 2. The theoretical predictions are compared with the Isotropic Gamma-Ray Background (IGRB) measured by Fermi-LAT [24], the 6-year cascade neutrino flux [25] and 7.5-year HESE data [52] measured by the IceCube neutrino Observatory.The fluxes are dominated by distant sources with a contribution peaking at z ≃ 1.Furthermore, the bulk of On the left, the spectral index is fixed at 2.3 for each source.Right: the same but considering a spectral index blending.In both panels, the fluxes are compared with the Isotropic Gamma-Ray Background (IGRB) measured by Fermi-LAT [24], the 6 year Cascade neutrino flux [25] and 7.5 year HESE data [52] measured by the IceCube neutrino Observatory.
the emissions come from ULIRGs saturating almost 51% of the emissions (see the appendix B for details).We find that the total contribution to the extra-galactic gamma-ray background (EBG) [24] between 50 GeV and 2 TeV is ≃ (12 ± 3)%, almost independent on the spectral index distribution considered.The neutrino spectrum, on the contrary, is strongly dependent on the spectral index distribution.Indeed, fixing a spectral index γ = 2.3 provides a soft diffuse spectrum which can explain only (4 +1 −2 )% of the 6-year cascade flux between 10 TeV and 1 PeV.By contrast, the spectral index blending hardens the spectrum and allows for the neutrino spectrum to explain (18 +3 −5 )% of the 6-year cascade IceCube flux.This result is mainly driven by sources with γ ≲ 2 which contaminate the overall distribution of 10%.Indeed, if we only considered the distribution with γ ≥ 2, the neutrino spectrum would be at level of ∼ 7% of the 6-year cascade IceCube flux, reducing the observable signature of the spectral index blending.We notice that, at the moment, some observed γ-ray spectra of young SNRs might point to very hard injected proton spectra [53,54], although it is still controversial if this is a true signature given by hard hadronic spectra or leptonic processess.We also underline that, given the limited number of discovered sources, it is not possible to derive a robust distribution for the spectral indexes and its impact might vary also with respect to the the statistical treatment of the data [4].The neutrino flux is also sensitive on the chosen high-energy cut-off for CRs.Indeed, Ref. [2] has argued that the highly dense environment of SBGs might cause turbulent amplification of the magnetic field, leading to an E max ≃ 50 − 100 PeV.Furthermore, it is possible that since E max is correlated to the magnetic field value in SBGs, it might have a non-trial dependence on the SFR leading to a further signature, which we leave for future explorations.
Our results are completely consistent with previous works [12,14] which employ the same technique and also with our previous multi-messenger analysis [4].We find that SFGs and SBGs contribute significantly less to the EGB than the limits imposed on non-blazar sources [55], which has sensibly reduced the possible role of SFGs and SBGs to the EGB suggested by earlier works [56].On the other hand, we find a slightly lower contribution than [2] due to several reasons.Firstly, we assume a steeper injected spectrum; secondly, the authors of Ref. [2] have assumed the background photon energy density to be equal to the Left: Diffuse 2σ γ-ray (dark red) and neutrino (cyan) bands predicted with the fit over discovered sources.On the left, the spectral index is fixed at 2.3 for each source.Right: the same but considering a spectral index blending.In both panels, the fluxes are compared with the Isotropic Gamma-Ray Background (IGRB) measured by Fermi-LAT [24], the 6 year Cascade neutrino flux [25] and 7.5 year HESE data [52] measured by the IceCube neutrino Observatory.
M82 value in order to estimate the contribution of internal absorption, while we take into account the fact that the energy density of background photon linearly scales with R SN [3].This leads to a further suppression of photons for high SFR sources.Thirdly, we assume a lower value for the high-energy cut-off for CRs.Moreover, there is a different assumption on F cal , since the authors have estimated only the contribution of high SFR sources assuming that they all had the same F cal as M82 which is an assumption mainly tuned on discovered sources.On this regard, we assess the impact of the undiscovered sources in the F cal , in Fig. 4 where we show the diffuse γ-ray and neutrino spectra obtained with the fit of the discovered sources only.In this case, we obtain that SFGs and SBGs may contribute more to the EGB and the diffuse neutrino flux (∼ a factor 2) compatible with the estimates given by [2].Therefore, undiscovered sources are not only important to correctly estimate the significance of the correlation between F cal and R SN , but they are also necessary to correctly extrapolate information to the whole source population [12].This is crucial because, typically, analyses which attempts to constrain the properties of SFGs and SBGs tune their models on the sources discovered in the γ-ray range, but there are a lot of sources with the similar astrophysical properties which have not been detected and they should be taken into account if the entire source population share the same properties.
However, SFGs and SBGs are still unable to completely saturate the IGRB between ∼ 1−100 GeV as recently obtained by [18].Also in this case, the difference with our approach is given by several factors such as the source count and the CR transport model.Moreover, their assumed CR transport allows for F cal being a function of redshift (see Fig. 4 in the extended data section in [18] ).Indeed, if distant sources, which dominate the diffuse background (see app.B), are more intrinsically luminous, then SFGs and SBGs are allowed to explain a higher portion of the diffuse fluxes.We emphasise that in our approach Eq. (6.1) is considered to be valid at each redshift even though only local sources has been used to constrain it.only future observation can challenge this assumption because at the moment Fermi-LAT is not sensitive enough to probe the calorimetric scenario for more distant sources.We point out that even though F cal calculated with Fermi-LAT data corresponds to average values of E CR between 10 − 10 4 GeV, we extrapolate this calorimetric fraction also to higher energies in order to estimate the neutrino contribution.This, from one hand, it may be pessimistic since in case of energy-independent escape timescales, F cal is logarithmically energy-increasing due to the energy behaviour of σ pp .From the other hand, at PeV energies, the diffusion process might not be negligible, leading to escape timescales might be energy dependent strongly suppressing the calorimetric fraction.On the whole, we find our approximation to be a reasonable trade-off, although it is difficult to quantify the uncertainty on the neutrino flux given the uncertainty on the nature of the diffusion process.

Conclusions
In this paper, we have analysed 70 local sources, classified as star-forming and starburst galaxies, using 15 years of Fermi-LAT data.In order to reduce contamination from possible AGN activity as well as to reduce the possibility of mis-identification of sources from limited PSF, we have searched for photons with E ν > 1 GeV.We have found evidence at ∼4σ for two nearby sources, M 83 and NGC 1365.On the contrary, even with 15 years of Fermi-LAT data, M33 still stands at 4σ due to an improved treatment of the background model.
We imposed strict upper limit at 95% CL fixing a spectral index γ = 2.3 for the other sources.Exploiting these findings, we have then revisited the correlation between the γ-ray luminosity and the star formation rate for local star-forming and starburst galaxies.For the first time, we have studied this correlation under a physically-motivated relation between the calorimetric fraction and the rate of supernova explosions.We have found that there is a good agreement between the measurements and the theoretical model and that undiscovered sources play an important role in constraining the calorimetric fraction.This is crucial in order to capture the shared properties of these sources.
Then, we have extrapolated this information to constrain the diffuse γ-ray and neutrino spectra of SFGs and SBGs, finding that they contribute about 12% to the EGB above 50 GeV.The corresponding neutrino flux is strongly dependent on the spectral index distribution along the source class.Indeed, if it is fixed at γ = 2.3 for the entire source spectrum, the contribution is negligible to the diffuse neutrino flux measured by ICeCube.
By contrast, if there is a continuous distribution of this parameter within the source class, the contribution to the diffuse neutrino flux could increase by up to 20% because of sources with hard spectra.Therefore, future measurements, which aim to expand the sample of galaxies above the discovery threshold, will be essential to test how this parameter varies across the SFGs and SBGs population and to quantify its impact on the diffuse neutrino flux.
Finally, with current Fermi-LAT data we have obtained that high SFR sources have F cal ≳ 50%, which is theoretically expected but further data can challenge this concept leading to even a smaller calorimetric fraction.This is crucial because they mainly drive the diffuse γ-ray and neutrino fluxes.Hence, future analyses and data aiming at directly probing the degree of calorimetry of these sources are fundamental to further constrain the diffuse emission of SFGs and SBGs.

A Spectral energy distributions
Here we report the spectral energy distributions (SEDs) for the sources above the discovery threshold.We divide the analysed energy range [1 − 1000] GeV in 9 independent bins (3 per decade) and perform a likelihood analysis in each bin fixing the spectral index to γ = 2.0.If TS < 4, then we report the upper limit at 95% CL.Our results are shown in Figs. 5 and 6, where we divide the sources in the northern hemisphere (equatorial declination (δ > 0 • ) and in the southern hemisphere (δ < 0 • ).The red points correspond to the best-fit Fermi-LAT measurements with the 1σ uncertainty, while the black line and the grey band respectively represent the best-fit and the 1σ band for the fit over the entire energy range.For M 82, NGC 253 and NGC 1068, we also report the measurements (in blue color) taken by VERITAS [57], H.E.S.S. [58] and MAGIC [59], respectively.Finally, for each source (from Tab. 4 to Tab. 16), we report the obtained TS in each energy bin and if TS > 4, we report the best-fit value of the SED and its 1σ uncertainty, otherwise we report its 95% CL upper limit.).The red points corresponds to the SED points for each energy bin, while the black line and the grey bands respectively correspond to the best fit and 1σ band for the fit over the entire energy range.For NGC 253 and NGC 1068, we also report with blue data points the H.E.S.S. measurements [58] and MAGIC upper limits [59], respectively.Table 5. NGC 253: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.
Energy range TS E 2 dF/dE 1σ error 95% CL upper limit   Table 8.Circinus: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.
Energy range TS E 2 dF/dE 1σ error 95% CL upper limit            Table 16.LMC: the energy range corresponding to each bin, the TS obtained for the source, the the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.

B Redshift distribution and the role of ULIRGs
Here, we assess which are the redshift and star formation rate values corresponding to the largest contribution to the diffuse γ-ray and neutrino fluxes.We focus our attention to the diffuse neutrino flux, because neutrinos are not absorbed by the EBL, therefore they maintain the information of the redshift distribution.Furthermore, we fix γ = 2.3 and E ν = 1 TeV, since the results do not change either in terms of the spectral index or the energy.On this regard, we notice that even though the flux redshifting impacts the high-energy cut-off leading to different conclusions for energies near the cut-off, the final SED is maximum for E ≲ 10 TeV, making our approximation reasonable.
In the left panel of Fig. 7 we show the redshift distribution of the differential flux, once integrated over the luminosity.It represents the neutrino flux coming at different redshift from the sources of all luminosities.The maximum of the distribution stands for z ≃ 1, which represents the maximum of the cosmic star formation rate distribution [23].In the left panel of Fig. 7 we show the dependence of the differential flux over L IR , integrating over all the redshifts.Hence, it quantifies the contribution from all the sources having a given IR luminosity.We stress that even though the maximum of the differential flux stands for the lowest values of L IR , the ULIRGs are the ones which contribute most to the total flux.Indeed, the integration over sources with SFR > 100 M ⊙ yr −1 (L IR > 7.2 • 10 11 L ⊙ ) provides about 51% of the total spectrum.Sources with 1 M ⊙ yr −1 < SFR < 100 M ⊙ yr −1 contribute for 44% and the remaining 5% is due to sources with lower star formation rates.Therefore, correctly assessing the calorimetric budget of ULIRGs is fundamental in order to derive correctly the contribution of the entire source population.
C Impact of the systematic uncertainty on R SN Here we discuss on the impact of the systematic uncertainty on R SN .To this purpose, we assume that the systematic uncertainty on R SN is ∼ 45% instead of 20% as adopted in the main analysis.Once summed in quadrature with the 10% uncertainty on the distance, it leads to a systematic error of 50% on F cal .Using this systematic error, we perform again the fit of the relation in Eq. (6.1) and reports the results in Fig. 8.
For the discovered sources, we find A = 1.9 ± 1.1 and β = 0.54 ± 0.11, while for the combined sample A = 0.5 +0.3 −0.2 and β = 0.36 ± 0.11.The fits are totally consistent within the 1σ with the ones presented in the main text.Furthermore, even though larger uncertainties increase the statistical error of the fits, we still find that undiscovered sources are able to constrain F cal reducing its value especially in the range [0.2−1] yr −1 at 1σ level.We emphasise that future measurements will be able to reduce the uncertainty on the supernovae explosion rate and they will provide us with a much more constrained correlation function leading to a smaller uncertainty on the diffuse emissions of SFGs and SBGs.

D Power-Law Fit and Comparison with the physically-motivated expression
In this section, we discuss the fit using a simple power law expression.Considering only the discovered sources, we obtain A = 0.8 ± 0.14 and β = 0.39 ± 0.04, while for the whole sample, we obtain A = 0.45 +0.11 −0.08 and β = 0.32 ± 0.05.Fig. 9 shows the obtained F cal in terms of R SN for these two fits (best-fit and 1σ band).
Given the enormous uncertainties on the data, at the moment, there is not any statistical preference for the power-law fit or for Eq.6.1.In particular, we find the following the reduced chi squared: • Discovered sources χ 2 /d.o.f = 1.58 (Eq.6.1), χ 2 /d.o.f = 1.66 (power-law) • Combined sources χ 2 /d.o.f = 0.81 for both Eq.6.1 and the power-law case Hence, we conclude that from the data standpoint, it is not possible to exclude a fit from the other, but we advise against the power-law fit, since it is not physical and for high-SFR sources allows for F cal > 1.
E Impact of Sources Containing AGN contamination and Hint Sources Some of the sources considered in sample B are highly likely contaminated by the AGN activity hosted by the galaxies.These sources are: NGC 2403, NGC 3424 and Circinus Galaxy [14].Therefore, in this appendix, we discuss how much these sources impact the result of our fits.Furthermore, the three hint sources (M83, M33, NGC 1365) might impact our results in the sense that with more data and a better background description, their detection significance might reduce, in principle changing the results shown in the main text.Furthermore, the M31 galaxy is characterised by a pretty soft spectrum which might not completely be compatible with a CR spectrum dominated by advection.As a result, we remove these 7 sources from the fit and perform the same analysis as shown in the main text.For the discovered sources, we find A = 2.4 ± 1.1 and β = 0.58 ± 0.08 with a a reduced chi square at the best-fit value of χ 2 /d.o.f.= 1.3.For the combined sample, A = 0.7 +0.3 −0.2 and β = 0.38 ± 0.08 with χ 2 /d.o.f.= 0.67 at the best-fit value.These results are completely consistent with the ones presented in the main text, and consequently, our conclusions remain unchanged.

Figure 1 .
Figure1.The γ-ray luminosity as a function of the total infrared luminosity for the entire sample.Cyan points represent discovered sources with a significance level of σ ≥ 5, while orange points denote hint sources with a significance level of approximately σ ≃ 4. In both cases, the best fit scenarios along with their 1σ uncertainty are considered.For undiscovered sources (σ ≤ 4), represented by black triangles, we provide 95% confidence level upper limits, fixing γ = 2.3.

1 RFigure 2 .
Figure 2. F cal in terms of R SN for the whole sample.Cyan points denote discovered sources, whereas black points denote undiscovered sources.Specifically, for sources exhibiting a flux compatible with zero, we present 95% CL upper limits indicated by black triangles.For all the sources, F cal and R SN values are reported respectively with 30% and 20% systematic uncertainties.We also report the best-fit and the corresponding 1σ uncertainty band of the fit performed over the whole sample (orange) and over discovered sources (blue).

Figure 3 .
Figure 3. Left: Diffuse 2σ γ-ray (dark red) and neutrino (cyan) bands predicted with the fit over the whole sample.On the left, the spectral index is fixed at 2.3 for each source.Right: the same but considering a spectral index blending.In both panels, the fluxes are compared with the Isotropic Gamma-Ray Background (IGRB) measured by Fermi-LAT[24], the 6 year Cascade neutrino flux[25] and 7.5 year HESE data[52] measured by the IceCube neutrino Observatory.

Figure 4 .
Figure 4. Left: Diffuse 2σ γ-ray (dark red) and neutrino (cyan) bands predicted with the fit over discovered sources.On the left, the spectral index is fixed at 2.3 for each source.Right: the same but considering a spectral index blending.In both panels, the fluxes are compared with the Isotropic Gamma-Ray Background (IGRB) measured by Fermi-LAT[24], the 6 year Cascade neutrino flux[25] and 7.5 year HESE data[52] measured by the IceCube neutrino Observatory.

Figure 5 .
Figure 5. Collection of the SEDs for sources situated in the southern hemisphere (δ < 0 • ).The red points corresponds to the SED points for each energy bin, while the black line and the grey bands respectively correspond to the best fit and 1σ band for the fit over the entire energy range.For NGC 253 and NGC 1068, we also report with blue data points the H.E.S.S. measurements[58] and MAGIC upper limits[59], respectively.

Table 4 .
M 82: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.

Table 6 .
ARP 220: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.

Table 7 .
NGC 1068: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.

Table 9 .
SMC: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.

Table 10 .
M 31:the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.

Table 11 .
NGC 2146: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.

Table 12 .
ARP 299: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.

Table 13 .
NGC 4945: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.

Table 14 .
NGC 2403: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.

Table 15 .
NGC 3424: the energy range corresponding to each bin, the TS obtained for the source, the best-fit value of the SED and its 1σ error and finally the the 95% CL upper limit.If TS > 4, we report the SED best-fit value and its 1σ error; otherwise we report the 95% CL upper limit.