© The Authors 2025

1. Introduction

Star-forming galaxies (SFGs) radiate across a wide spectrum, from radio waves to γ rays (Condon 1992; Ackermann et al. 2012; Lacki & Thompson 2013). Their nonthermal radio emission is caused by synchrotron radiation from energetic electrons in regions of intense star formation. These particles might originate from cosmic-ray (CR) accelerators, such as supernova remnants (Axford et al. 1977; Bell 1978) or young massive stellar clusters (Bykov & Fleishman 1992), which scale in number with the star formation rate (SFR) of their host galaxy. This assertion is supported by the robust correlation observed between the SFR and radio luminosity at 1.4 GHz (Condon 1992; Yun et al. 2001; Bell 2003; Kornecki et al. 2022).

Furthermore, CRs may interact with ambient protons in the interstellar medium, leading to the diffuse emission of high- and very-high-energy (VHE; E > 100 GeV) γ rays across the GeV and TeV bands (Völk et al. 1996; Blom et al. 1999; Domingo-Santamaría & Torres 2005; Persic et al. 2008; Yoast-Hull et al. 2013; Peretti et al. 2019; Krumholz et al. 2020; Roth et al. 2023). To date, the Fermi Large Area Telescope (LAT; Abdo et al. 2010a) has detected over a dozen such galaxies in the 0.1 − 100 GeV energy range (Ballet et al. 2023). This diverse group of galaxies includes starburst galaxies (SBGs) and luminous to ultraluminous infrared galaxies (ULIRGs). Additionally, SBGs such as NGC 1068 (Antonucci & Miller 1985) can host an inner active galactic nucleus (AGN) on parsec scales. These nonjetted AGNs, often found in spiral galaxies, are classified as Seyfert nuclei.

Current research generally supports the idea that emissions from SFGs in the GeV range predominantly stem from hadronic processes (e.g., Krumholz et al. 2020; Shimono et al. 2021; Kornecki et al. 2022; Roth et al. 2023). These processes may occur due to the interaction of energetic protons within the starburst nucleus on kiloparsec scales or within star formation environments of the surrounding galactic medium. The energy injected into CRs by highly active SFGs, such as Arp 220, appears to be almost entirely reprocessed into γ rays. In contrast, a considerable portion of the CR energy might be dissipated by advection or diffusion in more quiescent SFGs (Kornecki et al. 2020, 2022; Werhahn et al. 2021; Roth et al. 2021). The injection and physical processes that CRs undergo in these galaxies are expected to determine the spectral slope in the GeV−TeV range. However, the transport mechanisms and the true properties of the CR spectrum in these galaxies remain under debate, partly because of the limited number of detections at the highest energies.

Most SFGs observed at GeV energies with Fermi-LAT have not been detected in the TeV regime, except for the two well-known SBGs M82 and NGC 253. These galaxies, being the closest and most active SBGs, have been detected in the VHE range by the Very Energetic Radiation Imaging Telescope Array System (VERITAS; VERITAS Collaboration 2009) and the High Energy Stereoscopic System (H.E.S.S.; Abramowski et al. 2012; H.E.S.S. Collaboration 2018) observatories, respectively. Their emission in the TeV range is shown to essentially be an extension of the hadronic component observed in the GeV range by Fermi-LAT. This provides evidence of proton acceleration up to multi-TeV energies in extragalactic star-forming environments.

A correlation between SFR and γ-ray luminosity, measured with Fermi-LAT (Ackermann et al. 2012; Ajello et al. 2020; Kornecki et al. 2020), has been observed in SFGs, supporting the connection between CRs and the star formation process on galactic scales. However, this correlation in the GeV γ-ray range exhibits greater dispersion than the one found in the radio band at 1.4 GHz (Ackermann et al. 2012; Kornecki et al. 2022), which may arise from several factors. Some CR protons might lose their energy through nonradiative processes or escape the galaxy before radiating, leading to lower γ-ray emission levels than expected in the calorimetric limit of proton-proton processes (Martin 2014; Pfrommer et al. 2017; Peretti et al. 2019; Werhahn et al. 2021; Roth et al. 2021; Ambrosone et al. 2024). Additionally, determining the location of the γ-ray emission within the galaxy is challenging due to the limited angular resolution of GeV−TeV instruments.

Alternative studies suggest that γ-ray emission from SBGs comes from the interaction of particles accelerated within galactic superwinds, driven by the combined effects of the stellar winds of young stars and supernova explosions. The accelerated particles then cool down, typically a few kiloparsecs away from the galactic disk (Dorfi & Breitschwerdt 2012; Romero et al. 2018; Peretti et al. 2022). Seyfert nuclei can also launch powerful winds (Tombesi et al. 2010) known as ultrafast outflows, which could accelerate CRs to hundreds of PeV (Peretti et al. 2023). Investigations such as the studies of Lenain et al. (2010), Wojaczyński & Niedźwiecki (2017) and Lamastra et al. (2019) also argue that the nonthermal emission of some of these galaxies may be due to leptons accelerated as a result of AGN activity. Finally, a contribution to the diffuse γ-ray emission of SFGs is expected from nonresolved point sources such as pulsar wind nebulae (PWNs; see Vecchiotti et al. 2022), in particular in the TeV energy range (Mannheim et al. 2012; Ohm & Hinton 2013; Do et al. 2021; Chen et al. 2024). Consequently, the observed luminosity could be a combination of diffuse emission, AGNs, the population of PWNs, and emissions outside the galactic disk, such as from starburst superwinds.

The γ-ray energy range plays a pivotal role in elucidating the physical mechanisms governing the distribution of CRs within galaxies. Presently, alongside established observatories such as VERITAS, MAGIC, and H.E.S.S., a suite of new experiments is set to inaugurate a transformative era in γ-ray astronomy. Given the small number of SFGs detected at TeV energies by ground-based pointed telescopes, a priori promising instruments for establishing a larger population, owing to their wide fields of view, are the Large High Altitude Air Shower Observatory¹ (LHAASO; Cao et al. 2019) in the Northern Hemisphere and the forthcoming Southern Wide-field Gamma-ray Observatory² (SWGO; Albert et al. 2019). The sensitivity of these observatories up to a few hundred TeV could open the door to the study of CR acceleration mechanisms and γ-ray absorption effects beyond 10 TeV. The new-generation array of pointed γ-ray telescopes, the Cherenkov Telescope Array Observatory³ (CTAO; Hofmann & Zanin 2023), promises to revolutionize our understanding of SBGs with its extensive energy coverage from 20 GeV to 300 TeV, its exceptional sensitivity, and unparalleled angular resolution in the γ-ray band. The CTAO will help determine whether the observed γ rays originate from diffuse emissions, similar to those detected at radio wavelengths, or from additional contributions such as AGNs or other point sources. Thus, the CTAO has the potential to precisely delineate emission regions within some SBGs and elucidate the mechanisms governing CR transport.

