Fermi-LAT detection of A new starburst galaxy candidate: IRAS 13052-5711

For this analysis, the region of interest (ROI) was a 20^∘.0 $\times$ 20^∘.0 square centered at (R.A. = 197^∘.078, Decl. = -57^∘.458; labeled as 4FGL J1308.9-5730 in 4FGL). The energy range of the photon events was 100 MeV to 500 GeV; the observation period from 2008-08-04 to 2023-01-12 was selected to perform this analysis. Moreover, we selected the Pass 8 data with the “Source” event class (evtype = 3 & evclass = 128) and excluded events with zenith angles above 90^∘ to minimize the contamination from the Earth Limb. Version v11r5p3 package¹¹1https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/from Fermi Science Tools (Fermitools) was selected for this analysis. For the instrumental response function, we adopted “P8R3_SOURCE_V3”. The data analysis thread provided by the Fermi Science Support Center²²2http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/ was strictly followed. All the data were divided into 37 logarithmic energy bins, and each spatial pixel size was set to be 0^∘.1 $\times$ 0^∘.1. Using the binned maximum likelihood method, we selected all known sources of the ROI listed in 4FGL and the two diffuse background templates, which are the galactic diffuse emission (gll_iem_v07.fits) and the isotropic extragalactic emission (iso_P8R3_SOURCE_V3_v1.txt)³³3http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html, for this analysis (Abdollahi et al., 2020). All spectral parameters of sources within 5^∘.0 from the center of the ROI and the normalizations of the two diffuse backgrounds were left free.

First, we investigated the location of IRAS 13052-5711 is (R.A., decl. = 197^∘.078, -57^∘.458) utilizing SIMBAD⁴⁴4http://simbad.u-strasbg.fr/simbad/. We then investigated the $\gamma$-ray emission at the location of IRAS 13052-5711 in the 4FGL and found that the distance between IRAS 13052-5711 and 4FGL J1308.9-5730 is small, with a value of 0^∘.094. Using the gttsmap command, we generated a 5^∘.0 $\times$ 5^∘.0 test statistic (TS) map of IRAS 13052-5711 in the 0.5-500 GeV band, which was selected to avoid an excessively large point spread function, as shown in the left panel of Figure 1. We found significant $\gamma$-ray radiation with TS value=43.94 and 6.63$\sigma$ in the center and no significant $\gamma$-ray residual radiation in the surrounding region.

To check whether the $\gamma$-ray radiation originates from IRAS 13052-5711, we used the SIMBAD position of IRAS 13052-5711 to replace the position of 4FGL J1308.9-5730 and used the Fermitools to deduct the $\gamma$-ray radiation at this position. We found that the $\gamma$-ray emission at this location can be completely subtracted, as shown in the middle panel of Figure 1, which indicates that this $\gamma$-ray emission likely originates from IRAS 13052-5711.

Subsequently, the gtfindsrc command was used to generate the best-fit position of the target object (R.A., decl. = 197^∘.150, -57^∘.495), and its 1$\sigma$ and 2$\sigma$ error circles were 0^∘.036 and 0^∘.059, respectively. We found that the SIMBAD position of IRAS 13052-5711 was within a 2$\sigma$ error circle, indicating that they are related in spatial position. We used the best-fit position to replace the position of 4FGL J1308.9-5730 in all subsequent analyses. Other possible $\gamma$-ray candidates were also investigated using SIMBAD, NASA/IPAC Extragalactic Database (NED)⁵⁵5https://ned.ipac.caltech.edu/ and Aladin⁶⁶6http://aladin.cds.unistra.fr/aladin.gml. However, except for IRAS 13052-5711, we did not find other possible $\gamma$-ray candidates, like AGNs, $\gamma$-ray bursts, supernova remnants, pulsars or other galaxies. Therefore, we deduced that IRAS 13052-5711 remains the most likely candidate for the detected $\gamma$-ray emission at this location.

We tested the $\gamma$-ray spatial distribution of IRAS 13052-5711 in the 1-500 GeV energy band, which was selected to avoid a large point spread function in the low energy range. We chose two-dimensional Gaussian and uniform disk templates to separately test the target source. The independent variables $\sigma$ and radius of the two spatial templates ranged from 0^∘.1 to 1^∘.0, and the increment was set to 0^∘.05. We found that the TS value of the target source showed a gradually decreasing trend with the increase of $\sigma$ and radius. When the values of $\sigma$ and radius were 0^∘.1, the TS value of the target source reached its maximum. We then calculated $\rm TS_{ext}=2log(L_{ext}/L_{ps})$ (Lande et al., 2012) and found that the $\rm TS_{ext}\approx 8$ for the two models. Because the $\rm TS_{ext}<$16, the $\gamma$-ray emission of IRAS 13052-5711 had no significant extension features (Acero et al., 2016).