The study of the most intense SFGs in the GeV range has been identified as one of the key science projects of the CTAO (Cherenkov Telescope Array Consortium 2019). Fainter γ-ray emitting galaxies, including those not yet detected by Fermi-LAT, have so far been relatively neglected despite being crucial for establishing a relevant sample of SFGs at VHE. Several studies in recent years attempt to evaluate the detectability of SFGs in the VHE range, albeit with a primary focus on the Fermi-LAT energy range (Shimono et al. 2021; Xiang et al. 2021). These former studies did not aim for a census of SFGs observable at TeV energies.

In the present work, we investigate the chances of detection of SFGs in the TeV range, considering the up-to-date sensitivity of new-generation ground-based telescopes. Our analysis incorporates recent Fermi-LAT data and integrates information from full-sky optical and infrared surveys mapping galaxies in the nearby universe. To evaluate their observability, we also considered their sky positions to select the corresponding arrays for visibility and the zenith angle to choose the most accurate response available. We worked under the assumption that the GeV emission comes mainly from proton-proton interactions and that the VHE spectrum is an extension of the emission observed in the GeV range. Our study encompassed two groups of SFGs: those previously detected with Fermi-LAT in the GeV range and a selection of galaxies that have not been detected at high energies so far. In Sect. 2, we update the study of the correlation between SFR and γ-ray luminosity for GeV-detected galaxies and search for such emission in a larger sample of SFGs. In Sect. 3, we evaluate the detectability at VHE of all the galaxies in the sample and provide a list of the best candidates to be observed with current and upcoming observatories. Our conclusions are outlined in Sect. 4.

2. Candidate star-forming galaxies in the GeV energy range

2.1. Sample of galaxies of interest

We explored the prospect of discovering new SFGs in the GeV energy range and subsequently assessed their potential detection at TeV energies. To this end, we employed the Revised Mass AssociatioN for GRavitational waves ObserVations Efficiency (MANGROVE) Sample from Biteau (2021), which represents a comprehensive collection of SFGs over more than 90% of the sky and out to 350 Mpc (z < 0.08). This near-infrared flux-limited sample comprises around 400 000 galaxies, the distances of half of which were estimated photometrically and half spectroscopically (or from the cosmic distance ladder for the closest galaxies). Stellar masses were estimated from emission in the W1 band of the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), using the mass-to-light ratio inferred for a Chabrier initial mass function (Chabrier 2003). Galaxies with strong near-infrared emission associated with an AGN were excluded from the initial MANGROVE based on their W1−W2 color (Ducoin et al. 2020). The SFR of the galaxies kept in the sample was estimated from their Hα emission for galaxies in the Local Volume (d < 11 Mpc) studied by Karachentsev et al. (2018). These authors accounted for dust attenuation through the apparent orientation of the galaxy. For other galaxies, scaling relations are used by Biteau (2021) to infer SFR from stellar mass, albeit with a dispersion ranging from 0.3 dex to 0.8 dex, depending on the morphological class of the galaxy.

Our sample was composed of galaxies with γ-ray emission from star formation that may be challenging to detect at TeV energies, either due to low SFR or a large distance from Earth. To ensure the relevance of the subsample to be examined, we first selected galaxies that are a priori bright enough to be detectable by new-generation TeV γ-ray observatories, as determined by the ratio SFR/4πd². The new-generation observatories target sources with fluxes down to ∼1 mCrab in ∼50 h, to compare with the flux of the SBG NGC 253 of about ∼3.5 mCrab at 1 TeV, as measured with previous-generation observatories in 160 hours of observations (H.E.S.S. Collaboration 2018). We thus selected all galaxies in the Revised MANGROVE sample that satisfied the relation SFR/4πd² > αSFR_NGC 253/4πd_NGC 253², where α = 1/3.5. To ensure the completeness of the sample, we applied the same selection criterion to galaxies observed in the infrared with IRAS (Sanders et al. 2003) and selected by the Fermi-LAT Collaboration (Ackermann et al. 2012), using their radio, infrared, and HCN (J = 1 − 0) luminosities as a rough tracer of the SFR. This verification, illustrated in Appendix B, confirms the greater completeness of the Revised MANGROVE sample in the Local Volume and allowed us to include NGC 1068, otherwise excluded due to its powerful Seyfert nucleus. We excluded galaxies with dwarf elliptical, lenticular, or elliptical morphology, namely the Sagittarius Dwarf Galaxy, Maffei 1, Centaurus A, and M49. After this exclusion, our sample comprised 21 galaxies, including 8 already detected by Fermi-LAT (LMC, SMC, NGC 253, M82, NGC 4945, NGC 2403, Circinus, and M31). Table 1 lists the remaining 13 galaxies that have not been yet detected by Fermi-LAT.

We collected the luminosity distance and associated uncertainty for each galaxy from the Extragalactic Distance Database⁴, based on the recent update of Cosmicflows-4 (Tully et al. 2023). For SFR values, discrepancies arise when they are derived from different spectral indicators (such as UV, Hα, [OII] lines, and far-infrared fluxes). Many of these disparities are primarily attributed to uncertainties regarding dust attenuation. To address this issue, we refined the estimation of the SFR following Kornecki et al. (2020) and employed multiwavelength tracers (in particular, FUV+IR_25 μm or Hα+IR_25 μm and total IR emission in 8−1000 μm; Kennicutt 1998; Kennicutt & Evans 2012). We took the IR_25 μm fluxes from Sanders et al. (2003) and Rice et al. (1988), FUV from Cortese et al. (2012) and Gil de Paz et al. (2007), and Hα from Kennicutt et al. (2008). For the nine galaxies with data available in both FUV and Hα, we checked that the two SFR estimates aligned within measurement uncertainties. In such cases, we chose to use the FUV values for further analysis. We verified that the refined estimates of the SFR were statistically consistent with the initial estimate from the revised MANGROVE sample and had lower uncertainties. The uncertainties in the quoted SFR values account only for the uncertainties in the distances, FUV or Hα, and IR flux measurements, and not for uncertainties in the calibration constants, nor for dispersion around the calibration relations. The estimated distances and SFRs are listed in Table 1, together with the relevant references.

2.2. Fermi-LAT data analysis for non GeV-detected galaxies

To search for high-energy γ-ray emission associated with galaxies not included in the fourth Fermi Gamma-ray LAT (4FGL) catalog, we used Fermi-LAT PASS 8 archival data⁵, spanning almost 15 years of observations from 2008 August 4 to 2023 July 12. For each galaxy, we downloaded γ-ray events in a region of interest (ROI) of 15-deg radius, centered on coordinates from SIMBAD⁶. The data in the energy range 0.1 − 100 GeV were analyzed using the publicly available Fermi Science Tools v2.2 interfaced with Fermipy v1.2.0 (Wood et al. 2017), in combination with the latest instrument response functions (IRFs; P8R3_SOURCE_V3).