In this analysis, three frequently used spectral models were employed: powerLaw (PL), LogParabola (LOG), and PLSuperExpCutoff (PLEC)⁷⁷7https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html. The binned likelihood analysis method was used to globally fit the photons of IRAS 13052-5711 in the energy range of 0.1-500 GeV. We found that the $\rm TS_{curve}$⁸⁸8$\rm TS_{curve}$ is defined as 2(log $L$(curved spectrum) - log $L$(powerlaw)), and $\rm TS_{curve}>$16 suggests a spectrum has a significant curve referring to Nolan et al. (2012). values of LOG and PLEC were approximately equal to 2. The fitting results showed no significant curvature variability in the spectrum of the target source. Therefore, the PL was selected as the best spectral model for all subsequent analyses. The best-fit spectral index of the PL was 2.09$\pm$0.13, flux of photons was (2.98$\pm$1.04) $\times\ 10^{-9}$ cm^-2 s^-1, and TS value was 42.89 with a significance level of 6.55$\sigma$. The spectral energy distribution (SED) of IRAS 13052-5711 was then divided into six equal logarithmic energy bins, and each energy bin was fitted independently by using the binned likelihood analysis method. For an energy bin of TS$<$4, the upper limit was set at a 95% confidence level, and the relevant results are shown in Figure LABEL:Fig2.

To test whether a significant variability exists in the photon flux of IRAS 13052-5711, we analyzed the LC of IRAS 13052-5711 for 14.4 years. To avoid pollution caused by a large point spread function in the low-energy band, we deducted the data below 200 MeV and used Fermitools to generate the LC between 200 MeV and 500 GeV, as shown in Figure 3. For the LC with 20 time bins, the variability index, $\rm TS_{\rm var}\geq$ 36.19, of the target source is determined to be a variable source with a 99% confidence level (Nolan et al., 2012). We calculated $\rm TS_{\rm var}$ to be 32.06, with a variability significance level of 2.16$\sigma$. This result indicated that IRAS 13052-5711 has no significant flux variability. To further verify this result, a constant-flux model was selected to fit the data points with TS$>$4 in the LC. Finally, we found that the reduced $\chi^{2}$ (=$\chi^{2}/\rm N_{dof}$) of a constant-flux model was approximately 0. This result further demonstrated that IRAS 13052-5711 is a stable source without significant variability.

Thompson et al. (2007) predicted a linear relationship between the total infrared luminosity and the ${\gamma}$-ray luminosity for dense starbursts. Ackermann et al. (2012) established an empirical relation for local group galaxies and SFGs between the total infrared (8-1000 $\mu$m) and the ${\gamma}$-ray luminosity with the detection of high-energy ${\gamma}$-ray emission of star-forming systems. The scalar relationship between $L_{\rm 0.1-500\ GeV}$ and $L_{\rm 0.8-1000\ \mu m}$ can be expressed as follows:

Using the BinnedAnalysis function from the Fermitools⁹⁹9https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/python-usage-notes.html, we calculated the luminosity of IRAS 13052-5711 to be (3.28$\pm$0.67)$\rm\times 10^{42}\ erg\ s^{-1}$ ($d$=93.8 Mpc, from Crook et al. (2007)) in the 0.1-500 GeV energy band. Furthermore, the $L\rm_{\rm 0.1-500\ GeV}$ of other known 11 SFGs and SBGs were recalculated using the same period and method as those for IRAS 13052-5711, and the results are shown in Table 1. We compared the $L\rm_{\rm 0.1-500\ GeV}$ values of each SFG and SBG with those of Guo et al. (2019) and Ajello et al. (2020) and found no statistical difference.

Subsequently, formula (1) and the curve-fit() function was used in the scipy software package (Virtanen et al., 2020) to fit the $L\rm_{\rm 0.1-500\ GeV}$ data of 11 known SFGs and SBGs to obtain $\alpha=1.41\pm 0.07$ and $\beta=39.17\pm 0.08$. The related fitted results are shown in the left panel of Figure 4. The luminosity of IRAS 13052-5711 is marked in the left panel of Figure 4; it deviates from the correlation in Equation (1). To visually verify whether IRAS 13052-5711 exceeds the limit of the proton calorimeter, we performed a correlation analysis between $L_{\rm 0.1-500\ GeV}$/$L_{\rm 8-1000\ \mu m}$ and $L_{\rm 8-1000\ \mu m}$, similar to Ackermann et al. (2012) and Guo et al. (2019), then we found it exceeds this limit, as shown in the right panel of Figure 4.