We adopted a source-type event selection (evclass = 128, evtype = 3) and applied standard data quality selection criteria (“DATA_QUAL> 0 & & LAT_CONFIG= = 1”). Zenith angles were limited to 90 deg to minimize contamination from the Earth limb. We spatially binned the data with a scale of 0.1 deg per pixel and used eight logarithmically spaced bins per energy decade. A model file was built that included all sources of the 4FGL catalog Data Release 4 (4FGL-DR4; gll_psc_v33.fit, Ballet et al. 2023) that fell into the ROI, as well as the Galactic diffuse emission and extragalactic isotropic emission components. A standard binned likelihood analysis was applied in an iterative way. We fixed the spectral shape parameters of 4FGL-DR4 sources more than 10 deg away from the galaxy to account for event leakage in the ROI due to the larger point spread function at lower energies. In a second step, the sources contributing to less than 5% of the total number of counts in the ROI (N_pred/N_count < 0.05) and/or with low significance, TS < 9, had their parameters frozen. The detection significance of each source was evaluated with a test statistic (TS) defined as TS = 2 Δlog(likelihood) between models with and without the source. In the end, the only free parameters were those of sources less than 3 deg away from the target, if not frozen in the previous step, and the normalizations of the two diffuse background components. After the fit, the residuals and TS maps were obtained by subtracting the best-fit model from the data. Whether or not an excess was found in the TS map at the source location, in order to check for any detection, we refitted the data adding a new point source at the galaxy position, leaving free the photon index and normalization of the source of interest.

All the galaxies show significance considerably lower than 5σ (TS < 25) when the photon index is left free, except for M83 and M33 with TS values of 18 and 22, respectively. Notably, the excesses from M83 and M33 correspond to photon indices of 2.0 ± 0.2 and 2.3 ± 0.2, respectively, in line with the expected value for an SFG.

We calculated upper limits on the flux for each galaxy, assuming a fixed photon index of 2.2 for each source. This value of the photon index is expected in theoretical models describing the diffuse γ-ray emission from these galaxies, under the assumption that the emission mainly comes from proton-proton interactions (Yoast-Hull et al. 2013; Wang & Fields 2018; Peretti et al. 2019; Kornecki et al. 2020). We provide 95% confidence level upper limits on the energy flux integrated over the entire energy range of the analysis (${\mathcal{F}}_{0.1-100\mathrm{GeV}}^{\mathrm{UL}}$ in GeV cm⁻² s⁻¹), as shown in Table 1. We also calculated the spectral points and upper limits using four logarithmically spaced bins across the energy range.

2.3. GeV detected galaxies and the γ-ray luminosity – star formation rate correlation

After more than 15 years of surveying in the 50 MeV−1 TeV energy range, Fermi-LAT has detected more than 7000 γ-ray sources. Their properties are accessible in the 4FGL-DR4. Only a small fraction of these sources have been associated with galaxies with active star formation, either normal galaxies, SBGs, or Seyfert galaxies. We only considered galaxies that are included in the 4FGL-DR4. We did not include other possible sources that may have shown evidence of detection in the literature but that did not meet the selection criteria discussed in Section 2.1⁷. The 4FGL-DR4 provides two types of association (ASSOC1/2) for each source. The first corresponds to the name of the identified or likely associated source, and the second to a lower-confidence association. For each association, the catalog also provides an object type designation (CLASS1/2; see Table 6 in Ballet et al. 2023, for further details). As we aimed to compile a comprehensive sample of galaxies with active star formation, we included sources with class designations 1 and 2 that aligned with normal galaxies (“GAL”), starburst galaxies (“SBG”), or Seyfert galaxies (“SEY”). We found 29 γ-ray sources in the 4FGL-DR4 with at least one of these three classes.

Associations of Fermi-LAT sources are established when there is a statistically significant positional coincidence with a counterpart at another wavelength that is an a priori candidate for γ-ray emission (Abdo et al. 2010b). The automated source association carried by the Fermi-LAT Collaboration relies on curated catalogs containing potential counterparts. The classification of the same source may vary depending on the wavelength and selection criteria used in catalog generation. Because of the limited accuracy of source localization, additional data, such as temporal variability and/or spectral characteristics, are taken into account to strengthen the associations. This proves to work well for AGN sources but is trickier for quieter sources like nonjetted AGNs or SFGs.

To address this issue, we explored each candidate association. Stellar light is particularly prominent in the broad-band spectrum of SFGs (Kennicutt & Evans 2012). Such light reaches us directly from the emitting stars or is reprocessed in the infrared by dust present in the surrounding interstellar medium (ISM). Thus, SFGs show a continuum spectrum with two prominent peaks in the optical and infrared bands, respectively (da Cunha et al. 2008). This shape serves as a valuable indicator of active star formation within a galaxy. The presence of a star-forming component in the optical and infrared spectrum of the Fermi-LAT candidates was thus evaluated using the Spectral Energy Distribution (SED) Builder tool⁸.

Out of the 29 Fermi-LAT candidates, we excluded 15 galaxies that did not show clear evidence of a star-forming component at low energies. Among these were NGC 5380 and IC 678, which are classified as normal galaxies in 4FGL-DR4 but with spectra that are better represented by those of a giant elliptical galaxy (see Table A.1 in Appendix A for further details). Our study included the γ-ray source 4FGL J0737.4+6535, classified as an SBG and labeled as the infrared source WISEA J073707.21+653623.0 in the 4FGL-DR4 catalog. Bruzewski et al. (2023) associate this infrared source with a blazar. The 4FGL source is also coincident with the disk of the SBG NGC 2403 detected by Ajello et al. (2020). As a starting point, we assumed that 4FGL J0737.4+6535 is associated with the galaxy NGC 2403.

The correlation between SFR and the γ-ray luminosity (L_γ) in different energy ranges has been explored by several authors (Lacki et al. 2011; Ackermann et al. 2012; Ajello et al. 2020; Kornecki et al. 2020; Roth et al. 2021). As it is commonly accepted that the emission in the GeV range comes from proton-proton interactions, the spectrum in this energy range is expected to follow a power law with an index close to that of CRs injected in the star-forming environment. The γ-ray luminosity is expected to increase with increasing SFRs, so that the normalization of the γ-ray spectrum depends mostly on the SFR value and on the distance to the source. It should be noted that the γ-ray slope could differ from the slope of the injected CRs depending on energy due, for example, to CR escape by diffusion or to absorption. Thus, the correlation between γ-ray luminosity and SFR provides insight into the mechanisms and physical properties that govern the transport and dissipation of CRs within SFGs (Kornecki et al. 2022).

Diffusion may become dominant at sufficiently high energies (depending on the turbulence level and the size of the emission region), leading to a softening of the spectrum and inducing a transition regime in which protons can escape the galaxy before they interact and radiate. Particles may also escape due to the winds of SFGs (Veilleux et al. 2005), which can advect protons across the entire energy range, diminishing the overall normalization of the γ-ray spectrum without altering the slope (Peretti et al. 2019; Kornecki et al. 2022). This effect has been carefully studied in NGC 253 and its impact, which depends on wind speed, seems to be minor (H.E.S.S. Collaboration 2018). Nevertheless, it could be more pronounced for galaxies with stronger winds, such as M82 (Strickland & Heckman 2009).