In addition, the correlation between $L_{\rm 1.4\ GHz}$ and $L_{\gamma}$ is important (Yun et al., 2001). However, investigating the current literature and the NED database looking for 1.4 GHz luminosity of IRAS 13052-5711, we haven’t found any measurement. Therefore, we did not analyze the correlation between $L\rm_{\gamma}$ and $L\rm_{1.4\ GHz}$ in this study.

Through cross-identification of spatial locations, we excluded other common $\gamma$-ray candidates, and found that IRAS 13052-5711 and the $\gamma$-ray radiation at this position are interrelated. Based on this, we analyzed the spectral characteristics. We found that its spectrum had no significant curvature variability with $\rm TS_{curve}\approx 2$ and could be fitted by a simple PL spectral function. In addition, its PL spectral index was approximately 2.1, which is close to 2.2 of the SBGs reported thus far. Moreover, we explored the variability characteristic of the LC over 14.4 years. We found that the LC has no significant variability, with a variability significance level of 2.16$\sigma$, and can be well fitted by employing a constant-flux model with a reduced $\chi^{2}\sim$0. This stable characteristic of the photon flux is consistent with those of the SBGs found to date. Additionally, the SFR is regarded as an important evaluation indicator for identifying SBG, which is calculated using the formula (2) from Kennicutt (1998). The reported relation is the following:

where $\epsilon$ is a function dependent on the initial mass and its value is taken as 0.79 (Ackermann et al., 2012); the solar bolometric luminosity $L_{\bigodot}=3.83\times 10^{33}\ \rm erg\ s^{-1}$ (Sanders et al., 2003). We obtained SFR $\approx$ 29.38 $\rm M_{\odot}$ $\rm yr^{-1}$ as per Equation (2); this value is within the range of the SFR ($>$2 $\rm M_{\odot}\ yr^{-1}$) calculated by $L_{\rm 8-1000\ \mu m}$ for the currently discovered SBG candidates (Ackermann et al., 2012). Based on the foregoing discussion, we believe that IRAS 13052-5711 is likely an SBG.

We analyzed whether IRAS 13052-5711 conformed to the well-known correlation between the $\gamma$-ray luminosity and total infrared luminosity. By using 14.4 years of Fermi-LAT data, we reanalyzed the luminosity correlation from 11 SBGs and SFGs and found that values of 14.4 years of GeV luminosity and the best-fit result of $\alpha=1.41\pm 0.07$ and $\beta=39.17\pm 0.08$ are approximately consistent with those obtained by Guo et al. (2019) and Ajello et al. (2020). Subsequently, IRAS 13052-5711 was found to be beyond the relationship between $L_{\rm 0.1-500\ GeV}$ and $L_{\rm 8-1000\ \mu m}$ and to lie above the proton calorimeter limit. For the current luminosity analysis results, we considered the following two scenarios:

Scenario 1: The $\gamma$-ray radiation is produced by the AGN hosted inside the SBG. When a galaxy harbors an active nucleus, a supermassive black hole powers $\gamma$-ray emission. The radiation provided by the AGN and the inner star-forming region exceeds the limit of the proton calorimeter (Wang & Fields, 2018). Currently, exploring the time-varying characteristic of LC as a common means can help us identify whether there is a contribution from the AGN. This method has been used in previous studies, such as for Circinus (Guo et al., 2019), NGC 3424, UGC 11041 (Peng et al., 2019), and Arp 299 (Xi et al., 2020a). However, by investigating the time-varying characteristic of the 14.4 years of the LC of IRAS 13052-5711, we did not find any significant variability. Hence, our analysis does not support the AGN scenario.

Scenario 2: IRAS 13052-5711 might be a proton calorimeter, even if the $\gamma$-ray luminosity is above the Kennicutt relation as shown in Figure 4. This may occur, in fact, under the assumption that SNe inside IRAS 13052-5711 have a high particle acceleration efficiency $\eta$ than what is commonly assumed for Milky Way SNR. According to Equation (3) derived from Guo et al. (2019), a larger value of $\eta$ increases the ratio $L_{\rm 0.1-500\ GeV}/L_{\rm 8-1000\ \mu m}$ causing a deviation from the empirical relation given in formula (1).