Since SFGs can host a strong infrared field, primarily emitted by hot dust, γγ absorption may occur between these target photons and γ-ray photons. The relevance of this absorption mechanism starts to become significant at energies close to ∼10 TeV (Peretti et al. 2020) and increases at higher energies. It can be even more pronounced in galaxies with intense radiation fields, such as ULIRGs. At energies beyond 20 TeV for nearby galaxies (the most distant galaxy in our sample was Arp 220 at z ∼ 0.019), the extragalactic background light (EBL) and cosmic microwave background become the dominant target photon fields for γγ absorption (e.g., Biteau & Meyer 2022, for a review).

Here we revisited the L_γ–SFR correlation presented by Ackermann et al. (2012). Instead of using the luminosity integrated over the energy range of Fermi-LAT, we used the luminosity at a reference energy. We set the reference energy at 2 GeV as, for the γ-ray sources in our sample, this was close to the pivot energy of the best-fit power law, where the uncertainty on the γ-ray flux was the smallest. We obtained the 2 GeV energy flux and associated uncertainty using the best-fit power-law parameters from the 4FGL-DR4. For NGC 2146, Arp 220, and Arp 299, luminosity distances not listed in Cosmicflows-4 (Tully et al. 2023) were adopted from Kornecki et al. (2020). Except for NGC 7059, SFRs were recalculated following Kornecki et al. (2020) and references therein, using updated distances. For NGC 7059 we took the IR_25 μm flux from Rice et al. (1988) and the FUV flux from Bouquin et al. (2018) to estimate the SFR. We report these values in Table 2.

In the upper panel of Fig. 1, we show the γ-ray luminosity of SFGs at 2 GeV, L_{γ,  2 GeV}, as a function of their SFR, Ṁ_∗. We fitted the data in Table 2 with the orthogonal regression method to include the errors in both L_γ and SFR, using the Orthogonal Distance Regression Package (ODRPACK)⁹. We also show in Fig. 1 the upper limits obtained in Section 2.1, which we did not include in the fit¹⁰. A power-law fit, $L_{\gamma ,2\mathrm{GeV}}=C{\dot{M}}_{\ast }^m$, where L_{γ,  2 GeV} is in erg s⁻¹ and Ṁ_∗ is in M_⊙ yr⁻¹, to all the data (pink and blue dots in Fig. 1) yields a slope value of m = 1.27 ± 0.10 and a constant log C = 38.28 ± 0.10, with a residual variance of 25. The residual variance is calculated as the average squared orthogonal distance from the data points to the fitted model, accounting for errors along both axes and normalized by the associated number of degrees of freedom. We show this fit as a dotted green line in Fig. 1. Most galaxies follow a clear trend, except NGC 3424, NGC 4945, Circinus, NGC 2403, and NGC 7059.

Fig. 1. Upper panel: Lγ,  2 GeV – SFR observed correlation. The dotted green line corresponds to the best-fit model for all the data points of detected SFGs and the shaded region represents the 68% confidence level band. The dashed pink line and associated shaded region show the best-fit model excluding blue markers. The upper limits obtained in Section 2.1, shown as downward-pointing triangles, are not taken into account in the fit. Lower panel: Residuals to the fit excluding blue markers normalized to their standard deviation. The dashed black lines correspond to 4 standard deviations.

If we exclude these five galaxies from the sample, the correlation is tighter, with m = 1.31 ± 0.08, log C = 38.19 ± 0.08, and a smaller residual variance of 14. We show the corresponding fit as a dashed pink line in Fig. 1. The lower panel of Fig. 1 shows the pulls, i.e., residuals for this fit normalized to the standard deviation. The figure shows that NGC 3424, NGC 4945, Circinus, NGC 2403, and NGC 7059 stand out by more than 4σ from the best fit (at 19, 14, 8.7, 6.8, and 22σ, respectively).

The deviation of these galaxies from the correlation followed by the other SFGs could potentially stem from an additional contribution associated with an AGN, as was discussed by Gavazzi et al. (2011) for NGC 3424 and NGC 2403, Lenc & Tingay (2009) for NGC 4945, and Prieto et al. (2004) for Circinus. Such a deviation could also be attributed to a misestimation of their SFRs (see discussion in Kornecki et al. 2020). The emission from NGC 2403 could also be contaminated by a supernova explosion, as discussed by Xi et al. (2020a), or it could just be misassociated and, as suggested by Bruzewski et al. (2023), be a blazar. Finally, the γ-ray source 4FGL J2127.6−5959, associated with the galaxy NGC 7059 in the 4FGL-DR4, is the one that deviates the most from the correlation. Interestingly, the γ-ray source shows a power-law index of 1.8 ± 0.1, which is marginally outside the expected and observed values for the rest of the sample.

3. Candidate star-forming galaxies in the TeV energy range

3.1. TeV gamma-ray spectra and sensitivity

The best-fit model for the correlation discussed in Section 2.3 provides an estimate of the γ-ray luminosity of an SFG at 2 GeV, based on its SFR. Using the correlation and the expected index value of 2.2, we developed a simple empirical model to estimate the flux of SFGs. We assumed a power-law SED with a constant slope of α = 2.2 from GeV to TeV energies and a γ-ray luminosity provided by the best-fit model to the pink markers in Fig. 1. It should be noted that the γ-ray emission from some of these sources could deviate from a power law. This is the case, for example, for M31, NGC 253, and M82, the γ-ray spectra of which are better modeled by a log-parabola in the 4FGL. However, the curvature significance is not substantial (approximately 3σ). Therefore, the power-law hypothesis remains appropriate. We used the EBL model developed by Domínguez et al. (2011) to account for γγ absorption on megaparsec scales in the intergalactic medium. The γ-ray energy spectrum was thus estimated as

$$\begin{array}{c}\hfill \mathcal{F}[\mathrm{erg}{\mathrm{s}}^{-1}c{m}^{-2}]=\frac{L_{2\mathrm{GeV}}(\mathrm{SFR})}{4\pi d^2}\times {\left(\frac{E}{2\mathrm{GeV}}\right)}^{2-\alpha }{\mathrm{e}}^{-\tau (E,z)},\end{array}$$(1)

where τ(E, z) is the optical depth from the model of Domínguez et al. (2011). To compute z from the luminosity distance, d, we considered a Hubble constant of H₀ = 67.4 km s⁻¹ Mpc⁻¹ and a matter density of Ω_m = 0.315 (Planck Collaboration VI 2020).

We also evaluated an even simpler model for the sources detected at GeV energies, namely the extrapolation of their spectrum measured by Fermi-LAT with EBL attenuation at the highest energies. For SFGs not included in the 4FGL, for which we did not have an observed spectrum, we provided an upper bound on the energy flux at 2 GeV, ${\mathcal{F}}_{2\mathrm{GeV}}^{\mathrm{UL}}$, as estimated from the upper limits on the integrated energy flux, ${\mathcal{F}}_{0.1-100\mathrm{GeV}}^{\mathrm{UL}}$, derived from Fermi-LAT data, assuming a photon index of 2.2.

The relations in Eq. (1) or other extrapolations to the highest energies are strictly valid in the energy range where internal absorption is subdominant. As discussed in the previous subsection, internal absorption effects are expected to occur at energies on the order of tens of TeV and to depend on the intensity of the galaxy’s infrared field. We did not consider the modeling of internal absorption within the galaxy, as we aimed to maintain a simple, geometry-independent model in our study; however, we discuss the implications this may have when drawing conclusions. Spectral effects induced by CR escape could be noticeable at energies close to 10 TeV, but they would not completely change the spectral slope, except in the extreme case of fast diffusion as in Milky Way-like galaxies (Strong et al. 2011; Do et al. 2021). Therefore, the results from our simple model should be considered optimistic beyond 10 TeV.

As mentioned in the introduction, the best new-generation telescopes for the detection and study of SFGs at TeV energies are expected to be the CTAO, LHAASO, and SWGO. The CTAO, which is currently under construction, will be situated at the Paranal Observatory in Chile and at the Instituto de Astrofísica de Canarias, in La Palma, Spain (Cherenkov Telescope Array Consortium 2019). Using Gammapy¹¹ (Donath et al. 2023), we estimated the point-source sensitivity for 50 h of CTAO observation assuming 20 logarithmic energy bins across the energy range 0.03 − 30 TeV. The significance in each bin was computed for a 1D analysis (ON-OFF regions) assuming a fixed source offset of 0.5° from the pointing position. As we investigated which galaxies may approach the detection threshold of the CTAO, we imposed a minimal number of 5 signal counts per bin and a minimal significance of 3σ per bin, for a background fraction of 0.1. We verify with dedicated simulations that comparing the spectra of SFGs with these differential sensitivity curves provides a good estimate of detectability: for a photon index of 2.2, a spectrum above the sensitivity curve corresponds to an expected significance of the source, combined over all energy bins, above 5σ. We computed the sensitivity for an azimuth-averaged pointing using the latest version of the CTAO instrument response functions (IRFs), prod5 v0.1¹². The IRFs of the northern and southern arrays (CTAO-N and CTAO-S, respectively) were calculated for three different zenith angles (zen): 20, 40, and 60°. We determined the number of hours of darkness for each array using the TeVCat Object Visibility Tool¹³ and kept the IRF with the smallest zenith angle for which more than 50 hours of observation were available over a year for each galaxy. We discarded the low energy points on the sensitivity curve where the IRF was not well defined (implied energy thresholds: CTAO-S, zen = 40 deg, E > 0.05 TeV; CTAO-S, zen = 60 deg, E > 0.2 TeV; CTAO-N, zen = 60 deg, E > 0.07 TeV).

LHAASO is a multicomponent facility located in Daocheng, China that serves as a γ-ray telescope operating in the energy range from 0.1 TeV to 1 PeV. Due to its geographical location, most of the LHAASO field of view is in the Northern Hemisphere, with access to sources of declination larger than −30° (Cao et al. 2019; Yang & Razzaque 2019). Therefore, we compared the expected γ-ray flux of SFGs with the LHAASO sensitivity (with 20″ photomultiplier tubes, LHAASO Collaboration 2021) in cases of galaxies with such declinations. SWGO (Albert et al. 2019) is a forthcoming γ-ray observatory in South America, featuring ground-level particle detection technology with a wide field of view spanning several steradians. Focused on energies from hundreds of GeV to the PeV scale, SWGO will primarily use water Cherenkov detectors. The observatory is planned for construction between 10 and 30° South latitude to ensure a clear view of the Galactic Center. We then compared the SWGO optimistic sensitivity to sources with negative declination. This sensitivity was estimated for a steady point source at a zenith angle of 20 deg and under the assumption that it can be observed for 6 hours per day (Albert et al. 2019). The differential sensitivity curves of LHAASO and SWGO typically allow for 5 to 10 bins per decade in energy. To claim a detection, the number of excess events (over the background) is typically required to be above ten events, and the significance threshold for detection is generally set at 5σ.

As the calculation of the CTAO’s sensitivity to extended sources is more complex than for point sources and falls beyond the scope of this study, five galaxies in the Local Group deserve a specific discussion due to their proximity, namely LMC, SMC, NGC 6822, M31, and M33. Assuming that most of the emission is located near the central molecular zone of radius R < 1 kpc, the star-forming nucleus would appear “point-like” if it were viewed within a cone of half-opening angle θ ≈ R/d < θ₆₈(0.2 TeV), where θ₆₈(0.2 TeV)≈6′ is the 68% containment radius at 0.2 TeV for the CTAO¹⁴. Point-like emission would then be expected for SFGs at d > R/θ₆₈ ≈ 600 kpc × (R/1 kpc)×(θ₆₈/6′)⁻¹. With distances of 740 kpc and 850 kpc, M31 and M33 could potentially be detected as point-like sources and so are considered in this work for comparison with the sensitivity of new-generation γ-ray observatories. In contrast, the central molecular zones of the SMC and LMC, with distances of 62 kpc and 50 kpc, respectively, could appear as extended sources and require dedicated studies (e.g., Acharyya et al. 2023). Therefore, we omit direct comparison of the SMC and LMC spectra with the point-source sensitivity of γ-ray observatories. At a distance of 460 kpc, NGC 6822 could show a mild extension. We assumed that the source was point-like and will show that even in this optimistic scenario the source is difficult to detect.

3.2. Results: Detectability in the TeV energy range

By comparing the CTAO point-source sensitivity with the extrapolated best-fit power-law spectrum from 4FGL-DR4, or the upper limits for SFGs not included in 4FGL, we can identify the most promising candidate point sources at TeV energies. In Fig. 2, we show the SEDs of the SFGs detected at GeV energies that are the best candidates for detection in the TeV range by the CTAO, LHAASO, and SWGO. The extrapolated γ-ray spectra are above the CTAO point-source sensitivity for NGC 253, M82, NGC 1068, NGC 4945, and Circinus. Figure 3 shows the interesting candidates for which GeV γ-ray observations could still allow for TeV detection. The next best candidates are shown in Fig. 4. The SEDs of the other SFGs that were not selected as strong candidates for detection at TeV energies are provided in Appendix A.