Previous studies have shown that Arp 220 is a probably hadronic calorimeter and that the internal proton and proton interaction (PP) dominates its internal $\gamma$-ray radiation (Wang & Fields, 2018; Peng et al., 2016; Torres, 2004; Lacki & Thompson, 2013; Yoast-Hull et al., 2015). Wang & Fields (2018) considered Arp 220 to harbor a large number of SNe that could systematically and more efficiently accelerate the internal CRs than in the Milky Way under the assumption of a proton calorimeter. SNe inside the Milky Way are usually assumed to have an $\eta$ of 3-10% (Strong et al., 2010). Wang & Fields (2018) assumed that the SNe in the proton calorimeter have a higher $\eta\sim$30% than the one in the Milky Way, and the proton injection spectral index is between 2.0 and 2.4. However, Wang & Fields (2018) found that Arp 220 was still above the inferred proton calorimeter limit, suggesting that the $\eta$ of protons in CRs per SNR inside Arp 220 could be higher than those of previously discovered SBGs, excluding the contribution of the AGN. As the LC of IRAS 13052-5711 does not exhibit significant variability and the $\gamma$-ray luminosity and distance of IRAS 13052-5711 are close to those of Arp 220 (Peng et al., 2016; Xiang et al., 2021a), the scenario in which IRAS 13052-5711 is a proton calorimeter cannot be ruled out.

Based on the assumption of proton calorimeter, we selected a hadronic model to simulate the spectrum of IRAS 13052-5711. Two types of particle energy distribution models were selected, and their expressions are as follows:

Model 1: a power law function (PL):

Model 2: a power law with an exponential cutoff function (ECPL):

where $N_{\rm 0}$ is the amplitude, $E$ is the particle energy, $\alpha$ is the spectral index, and $E_{\rm cutoff}$ is the cut-off energy, and $E_{\rm 0}$ is adopted to be 1 TeV.

In the case of PL, Abdalla et al. (2018) found that NGC 253 exhibited a hard spectral signature with a spectral index$\sim$2.1 of Fermi-LAT. They used the PL particle energy spectrum to fit the observed data. As IRAS 13052-5711 also presents a hard spectrum with the same spectral index of $\sim$2.1, we adopted the PL model for this analysis.

We also decided to test the ECPL looking for the possible presence of a cutoff. The presence of a cutoff is indeed expected in the scenario where particles are accelerated at SNRs inside the starburst nucleus, however, if the process is similar to the Milky Way, we do expect a proton cutoff at energies $>$ 100 TeV, hence for the $\gamma$-ray spectrum the cutoff energy should be located at E$>\sim$10 TeV. Actually, given the larger number of SNR in SBGs, it is probable that the maximum energy of accelerated protons is even larger than the one in the Milky Way (see Peretti et al., 2019). Hence the possible presence of a cutoff in the Fermi-LAT energy range would imply a physical mechanism rather different from the standard SNR acceleration.

A spectral analysis tool named Naima (Zabalza, 2015, and references therein), was used to perform this analysis. Here, the PP differential cross-section of PYTHIA 8 (Kafexhiu et al., 2014) was selected to fit the spectrum of IRAS 13052-5711. For the absorption of extragalactic background light (EBL), here we considered the EBL model from Domínguez et al. (2011) and used redshift z=0.02 (Strauss et al., 1992) of IRAS 13052-5711 in all the following analyses. Considering that the distance of IRAS 13052-5711, the spectral index, GeV photon flux, and gamma-ray luminosity are close to those of Arp 220 (Peng et al., 2016; Xiang et al., 2021a), we assumed that they have a close interstellar medium density $n_{\rm ISM}$=3500 $\rm cm^{-3}$, referring to Peretti et al. (2019). The PL and ECPL were selected to fit the Fermi-LAT observed data points in Figure LABEL:Fig2. The Bayesian information criterion (BIC) was used to determine the goodness of fit between the two models (Schwarz, 1978; Ambrogi et al., 2019). The related BIC formula is $\rm log(n)k-2log(L)$, where n represents the number of the observed data points, k represents the number of model parameters, and L is the maximum likelihood value. The model with the lowest BIC value is preferred. The relevant parameters obtained are listed in Table 2.

By conducting numerical simulation, we found that the values of reduced $\chi^{2}$ of PL and ECPL models are approximately 0. By comparing the BIC values obtained from the fitting of two models, we found no significant difference, as shown in Table 2. These results show that both PL and ECPL can both well describe the $\gamma$-ray emission from IRAS 13052-5711, as shown in Figure 5. However, the cutoff value of $\sim$ 11 TeV should only be considered as a lower limit for a possible cutoff, given that the data does not show any clear departure from a simple power-law.