Fig. 2. SEDs of the top candidates selected from the GeV-detected sample of SFGs. In each panel, pink dots denote spectral points from Fermi-LAT, the gray line and gray shaded region represent the best power-law fit to the γ-ray data, as provided in the 4FGL-DR4 catalog, including absorption on the EBL at multi-TeV energies. The orange dotted line and associated shaded region depict the empirical model scaled to the SFR of the galaxy, along with its associated uncertainty. The line defined by the triangles indicates the 3σ-level point-source sensitivity of the CTAO in 50 h. The dashed gray line and dotted-dashed black line show the LHAASO and SWGO sensitivities, respectively. VHE data and upper limits for NGC 253, M82, and NGC 1068 were taken from H.E.S.S. Collaboration (2018), VERITAS Collaboration (2009) and Acciari et al. (2019), respectively.

Fig. 3. SEDs of interesting candidates detected at GeV energies, namely M31 (upper panel), or on the verge of GeV detection, namely M33 (mid panel) and M83 (lower panel). The line and color codes match those of Fig. 2. The solid black line with arrows represents the upper bound on the expected γ-ray emission of each galaxy including the absorption on the EBL at multi-TeV energies.

Fig. 4. SEDs of possibly interesting candidates from the GeV detected sample. The line and color codes match those of Fig. 2.

Among the SEDs of GeV-detected SFGs listed in Table 2, the empirical model scaled with SFR is consistent with the extrapolated γ-ray spectrum from Fermi-LAT for all galaxies but NGC 4945, NGC 3424, NGC 2403, and M31. Such a discrepancy was expected for the first three galaxies, as their γ-ray luminosity expected from the correlation in Fig. 1 differs by more than 4σ from the observation. Regarding M31, the discrepancy is due to the fact that its observed photon index (α ≈ 2.5) is softer than that assumed in the empirical model scaled to SFR. Normal SFGs (with low SFRs) exhibit less dense environments, potentially reducing the effective density of ambient gas encountered by CR protons and facilitating their escape (e.g., Kornecki et al. 2020). This context can enhance the relative contributions of leptonic emission and unresolved point-source populations to the integrated γ-ray emission (Do et al. 2021). Consequently, the final spectrum observed by Fermi-LAT may deviate from the slope that our simple model predicts.

For SFGs not detected in the GeV energy range, the upper-bound spectrum is shown as a black line with downward arrows in Fig. 3 and as colored lines in the Appendix (Fig. A.1). As shown in Fig. 3, the flux of M83 and M33 inferred from the scaling of the γ-ray luminosity with SFR is at the threshold of the sensitivity of the CTAO, while the upper limit inferred from GeV γ-ray observations could still allow for TeV detection provided there is sufficient exposure. We did not include NGC 7059 in the best-candidate sample, although the extrapolation of the spectrum of the associated Fermi-LAT source is above the CTAO sensitivity curve. As illustrated in Fig. 4, the empirical model scaled with SFR differs strongly from the extrapolated γ-ray spectrum for NGC 7059, a discrepancy which is amplified at the highest energies by the hard spectrum of the γ-ray source. This discrepancy suggests that the dominant γ-ray emission does not originate from the star formation in the disk of this galaxy (Kaur et al. 2019).

We discuss the sources in more detail on a case-by-case basis in the next section. A summary of the prospects for TeV detection is also given in Table B.1.

3.3. Case-by-case discussion

3.3.1. NGC 253 and M82 (see Fig. 2):

In the TeV energy range, NGC 253 and M82 show empirical model spectra and extrapolated spectra from Fermi-LAT data that are above the 50 hsensitivity of the CTAO over almost one and a half orders of magnitude in energy. These two SBGs have already undergone dedicated investigations by the CTAO Consortium (Cherenkov Telescope Array Consortium 2019). Furthermore, the extrapolated model for NGC 253 and both models for M82 are above the sensitivity of LHAASO beyond 30 TeV. The extrapolated and empirical models for NGC 253 are also above the sensitivity of SWGO. For these two observatories, the detectability of NGC 253 and M82 strongly depends on the level of internal absorption within the sources (see Peretti et al. 2020, for M82).

3.3.2. NGC 1068, NGC 4945, and Circinus (see Fig. 2):

Albeit with a narrower margin, the extrapolated spectra of NGC 1068, NGC 4945, and Circinus are above the CTAO sensitivity curve. The extrapolated models are consistent with the empirical models for NGC 1068¹⁵, and Circinus. This is not the case for NGC 4945: its empirical model falls below from the extrapolation from GeV observations. As mentioned above, this is not surprising as this galaxy is an outlier from the best-fit correlation in Fig. 1. A higher SFR or an additional nonthermal contribution from either an AGN or a superwind (Pérez-Beaupuits et al. 2011) could bring the two models into agreement. Previous studies excluded NGC 4945 and NGC 1068 from consideration due to the potential presence of AGNs (Shimono et al. 2021). Without clear evidence that the GeV γ-ray emission originates from the AGN, we studied these galaxies and highlight NGC 1068 and NGC 4945 as potential targets for CTAO observations. The extrapolated spectra of NGC 4945 and the Circinus galaxy are above the sensitivity of SWGO for energies greater than 10 TeV. As for NGC 253 and M82, observations above 10 TeV, together with those at TeV energies, will lead to constraints on the amount of internal absorption.

3.3.3. M83 and M33 (see Fig. 3):

These galaxies are marginally detected by Fermi-LAT with photon indices of 2.1 ± 0.2 (TS = 18) and 2.4 ± 0.2 (TS = 22), respectively. Previous excesses near the position of M83 are reported by Xing & Wang (2023) and Ambrosone et al. (2024), with TS values of 14 between 0.1 − 500 GeV and 15 between 1 − 1000 GeV, respectively. Xi et al. (2020b) and Ajello et al. (2020) claim the detection of M33 (photon index of 2.41 ± 0.16 in Ajello et al. 2020). Xi et al. (2020b) argue that the γ-ray excess from M33 does not originate from the center but rather from a region in the northeast that is positionally coincident with a supergiant HII region. Nevertheless, this source remains absent from the 4FGL-DR4 catalog. Further data from Fermi-LAT will enable a more accurate determination of the GeV spectrum of M33 and of its extrapolation to higher energies.

The empirical model and the upper bound on the extrapolated γ-ray emission are in line for M83 and M33. The spectral models for M83 and M33 are close to the sensitivity of the CTAO. In line with the predictions of Shimono et al. (2021), we identify M33 and M83 as interesting candidates for observation in the TeV energy range.

3.3.4. M31 (see Fig. 3):

The extrapolated spectrum of M31 falls below the 50h sensitivity of the CTAO. Despite being cataloged as a point-like source in the 4FGL catalog, M31 (also known as 4FGL J0043.2+4114) has been suggested to exhibit marginal γ-ray extension by Abdo et al. (2010a), Ackermann et al. (2017) and Ajello et al. (2020). Depending on the spatiospectral model employed, the best-fit photon index varies considerably, from 2.1 ± 0.3 to 2.5 ± 0.1 (Abdo et al. 2010a; Ballet et al. 2023). McDaniel et al. (2019) propose that M31 γ-ray emission is a mix of pionic and Compton scattering from primary and secondary electrons on ambient radiation fields. Alternatively, Xing et al. (2023) suggest it may arise from two different point sources: one at the galaxy center, likely from unresolved objects such as millisecond pulsars (Eckner et al. 2018), and another 0.4° away, with an unclear origin. Additionally, Persic et al. (2024) propose a more complex scenario, involving a combination of diffuse pionic, pulsar, and nuclear black hole emissions in the M31 core. Thus, M31 remains an exceptional candidate for exploring the contribution of additional source populations and is also crucial for studying the L_γ–SFR correlation, as it is one of the few candidates with a low SFR.

3.3.5. NGC 2146, NGC 2403, NGC 3424 (see Fig. 4):

The extrapolated spectrum is consistent with the empirical model only for NGC 2146, while a tension is observed for NGC 2403 and NGC 3424, suggesting that an additional emission may be needed (see Sect. 2.3). The extrapolated spectra of these three sources are below the 50 h sensitivity of the CTAO, although NGC 2403 could be marginally detectable with the CTAO, and with LHAASO if not prevented by internal absorption. With updated Fermi-LAT data, contrary to what Cao et al. (2019) anticipate, we do not find sufficient emission from NGC 2146 at the highest energies for detectability by LHAASO in 1 yr.

3.3.6. Arp 299 and Arp 220 (see Fig. 4):

The empirical model closely matches the Fermi-LAT extrapolation for these galaxies. Both models are below the 50 h CTAO sensitivity and the 1 year sensitivity of LHAASO. While Arp 220 has been proposed as part of the key science project of the CTAO Consortium dedicated to star formation, other galaxies such as Circinus, NGC 1068, and NGC 4945 (and even NGC 2403), emerge as more promising candidates for detection. The case of Arp 299 bears similarities to that of Arp 220.

Although they are not among the preferred candidates, it is crucial to keep an eye on these luminous infrared galaxies for further investigation, as improved measurements of their GeV fluxes may grant a TeV flux close to the sensitivity limit of the CTAO. Incorporating at least one of these objects in the observation program of the CTAO is pivotal for exploring a possible correlation between TeV emissions and SFR (Cherenkov Telescope Array Consortium 2019; Kornecki et al. 2023), as at least one object per decade at SFR is necessary to comprehensively assess the influence of CRs on SFGs. This could involve extending observation time, for example, to 100 hours, and conducting a more exhaustive analysis of theirγ-ray spectra, leveraging their multiwavelength emissions asconstraints.

3.3.7. NGC 7059 (see Fig. 4):

The source 4FGL J2127.6−5959 is discussed by Foschini et al. (2022), who question its association with the SFG NGC 7059 due to inconsistent sky locations and counterparts observed at radio and X-ray wavelengths. As mentioned previously, NGC 7059 not only shows an unexpected normalization with respect to the L_{γ, 2 GeV}–SFR correlation, but also exhibits a hard photon index more aligned with blazar emission (Kaur et al. 2019). Nevertheless, the association of this source with NGC 7059 cannot be ruled out, necessitating further studies to evaluate its true counterpart (Foschini et al. 2022). Assuming absorption on the EBL for a distance consistent with that of NGC 7059, the model extrapolated from Fermi-LAT data exceeds the sensitivity curves of both the CTAO and SWGO. Observing this source at VHE would be crucial for measuring its cutoff energy, providing key insight into the association of 4FGL J2127.6−5959.

3.3.8. M106, M51, NGC 1569 and NGC 6822 (see Fig. A.1):

These galaxies show one energy bin with TS > 4 in in the GeV energy range. This could suggest a simple and expected background fluctuation or, alternatively, it might be a remote indication of an excess linked to the galaxy emission. Theaforementioned scenario is more likely to apply to M106 and possibly M51 (TS = 13 and 7, respectively), due to their high Galactic latitude. The excess associated with M106 is dominated by low-energy photons (best-fit photon index close to a value of five), increasing the likelihood that these photons originate from nearby regions rather than being directly connected to emission from the galaxy. Although M51 is identified as an interesting candidate at GeV energies by Rojas-Bravo & Araya (2016), the empirical model of its VHE flux and the extrapolated upper limit do not favor its observation in the TeV range.

The extrapolated upper limit of NGC 6822 (TS = 8) is close to the CTAO detection limit. This galaxy, as NGC 1569 (TS = 10), is predicted to have a low γ-ray flux according to its SFR. Neither of these two galaxies appears to be a strong candidate for VHE detection.

3.3.9. NGC 3621, NGC 55, M81, NGC 6946, M101, IC 342, and NGC 300 (see Fig. A.1):

For these galaxies, the extrapolated upper limit and the empirical models are consistent, with the exception of NGC 6946, where the empirical spectrum prediction is marginally above the upper bound. This discrepancy between the two models for NGC 6946 may be attributed to its location near the Galactic plane (b ≈ 11.7 deg), where the Fermi-LAT background model is less accurate, or to the relatively large uncertainty regarding its SFR estimate. In all cases, the spectra are below the sensitivity curves of TeV observatories. Earlier studies of SFGs (Shimono et al. 2021) suggest that IC 342 and NGC 6946 are likely to be detected at VHE. This is at variance with our conclusions, which highlights the importance of our updated study of Fermi-LAT data and of our comparison with the appropriate instrument response functions at higher energies.

4. Conclusion

Star-forming galaxies are of significant interest due to their association with young stars and supernovae, which enable particle acceleration and CR production. Their high target density and strong magnetic fields foster CR interactions, resulting in nonthermal radiation, including γ-rays. However, current instruments have not clearly unveiled the origin of the γ-rays from SFGs. In this paper, we explored SFGs detectable at VHE, using the latest Fermi-LAT data and the expected sensitivity of current and future VHE observatories. We used updated data to gather distances, SFRs, and γ-ray fluxes or upper limits for 27 SFGs. Using the 14 SFGs detected in the 4FGL-DR4, we present a revised L_{γ, 2 GeV}–SFR correlation. We employed a simple empirical emission model to predict the γ-ray flux of SFGs based on this correlation. This model assumes diffuse emission from proton-proton interactions, with a photon index of 2.2 from GeV to TeV energies, and a γ-ray luminosity proportional to the SFR. Additionally, we constrained the γ-ray emission of these galaxies using a simple extrapolation of the Fermi-LAT data, which allowed us to constrain their TeV emission more precisely. In both cases, we accounted for the EBL absorption at the highest energies, but did not account for internal absorption effects, as the latter depend on the specific geometry of both the CR emission region and the target photon field.

We obtain promising results regarding the prospects for detection with the CTAO. Among the 14 GeV-detected SFGs, we identify 3 new candidates, NGC 1068, NGC 4945, and Circinus – in addition to M82, NGC 253 and the Magellanic Clouds – that we recommend for observation in the VHE domain. Notably, two galaxies that are not currently cataloged at GeV energies, M83 and M33, may also be within the CTAO detection range, highlighting their potential as targets for further investigation.