In addition, we find that the proton slope $\alpha$ for the PL model is approximately 2.22, which is in close agreement with that of NGC 253 with $\alpha\sim$2.2 from the work of Abdalla et al. (2018). Such a hard particle spectrum is consistent with the proton calorimetric scenario in which the equilibrium spectrum inside the starburst nucleus is close to the spectrum injected by SNRs. In fact, the spectral slope of CRs injected by SNRs in the Milky Way is between 2.2 and 2.3 in fully agreement with our result of $\sim$2.2. The slope inferred for IRAS 13052-5711 is also similar to other SBGs like NGC 253 and M 82, while the other famous example, Arp 220, shows a steeper spectrum (2.45 - see Peretti et al. (2019)). The latter case, Arp 220, has parameters very similar to IRAS 13052-5711, hence we compare the proton energy content of the two SBGs: from Table 4 of Peretti et al. (2019) the energy density of the protons is $U_{\rm p}$ = 1324 eV cm^-3 and the radius is $R$ = 250 pc; hence, $W_{\rm p}=\frac{4}{3}\pi R^{3}U_{\rm p}\approx$4e54 erg. Such a value is only twice smaller that the value we estimated for IRAS 13052-571 (see Table 2), which reinforces the idea that IRAS 13052-5711 could be a proton calorimeter. Following Wang & Fields (2018), we have tested the scenario with an acceleration efficiency as large as $\eta$= 30%, finding that such a value is not enough to explain the observed value of $L_{\rm 0.1-500\ GeV}$/$L_{\rm 8-1000\ \mu m}$ in IRAS 13052-5711. This is in tension with results from the theory of particle acceleration at SNRs, where acceleration efficiency larger than 30% is difficult to achieve. Possible explanations of such a discrepancy may be: i) the star formation rate is larger than what estimated in this work or ii) there is an additional source of particle acceleration beyond SNRs. One possibility for ii) is the acceleration at the termination shock of stellar winds [see e.g. Bykov et al. (2020); Morlino et al. (2021)]. Further data below 400 MeV, that may be collected in the future, will help to confirm the hadronic origin of the $\gamma$-ray emission.

Finally, we considered the inverse Compton scattering (IC) and bremsstrahlung radiation (Brem) from relative electronics to interpret the GeV SED of IRAS 13052-5711. For Brem, we took ionized gas density $n_{\rm ion}=87.5\ \rm cm^{-3}$, referring to Peretti et al. (2019). For IC, the cosmic microwave background radiation field was considered. Through simulation with PL and ECPL, we found that IC and Brem can effectively explain the GeV SED of IRAS 13052-5711 with ${\chi}^{2}/N_{dof}\sim 0$. The energy contents of leptons are $W_{\rm e,Brem}\sim 4e55$ erg and $W_{\rm e,IC}\sim 4e60$ erg.

Now, let us assume that electrons are produced by acceleration in SNRs. The maximum energy that can be accumulated in accelerated electrons is given by

Where $\xi_{\rm e}$ is conversion efficiency of SN energy into electrons ($<$1 but usually $\sim$1e-3 for SNR in the Milky Way), $E_{\rm SN}$ = 1e51 erg is the average SN explosion energy, $\Gamma_{\rm SN}$ is the SN explosion rate that can be easily estimated from the SFR provided in the paper just below Eq.(2) (SFR=29.38 $M_{\odot}$/yr)). In fact assuming a Salpeter initial mass function (Salpeter, 1955), the rate of SN can be estimated as the rate of production of stars with mass $>$ $M_{\rm c}$ = 8$M_{\odot}$, hence:

Where $M_{\rm min}$ and $M_{\rm max}$ are the minimum and maximum stellar mass, respectively assumed to be 0.08 and 150$M_{\odot}$. Finally $\tau_{\rm loss}$ is the loss time for electrons inside the nucleus. This should be properly calculated in a self-consistent model, however, given the already mentioned similarity with Arp 220, we will use the maximum possible value estimated in Peretti et al. (2019) (see their figure 1 top panel), namely $\tau_{\rm loss}\sim 10^{5}$ yr. Using those numbers we get $W_{\rm e}=\xi_{\rm e}*$ 2e55 erg. Such a value is a very optimistic upper limit, using more realistic value of $\tau_{\rm loss}\sim 10^{3}$ yr and $\xi_{\rm e}$= 1e-2, we got $W_{\rm e}\sim$ 2e51 erg. Hence leptonic emission is ruled out in our opinion.