The case of M31 is more open in that, while its SFR suggests detectability at VHE, the extrapolation of its uncertain spectrum at lower energies does not. We find that galaxies such as Arp 220, Arp 299, NGC 3424, and NGC 2403 may be close to the detection limit of the CTAO. This suggests that, with an exposure time exceeding 50 hours, these galaxies have the potential to be detected in the TeV range.

With a similar aim, Shimono et al. (2021) predict the γ-ray emission properties of nearby galaxies based on the model of Sudoh et al. (2018) and discuss the prospect for detectability with the CTAO. With respect to this work, we confirm the relevance of NGC 253, M82, M83 (aka NGC 5236), and M33 but discard NGC 6946 and IC 342 as primary targets. We also highlight three new candidates suggested for observation with the CTAO (NGC 1068, NGC 4945, and the Circinus galaxy) and list four other galaxies that are near the CTAO sensitivity limit (NGC 2403, NGC 3424, Arp 220, and Arp 299).

This work is also one of the first to conduct a straightforward population analysis of SFG candidates for study by SWGO and LHAASO. NGC 253 and M82 have the best potential for detection by LHAASO, and NGC 2403 may not be too far from the detection limit, provided their GeV flux can be extrapolated to the highest energies. NGC 253, NGC 4945, and the Circinus galaxy could potentially be detected by the upcoming SWGO experiment. However, these galaxies may remain undetectable due to the effects of internal absorption at energies greater than 10 TeV. A detailed model of their radiation field intensity and spatial distribution, along with VHE observations, might offer significant insights into constraining the impact of internal absorption on their spectrum. We recommend detailed predictions of VHE spectra for the identified candidates to fully assess their detection potential. Additionally, improved measurements of the GeV spectra of M31, M33, M83, Arp 220, and Arp 299 would provide a better understanding of their emission and a sharper view of the prospects for detection at VHE.

In a multimessenger context, investigating nearby SFGs with forthcoming observations presents a unique opportunity to constrain their neutrino flux (Ambrosone et al. 2021). Comprehensive research on the individual neutrino emissions of SFGs and their collective impact on the neutrino background could shed light on the origin of TeV neutrinos from NGC 1068 (IceCube Collaboration 2022). Upcoming neutrino telescopes will potentially detect other point-like sources colocated with nearby SFGs, linking emissions to observed star-forming or black hole activity.

In conclusion, our study offers a new perspective on potential VHE sources among SFGs, paving the way for future observations and advancements in understanding high-energy astrophysical phenomena in these galaxies. Soon, the CTAO will be ready to investigate the production and propagation of CRs in various sources. The preparation of both theoretical models and a large sample of γ-ray emitting SFGs in the GeV–TeV range is essential to test these predictions. Having identified the most promising candidates, further investigation of these galaxies is warranted. More precise predictions of their spectra, achieved through CR transport modeling and multiwavelength spectral fitting, would enhance our understanding and predictive capabilities for their TeV emissions.

Acknowledgments

This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, the NASA’s Astrophysics Data System (ADS) and also the CTAO instrument response functions provided by the CTAO Consortium and Observatory, see https://www.ctao-observatory.org/science/cta-performance/ (version prod5 v0.1; Observatory & Consortium 2021, for more details). The authors would like to thank the reviewer, whose comments substantially improved the quality of this manuscript. PK acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033. JB gratefully acknowledges funding from ANR via the grant MICRO, ANR-20-CE92-0052.

References

Appendix A: Complementary tables and plots

In this appendix, we present additional information that complement the findings discussed in the main body of the paper. Table A.1 lists the Fermi-LAT sources catalog in the 4FGL-DR4 that could be associated with SFGs but were not identified as such in multiwavelength investigation. Figure A.1 provides upper limits on the spectra of candidate SFGs that are deemed of secondary interest for VHE observation.

Fig. A.1. SEDs of nonselected galaxies from the Southern (left panel) and Northern (right panel) hemispheres. The sensitivity lines shown in black match those in Fig. 3. The spectral points and upper limits are in different colors for each galaxy, as labeled in the figure.

Appendix B: Summary table

Table B.1 gives an overview of the SFGs studied here and the level of interest for the new γ-ray observatories. The first four columns provide information about the identification and location of the SFGs. The next three columns indicate the samples for which the SFGs verify the selection criteria defined in section 2.1. It should be noted that the five most distant sources do not satisfy these criteria for any sample, although the galaxies are detected by Fermi-LAT. The next three columns give the SFR and the constraints on the γ-ray flux. The last two columns qualify the level of interest of the observation of each SFG by the CTAO, LHAASO, and SWGO, thus summarizing the discussion provided in section 3.3.

All Tables

All Figures

Fig. 1. Upper panel: Lγ,  2 GeV – SFR observed correlation. The dotted green line corresponds to the best-fit model for all the data points of detected SFGs and the shaded region represents the 68% confidence level band. The dashed pink line and associated shaded region show the best-fit model excluding blue markers. The upper limits obtained in Section 2.1, shown as downward-pointing triangles, are not taken into account in the fit. Lower panel: Residuals to the fit excluding blue markers normalized to their standard deviation. The dashed black lines correspond to 4 standard deviations.

In the text

Fig. 2. SEDs of the top candidates selected from the GeV-detected sample of SFGs. In each panel, pink dots denote spectral points from Fermi-LAT, the gray line and gray shaded region represent the best power-law fit to the γ-ray data, as provided in the 4FGL-DR4 catalog, including absorption on the EBL at multi-TeV energies. The orange dotted line and associated shaded region depict the empirical model scaled to the SFR of the galaxy, along with its associated uncertainty. The line defined by the triangles indicates the 3σ-level point-source sensitivity of the CTAO in 50 h. The dashed gray line and dotted-dashed black line show the LHAASO and SWGO sensitivities, respectively. VHE data and upper limits for NGC 253, M82, and NGC 1068 were taken from H.E.S.S. Collaboration (2018), VERITAS Collaboration (2009) and Acciari et al. (2019), respectively.

In the text

	Fig. 3. SEDs of interesting candidates detected at GeV energies, namely M31 (upper panel), or on the verge of GeV detection, namely M33 (mid panel) and M83 (lower panel). The line and color codes match those of Fig. 2. The solid black line with arrows represents the upper bound on the expected γ-ray emission of each galaxy including the absorption on the EBL at multi-TeV energies.
In the text

	Fig. 4. SEDs of possibly interesting candidates from the GeV detected sample. The line and color codes match those of Fig. 2.
In the text

	Fig. A.1. SEDs of nonselected galaxies from the Southern (left panel) and Northern (right panel) hemispheres. The sensitivity lines shown in black match those in Fig. 3. The spectral points and upper limits are in different colors for each galaxy, as labeled in the figure.
In the text