Sub-GeV Gamma Rays from Nearby Seyfert Galaxies and Implications for Coronal Neutrino Emission

Recent observations of high-energy neutrinos by IceCube and gamma rays by the Fermi Large Area Telescope (LAT) and the MAGIC telescope have suggested that neutrinos are produced in gamma-ray opaque environments in the vicinity of supermassive black holes. In this work, we present 20 MeV–1 TeV spectra of three Seyfert galaxies whose nuclei are predicted to be active in neutrinos, NGC 4151, NGC 4945, and the Circinus galaxy, using 14.4 yr of Fermi LAT data. In particular, we find evidence of sub-GeV excess emission that can be attributed to gamma rays from NGC 4945, as was also seen in NGC 1068. These spectral features are consistent with predictions of the magnetically powered corona model, and we argue that NGC 4945 is among the brightest neutrino active galaxies detectable for KM3Net and Baikal-GVD. On the other hand, in contrast to other reported results, we do not detect gamma rays from NGC 4151, which constrains neutrino emission from the accretion shock model. Future neutrino detectors such as IceCube-Gen2 and MeV gamma-ray telescopes such as AMEGO-X will be crucial for discriminating among the theoretical models.


INTRODUCTION
Since the IceCube Collaboration announced the first evidence for high-energy cosmic neutrinos in 2013 (Aartsen et al. 2013a,b), the detection of many astrophysical TeV-PeV neutrinos has been reported.They are produced by hadronuclear (pp) and photohadronic (pγ) interactions, inevitably producing a similar amount of gamma rays.These gamma rays can be absorbed by the Breit-Wheeler process (γ+γ → e + +e − ) and reprocessed to lower-energy gamma rays, which can contribute to the cosmic gamma-ray background measured by Fermi Large Area Telescope (LAT; Ackermann et al. 2015).If all the neutrino sources were gamma-ray transparent, gamma rays accompanied by the large all-sky neutrino flux at neutrino energies of E ν ∼ 1 − 10 TeV (Aartsen et al. 2015(Aartsen et al. , 2020a) ) would overshoot the observed diffuse isotropic gamma-ray background.This necessitates the existence of hidden cosmic-ray (CR) accelerators that are opaque to GeV-TeV gamma rays, and the vicinity of supermassive black holes (SMBHs) is among the most promising sites (Murase et al. 2016).
The existence of hidden neutrino sources has been further supported by the recent IceCube observation of neutrino signals from the nearby Seyfert galaxy, NGC 1068 (Aartsen et al. 2020b;Abbasi et al. 2022).The GeV-TeV gamma-ray fluxes seen by LAT (Abdollahi et al. 2020;Ajello et al. 2023) and MAGIC upper limits (ULs; Acciari et al. 2019) are much lower than the TeV neutrino flux measured with IceCube.The region has to be compact for the system to be opaque to GeV gamma rays (Murase 2022), as expected in diskcorona models of active galactic nuclei (AGNs; Murase et al. 2020;Kheirandish et al. 2021).This is the main consequence of the GeV-TeV gamma-ray and TeV-PeV neutrino observations.It has been predicted that X-ray bright AGNs are promising targets for neutrino observations (see Murase & Stecker 2023, as a review), while TeV gamma rays are cascaded down to sub-GeV energies, which will be ideal targets for MeV-GeV gammaray telescopes.
In this work, we investigate the multimessenger implications of three nearby Seyfert galaxies, NGC 4151, NGC 4945, and the Circinus galaxy, using Fermi LAT data.We analyze 14.4 yr of Fermi LAT data to obtain the gamma-ray fluxes from these galaxies.Then, we compare the sub-GeV gamma-ray data with CRinduced electromagnetic cascade emission and discuss corresponding neutrino fluxes.We use Q x = Q/10 x in cgs units and assume cosmological parameters with Ω m = 0.3, Ω Λ = 0.7, and h = 0.7.

NEUTRINO AND GAMMA-RAY EMISSION FROM THE HEARTS OF AGNS
Multimessenger data from NGC 1068 suggest that high-energy neutrinos from Seyfert galaxies are likely to originate from compact regions with a size of R ∼ < 30 − 100 R S (Murase 2022), and X-ray emitting coronae or possible shock regions, located at the vicinity of SMBHs, are the most plausible site of the neutrino production (Murase et al. 2020;Eichmann et al. 2022).In this work, we consider neutrinos and gamma rays from such coronal regions with R = 10 − 30 R S .
The spectral energy distribution (SED) of AGNs is modeled empirically (Murase et al. 2020) based on the parameterization as a function of the SMBH mass M BH and the intrinsic X-ray luminosity L X .We use M BH = 1.0 × 10 7 M ⊙ and L X = 2.6 × 10 42 erg s −1 for NGC 4151 (Bentz et al. 2022;Koss et al. 2022), M BH = 1.4 × 10 6 M ⊙ and L X = 3.3 × 10 42 erg s −1 for NGC 4945 (Greenhill et al. 1997;Puccetti et al. 2014), and M BH = 1.7 × 10 6 M ⊙ and L X = 10 42.5 erg s −1 for the Circinus galaxy (Greenhill et al. 2003;Tanimoto et al. 2019), respectively.We checked that the adopted SED model (see Figure 2 of Murase et al. 2020) agrees with the observed SEDs of the AGNs considered in this work.
High-energy neutrinos are produced through meson production by pγ and pp interactions.We calculate neutrino and gamma-ray emission, considering CR-induced cascades, with the AGN module of the AMES (Astrophysical Multimessenger Emission Simulator; e.g., Murase 2018;Zhang & Murase 2023).Details of our steady-state corona model are described in Murase et al. (2020), but we use the detailed cross sections of the photomeson production and Bethe-Heitler pair production processes as in Murase (2022).The model predicts that TeV neutrinos must be accompanied by sub-GeV gamma rays, predicting gamma-ray and neutrino energy fluxes of E γ F Eγ ∼ (0.1−1)E ν F Eν within 1 order of magnitude.This is a consequence of the multimessenger connection, which is insensitive to model parameters.

Magnetically Powered Coronae
Magnetohydrodynamic (MHD) simulations for accretion flows revealed that the magnetorotational instability (MRI) and subsequent dynamos amplify the magnetic field and develop strong turbulence in accretion flows (e.g., Balbus & Hawley 1998;Kimura et al. 2019).It is believed that magnetic dissipation in the accretion flows heat the plasma and form a hot region called a corona.The observed X-ray emission with a power-law spectrum originates from the Comptonization of disk photons by thermal or bulk electrons.
The MHD turbulence scatters CRs, and the scattered particles randomly change their energies, and effectively gain energy from the turbulence.In the coronal region of AGNs, Alfvén and turbulent velocities are mildly relativistic, enabling the turbulence to accelerate CRs up to TeV-PeV energies.In this process, the resulting CR spectrum can be hard with a spectral index of s ∼ < 2 (for dN CR /dε p ∝ ε −s p , where ε p is the CR proton energy in the source frame), and the maximum CR energy typically lies in the range of ∼ 0.1−1 PeV with a reasonable acceleration efficiency (η tur ∼ 10 − 100).The stochastic acceleration by large-scale turbulence has been studied by test-particle simulations with MHD (Kimura et al. 2016(Kimura et al. , 2019;;Sun & Bai 2021) and large-scale particle-incell (PIC) simulations (Zhdankin et al. 2019;Comisso & Sironi 2019, 2022).Alternatively, magnetic reconnections accelerate nonthermal particles efficiently, especially in the relativistic regime (σ mag = B 2 /(4πnmc 2 ) > 1; Hoshino & Lyubarsky 2012;Guo et al. 2020).PIC simulations have demonstrated that the development of MRI turbulence triggers magnetic reconnections, leading to nonthermal particle production (e.g., Hoshino 2015; Kunz et al. 2016;Bacchini et al. 2022).Kheirandish et al. (2021) suggested such a magnetic reconnection model to explain the neutrino data of NGC 1068 (see also Mbarek et al. 2023;Fiorillo et al. 2023).
In this work, we consider the magnetically powered corona model with stochastic acceleration (Murase et al. 2020), where we solve the Fokker-Planck equation to obtain the CR spectrum and calculate the resulting neutrino and gamma-ray spectra taking into account electromagnetic cascades through γγ → e + e − , synchrotron and inverse-Compton emission processes.The coronal plasma is assumed to be geometrically thick and advected to an SMBH with the infall velocity V fall = α GM BH /R, where α = 0.1 is the viscosity parameter.The magnetic field is obtained through the plasma beta set to β = 1, where synchrotron cascade emission is important for magnetized coronae with β ∼ < 1 and the cascaded ∼ 0.1-10 MeV gamma-ray flux is a robust signature (Murase et al. 2020).Then, the remaining principal parameters of the model are for CR inputs regarding the CR normalization and maximum energy.We consider two specific cases.For model A we assume the CR injection efficiency1 and η tur = 70 that are used to explain the IceCube data of NGC 1068 (see Figure 4 of Abbasi et al. 2022) with R = 30 R S .For model B, the CR pressure P CR is set to 8% of the thermal pressure with the virial temperature (P vir ).Because the coronal pressure would be dominated by magnetic fields, this assumption is reasonable.The emission radius is assumed to be R = 10 R S , in which more neutrinos and gamma rays may be produced due to pγ interactions.For the three AGNs studied in this work, model B gives more optimistic predictions for the neutrino and gamma-ray fluxes.

Accretion Shocks
Accretion shocks have been discussed as a CR acceleration site for a long time (e.g., Stecker et al. 1991;Inoue et al. 2020;Anchordoqui et al. 2021).In this model, a spherical accretion flow, aside from the accretion disk flow, is assumed to form a shock with a velocity of V s ≈ V ff = GM BH /R in the vicinity of the SMBH, and s ∼ > 2 is expected based on diffusive shock acceleration theory (e.g., Drury 1983).However, the neutrino data of NGC 1068 suggest that the particle acceleration efficiency, η sho , has to be much larger than the canonical value in the Bohm limit (η sho ∼ 1 − 10) not to violate the constraint on the neutrino break/cutoff (see Figure 5 of Kheirandish et al. 2021).Thus, in this work, we ad hoc fix the proton maximum energy to 100 TeV (implying η sho ≫ 10) and consider the largest CR luminosity L CR allowed by the LAT data.Then, we calculate neutrino and gamma-ray spectra assuming the steady-state injection of CRs for given R and SEDs.
In the accretion shock model, magnetic fields are expected to be weaker than the magnetically powered corona model, and we set B = 30 G throughout this work (Inoue & Khangulyan 2023).For such low magnetizations, the results on gamma-ray spectra are largely insensitive to CR spectra as long as electromagnetic cascades proceed mainly through γγ → e + e − and inverse-Compton emission.Note that we only consider CR-induced cascade gamma rays without any primary nonthermal electron component, unlike Inoue & Khangulyan (2023), who focused on gamma rays from primary electrons only with attenuation (without cascades).The primary electron component would add an additional sub-GeV gamma-ray component, which would result in a lower neutrino flux in order to be consistent with the given sub-GeV gamma-ray data.

DATA ANALYSIS
We study three nearby Seyfert galaxies, NGC 4151, NGC 4945, and the Circinus galaxy.These objects were all reported as LAT sources and are also very bright in the hard X-ray band.In the magnetically powered corona model, the neutrino luminosity is roughly proportional to the X-ray luminosity (Murase et al. 2020;Kheirandish et al. 2021), and thus NGC 1068 and NGC 4151 (NGC 4945 and Circinus) are predicted to be among the brightest sources in the northern (southern) neutrino sky among Seyfert galaxies in the Swift BASS catalog with LAT detection2 .Prior to this work, the results of NGC 1068 (being the brightest in the IceCube sky) are presented in Ajello et al. (2023).
We analyze data collected by the LAT (Atwood et al. 2009(Atwood et al. , 2013;;Bruel et al. 2018) from 2008 August 04 to 2023 January 5 (14.4 yr).We select an energy range of 20 MeV-1 TeV, and bin the data using eight energy bins per decade.The pixel size is 0.08 • .We use a 10 • × 10 • region of interest (ROI) centered on each respective galaxy.The standard data filters are used: DATA QUAL>0 and LAT CONFIG==1.The analysis Note-These energy bins distinguish different components of the joint likelihood analysis; however, the underlying fits still use eight energy bins per decade, as described in the text.Z Max specifies the maximum zenith angle which controls the Earth limb contamination. is performed using Fermipy (v1.2) 3 , which utilizes the underlying Fermitools (v2.2.0).We select photons corresponding to the P8R3 SOURCE V3 instrument response.
In order to optimize the sensitivity of our analysis, we implement a joint likelihood fit with the four point spread function (PSF) event types available in the Pass 8 data set 4 .The data are divided into quartiles corresponding to the quality of the reconstructed direction, from the lowest quality quartile (PSF0) to the best quality quartile (PSF3).Each subselection has its own binned likelihood instance that is combined in a summed likelihood function for the ROI.This is easily implemented in Fermipy by specifying the components section in the configuration file.We include different event types depending on the energy interval.This is motivated by the energy-dependence of the LAT instrument response.Additionally, in order to reduce the contamination from the Earth's limb, we apply an energy-dependent cut on the maximum zenith angle.These selections are similar to those used for the 4FGL-DR3 catalog (Abdollahi et al. 2022), and they are summarized in Table 3. Note, however, that the 4FGL-DR3 did not use photon energies below 50 MeV.Each PSF type also has its own corresponding isotropic spectrum, namely, iso P8R3 SOURCE V3 PSFi v1, for i ranging from 0 to 3. The isotropic spectrum file is defined from 34 MeV to 878 GeV, and so we use a power-law extrapolation for energies outside this range.The Galactic diffuse emission is modeled using the standard component (gll iem v07).It is only defined between 50 MeV and 814 GeV, so we likewise use a power-law extrapolation.The point source emission is modeled using the 4FGL-DR3 catalog (gll psc v28; Abdollahi et al. 2022).In order to account for photon leakage from sources outside of the ROI due to the PSF of the detector, the model includes all 4FGL sources within a 15 • × 15 • region.The energy dispersion correction (edisp bins=-1) is enabled for all sources except the isotropic component.
Before calculating SEDs, we first perform an initial fit for each of the three ROIs, in which we optimize the model, freeing the normalization and index of the Galactic diffuse emission, the normalization of the isotropic emission, and all point sources (normalization and index) within 3 • .We then find new sources using the Fermipy function find sources, which generates test statistic (TS) maps and identifies new sources based on peaks in the TS.The TS maps are generated using a power-law spectral model with a photon index of 2.0.The minimum separation between two point sources is set to 0.5 • , and the minimum TS for including a source in the model is set to 16.Finally, we perform a second fit in a similar way as the first, with the exception that we free all point sources within 5 • having a TS> 25.
Among the three Seyfert galaxies, the LAT observations of NGC 4151 are complicated by the pres- ence of two nearby BL Lacs, 4FGL J1211.6+3901 and 4FGL 1210.3+3928.In order to determine the most likely origin of the corresponding gamma-ray emission, we relocalize the point source closest to NGC 4151.This is done following a similar procedure as the ROI optimization described above.We use the same data selection, with the exception that we use a pixel size of 0.04 • .In the first step of the optimization, we only include point sources from the 4FGL-DR3; i.e., we do not include a point source at the location of NGC 4151.Then, before finding new sources, we remove 4FGL 1210.3+3928from the model, thus allowing the fitting algorithm to find the best-fit location of the associated source.The resulting localization is shown in Figure 1.The pink contour shows the 95% localization uncertainty, based on the peak in the TS map.Red contours show the 95% localization uncertainty from the 4FGL-DR3, and the blue circles show the BL Lacs which are associated with the 4FGL sources.As can be clearly seen, the gamma-ray source near NGC 4151 is most likely associated with 4FGL 1210.3+3928(and hence the BL Lac), consistent with the 4FGL catalog.We also verified that the results remain consistent when using energies E γ > 1 GeV.Additionally, we checked with the preliminary version of the 4FGL-DR4, also finding consistent results.We therefore keep 4FGL 1210.3+3928 in the model when calculating the SED for NGC 4151.
After optimizing the ROIs, SEDs are calculated for each source using the Fermipy SED analysis.This method computes an SED by performing independent fits for the flux normalization in each energy bin.For the calculation, we combine the original energy bins into larger bins, using two bins per decade between 100 MeV and 100 GeV, as well as a single bin between 20 and 100 MeV and a single bin between 100 and 1000 GeV.The normalization in each bin is fit using a power-law spectral model with a fixed index of 2.0.The parameters of the background components are held fixed at their best-fit values from the baseline fit.Further details of the SED analyses for the three Seyfert galaxies are provided in Appendix A. The resulting SEDs of NGC 4151, NGC 4945, and Circinus are shown in Figures 2 and 3.

IMPLICATIONS AND DISCUSSIONS
GeV gamma rays mainly interact with X rays from the corona.The two-photon annihilation optical depth for a photon index of Γ X ≈ 2 is (Murase 2022) 1) where η γγ ∼ 0.1 is a numerical coefficient depending on Γ X , σ T ≈ 6.65 × 10 −25 cm 2 , εγγ−X = m 2 e c 4 /ε X ≃ 0.26 GeV (ε X /1 keV) −1 , ε X is the reference X-ray energy, and n X ≈ LX /(2πR 2 cε X ) is used.Here LX is the differential X-ray luminosity.The numerical results shown in Figures 2 and 3 are consistent with Equation (1), and we may expect ∼ < 0.1 GeV gamma rays.
For NGC 4945 and Circinus, while GeV emission is detected, its origin should be different, e.g., star-forming activities, as seen in NGC 1068 (Ajello et al. 2023).
Note that hadronic emission from the starburst region predicts a low-energy break in the gamma-ray spectrum due to the decay of neutral pions, although there could be leptonic contributions from inverse-Compton and bremsstrahlung emission.However, detailed studies of the GeV emission are beyond the scope of this work.NGC 4151 is known to be a Seyfert 1.5 galaxy, showing the characteristic features of both types 1 (dominated by broad-line components) and 2 (dominated by narrow-line components).NGC 4151 is Compton thin, unlike the other two Seyferts and NGC 1068.The escape of GeV gamma rays is easier than in NGC 1068 due to its lower intrinsic X-ray luminosity, which makes the strong limits in the GeV range important.As shown in the left panels of Figures 2 and 3, we find that the magnetically powered corona and accretion shock models are consistent with the ULs obtained in this work.Nevertheless, the data give interesting constraints on neutrino emission.In the magnetically powered corona model, model B predicts a neutrino flux that can explain the possible neutrino excess emission found in the IceCube data (Glauch et al. 2023;Liu et al. 2023;Neronov et al. 2023).In the case of NGC 4151, the effective maximum CR energy5 is limited by the diffusive escape rather than the Bethe-Heitler cooling process that is relevant for luminous AGNs, including NGC 1068.On the other hand, due to the LAT ULs in the GeV range, the neutrino flux is unlikely to be explained in the accretion shock model even if magnetic fields change.This demonstrates that observations of NGC 4151-like galaxies are useful for discriminating between the magnetically powered corona and accretion shock models.When we include not only 2-10 keV but also 14-195 keV data, the neutrino brightness of NGC 4151 is next to NGC 1068, where the ranking is higher than that in Kheirandish et al. (2021), and NGC 4151 will be an important target for IceCube-Gen2 which can reach a sensitivity of E ν F Eν ∼ 10 −9 GeV cm −2 s −1 (Aartsen et al. 2021), especially in the magnetically powered corona model (see left panel of Figure 2).
For NGC 4945, as shown in the middle panels of Figures 2 and 3, we find that GeV and higher-energy emission has a break at E γ ∼ 1 GeV, which is consistent with a pionic gamma-ray component from π 0 → γγ, as expected in models relying on wind-torus interactions (Inoue et al. 2022) and/or star-forming activities (Aguilar-Ruiz et al. 2021).In addition, the LAT data around ∼ 30 MeV shows a clear excess compared to the pionic gamma-ray component.However, it is important to note that the LAT performance below 50 MeV quickly deteriorates, and the instrument response at these low energies has not been well characterized.We have performed a number of tests to validate the robustness of the observed low-energy spectral feature in NGC 4945 (summarized in Appendix A), and none of them have shown any indication of the signal being spurious.But nevertheless, this feature should still be interpreted with caution due to the limitations of the LAT performance below 50 MeV.That said, this signal is very interesting because a similar excess flux over the pionic gammaray component was reported for NGC 1068 (Ajello et al. 2023).More intriguingly, the sub-GeV gamma-ray signatures of these two most promising neutrino-bright AGNs are consistent with predictions of the magnetically powered corona model, and the 20 − 100 MeV data of NGC 4945 can be explained if the CR pressure is ∼ (5 − 10)% of the thermal pressure with the virial temperature (see Figure 2 middle).Corresponding to the systematic uncertainty with a factor of 2 − 3 in the lowest-energy bin data, theoretical fluxes can be even lower by reducing P CR and/or increasing η tur .For Circinus, as shown in the right panels of Figures 2 and 3, we only have ULs at sub-GeV energies, which are also consistent with both of the corona and accretion shock models.In all the three Seyfert galaxies, the model B fluxes shown in Figure 2 can be regarded as ULs allowed by the LAT data, while they lie below the COMPTEL sensitivity of E γ F Eγ ∼ a few × 10 −7 GeV cm −2 s −1 (De Angelis et al. 2017).
The corona model predicts that NGC 4945 and the Circinus galaxy are Seyfert galaxies that are among the most promising for neutrino detection, and the all-flavor neutrino flux reaches E ν F Eν ∼ 10 −8 − 10 −7 GeV cm −2 s −1 as seen in Figure 2.They are promising targets for Northern Hemisphere neutrino detectors such as KM3Net (Adrian-Martinez et al. 2016), Baikal-GVD (Avrorin et al. 2020), P-ONE (Agostini et al. 2020), andTrident (Ye et al. 2022).For exam-ple, the sensitivity of KM3Net is E ν F Eν ∼ a few × 10 −9 GeV cm −2 s −1 .We also encourage IceCube searches for starting tracks and showers in the southern sky6 .

SUMMARY
We studied gamma-ray emission from three nearby Seyfert galaxies, NGC 4151, NGC 4945, and the Circinus galaxy.The theoretical models predict that the neutrino flux is approximately proportional to the intrinsic X-ray luminosity.NGC 1068 and NGC 4151 are the most promising neutrino sources in the IceCube sky, while NGC 4945 and Circinus are the best in the southern sky.We analyzed the 20 MeV -1 TeV spectra of the three Seyfert galaxies, NGC 4151, NGC 4945 and Circinus, using 14.4 yr of the Fermi LAT data.We derived ULs on gamma-ray emission from NGC 4151, which gives strong constraints on the models of highenergy neutrino emission from this source.The recent hint of a neutrino excess toward NGC 4151, reported by the IceCube Collaboration (Glauch et al. 2023;Liu et al. 2023), appears to be consistent with the magnetically powered corona model.Future observations of NGC 4151 and other Compton-thin AGNs with nextgeneration neutrino detectors such as IceCube-Gen2 are crucial for discriminating the different models.We also found evidence of sub-GeV gamma-ray emission in the direction of NGC 4945.If it originates from NGC 4945, a component other than pionic gamma-ray emission is needed.Intriguingly, this excess is consistent with CRinduced electromagnetic cascades expected in the magnetically powered corona model, although other interpretations such as high-energy emission from jets would also be possible.Cascade signatures expected in the MeV and sub-GeV gamma-ray bands are important for future MeV gamma-ray detectors such as AMEGO-X (Caputo et al. 2022)  NGC 4151 is a nearby Seyfert 1.5 galaxy, located at a distance of d = 15.8Mpc (Yuan et al. 2020).X-ray observations of NGC 4151 indicate that it contains an ultra-fast outflow (UFO) (Tombesi et al. 2010(Tombesi et al. , 2012)).Evidence for high-energy gamma-ray emission from UFOs, including NGC 4151, was first reported in Ajello et al. (2021), based on Fermi LAT observations.In that work, a stacking analysis was performed using a small sub-sample of the UFO population, resulting in a 5.1σ detection.Although NGC 4151 itself was below the LAT detection threshold, it had an individual significance of 4.2σ, which was the highest of the sample.It should be noted that when excluding NGC 4151 from the UFO sample, the population was still detected at the level of 3.5σ, with similar spectral properties.This showed that the signal was not due to NGC 4151 alone.It should also be noted that in the stacking analysis, the source locations are fixed at their optical positions.
Following the detection of the sub-threshold UFO sample, the individual detection of NGC 4151 was recently reported in Peretti et al. (2023).Importantly, NGC 4151 is located at a distance of only 0.18 • from the fourth most significant hot spot in the search for northern neutrino point sources (Abbasi et al. 2022), and the gamma-ray detection could potentially have significant implications for the nature of the IceCube neutrinos.
The LAT observations of NGC 4151 are complicated by the presence of two nearby BL Lacs.In fact, one of these sources (4FGL J1211.6+3901,located 0.43 • away) was carefully analyzed in Peretti et al. (2023), and it was determined that the detection of NGC 4151 was indeed robust against the presence of this nearby gamma-ray bright blazar.The point source model in Peretti et al. (2023) was originally based on the 4FGL-DR2 catalog; however, in the most recent catalog (4FGL-DR3) there is now another source (4FGL 1210.3+3928) that is located even closer, at a separation of 0.08 • , associated with the BL Lac, 1E 1207.9+3945.It therefore seems that the gamma-ray emission thought to be coming from NGC 4151 is in fact dominated by a nearby blazar.
Based on the localization results (presented in Section 3), we keep the nearby source (4FGL 1210.3+3928) in the model.We add another source at the optical center of NGC 4151 and refit the data.Unsurprisingly, NGC 4151 is no longer detected (TS ≈ 0), and so we calculate ULs, shown in the left panel of Figure 2, with corresponding data provided in Table A1.The bins have TS=0, and so we calculate the ULs using a Bayesian approach (as opposed to a frequentist approach).Specifically, we use the calc int method from the IntegralUpperLimit class, available in the Fermi Science Tools.This method calculates an integral UL by integrating the likelihood function up to a point which contains a given fraction of the total probability.We verified that the field is well-modeled, with good datamodel agreement seen in the fractional count residuals and spatial residuals.We note that in a revised version of their work, Peretti et al. (2023) exclude on the basis of SED modeling that 1E 1207.9+3945contributes to the emission of the source they associate to NGC 4151.In our work (and in the upcoming 4FGL-DR4) we clearly find that the position of the gamma-ray source 4FGL 1210.3+3928 is compatible with 1E 1207.9+394.

A.2. NGC 4945
We calculate the SED for NGC 4945, located at the distance d = 3.6 Mpc (Monachesi et al. 2016), using the same fitting approach as was used for NGC 4151 (in this case there is no need to relocalize).Results for this are shown in the middle panel of Figure 2, and the corresponding data are provided in Table A2.The source is detected with a TS=830.Overall, the SED is consistent with the values from the 4FGL-DR3.However, compared to the 4FGL, the two lowest-energy bins now have significant signals (TS>9), as opposed to ULs.For the second energy bin, this is quite reasonable given the additional exposure.Interestingly, extending the lowestenergy bin down to 20 MeV also results in a significant signal.
The low-energy detection of NGC 4945 has important implications for the corona model, and so we perform additional tests to check the robustness of the emission.The 68% containment angle at 20 MeV for PSF3 events is ∼ 13 • , and the 95% containment angle is ∼ 30 • .However, our ROI size is only 10  20 − 100 MeV, and below 100 MeV, the effective area quickly falls off.Specifically, the (total) effective area at 100 MeV is ∼ 4.4× higher compared to 20 MeV.More of the events will therefore be coming towards the upper part of the bin, where the 68% containment angle is ∼ 3.3 • .Moreover, for a Gaussian PSF, the distribution peaks towards the center.On the other hand, the behaviour of the LAT PSF below 50 MeV has not been well characterized.We therefore test the robustness of the emission with respect to the ROI size.
We repeat the analysis for NGC 4945 using two different regions, motivated by the approximate 68% and 95% containment radii.For one test we use a 20  ×60 • region.Note that for the latter test, the data cover the entire 68% containment region, and the model sources cover the entire 95% containment region.Figure A1 shows the resulting SEDs, significance ( √ TS), and fractional count residuals for all three fits.Additionally, we overlay the SED and corresponding TS from the 4FGL-DR3.We find consistent results for all three fit variations.The flux in the first bin is compatible with the 95% ULs from the 4FGL-DR3.The flux for the second energy bin is slightly lower for the 10 • ROI, but it is still consistent with the other fits within uncertainties.The fractional count residuals show that the data are well modeled, with agreement better than 5% for the full energy range.The data are specifically well modeled in the first two energy bins.The spatial residuals are shown in Figure A2 (calculated with the fermipy tool residmap) for both the entire energy range (top row) and the first energy bin (bottom row).For the full energy range, the 10 • ROI is very well modeled, whereas the larger ROIs show some excesses/deficits towards the edges of the field.For the first energy bin, all fit variations show good data-model agreement.
In order to better understand the reason for the discrepancy in the first energy bin with respect to the 4FGL-DR3, we divide it into two smaller bins to match the 4FGL binning, and recalculate the SED.The lower part of the bin spans the energy range 20 − 45 MeV, and the upper part of the bin spans the energy range 45 − 100 MeV.We make the split at 45 MeV (compared to 50 MeV for the lower edge of the 4FGL bin), as required by our initial energy binning.This test is performed for all three ROI variations, and the results are shown in Figure A1 with open markers, corresponding to the respective ROI size.
For our nominal ROI size of 10 • , we find a significance ( √ TS) of 5.8 and 2.6 for the two bins.The significance of the upper part of the bin is now reasonable given that of the 4FGL, and the flux is also compatible with the 4FGL ULs.Additionally, we find a significant signal in the lower portion of the bin.However, the corresponding count residuals in the lower portion of the bin are slightly elevated at the level of ∼ 5%.This can plausibly be attributed to some mismodeling of either the isotropic emission or the Galactic diffuse emission (which at these energies is dominated by inverse-Compton radiation).Note that the isotropic model uses a power-law extrapolation below 34 MeV, and the Galactic diffuse model uses a power-law extrapolation below 50 MeV.For the 20 • ROI, the flux of the two smaller bins is found to be consistent with that of the combined bin, and likewise the count residuals are close to 0. For the 26 • ROI, the upper portion of the bin is also compatible with the combined bin; however, the lower portion of the bin is over-modeled by ∼ 3%, and the SED calculation results in only an UL, albeit at the same level as the observed signal.
As another test, we have made the SED calculation using a spectral index of 3.0 (instead of the commonly used value of 2.0), consistent with our model predictions, and also freeing all point sources within 5 • , having a TS>16.This test was run for the 20 • ROI.We have found that the low-energy bin (E γ < 100 MeV) is still detected with a similar flux having a TS of 19.0.Finally, we have tested the SED calculation when freeing all sources within 7 • having a spectral index > 2.5.In this case, we have found consistent results for the low-energy bin having a TS of 27.5.
Overall, these tests show that the emission between 20 − 100 MeV appears to be robust, although there is clearly a systematic uncertainty in the measured flux of a factor of ∼ 2 − 3.However, this result should still be interpreted with caution, due to the generally poor performance of the LAT below 50 MeV.

A.3. Circinus
We also calculate the SED for Circinus, located at d = 4.2 Mpc (Freeman et al. 1977); the results are shown in the right panel of Figure 2, and the corresponding data are provided in Table A3.The source is detected with a TS=130.Overall, the SED is consistent with the values from the 4FGL-DR3.The fractional count residuals are very good, being close to 0 for all energies.

Figure 1 .
Figure 1.Left: TS map of the NGC 4151 field.Contours show the 95% localization uncertainty.Red contours are from the 4FGL-DR3 catalog, and the pink contour is calculated based on the peak in the TS map.Blue circles show the locations of the BL Lacs that are associated with the 4FGL-DR3 sources (4FGL 1211.6+3901towards the lower left, 0.43 • separation, and 4FGL 1210.3+3928closest to NGC 4151, 0.08 • separation).Right: Same as the left plot but zoomed in on NGC 4151.Overlaid are multiwavelength images of NGC 4151 in both optical (black and white image) and radio (H I, shown with the rainbow color map).

Figure 2 .Figure 3 .
Figure 2. Multimessenger SEDs of NGC 4151, NGC 4945 and Circinus galaxy in the magnetically powered corona model.Black points and ULs are Fermi LAT data in gamma rays (this work).The error bars on the LAT data points are 1σ, and the ULs are plotted for bins with TS<4, and they are shown at the 95% C.L. The sky-blue shaded region represents IceCube data in neutrinos(Glauch et al. 2023;Liu et al. 2023).Sensitivities for e-ASTROGAM (De Angelis et al. 2017) and GRAMS(Aramaki et al. 2020) with an effective exposure time of 1 yr are overlaid, together with the AMEGO-X sensitivity for the 3 yr mission(Caputo et al. 2022).For model A, we use parameters to explain NGC 1068 as inMurase et al. (2020), where the emission radius is set to R = 30RS with ηtur = 70.For model B, we assume R = 10RS and the CR pressure is set to PCR/Pvir = 8%, where ηtur = 100 for NGC 4151 and ηtur = 10 for NGC 4945 and Circinus are used.
e-ASTROGAM(De Angelis et al. 2017), and GRAMS(Aramaki et al. 2020).the support by KAKENHI Nos.22K14028, 21H04487, and 23H04899 and the Tohoku Initiative for Fostering Global Researchers for Interdisciplinary Sciences (TI-FRIS) of MEXT's Strategic Professional Development Program for Young Researchers.This work was supported by the European Research Council, ERC Starting grant MessMapp, S.B. Principal Investigator, under contract No. 949555.The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis.These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden.Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d' Études Spatiales in France.This work performed in part under DOE Contract DE-AC02-76SF00515.Work at NRL is supported by NASA.APPENDIX A. DETAILS OF GAMMA-RAY DATA ANALYSES A.1.NGC 4151

Figure A1 .Figure A2 .
Figure A1.NGC 4945 results for three different ROI sizes, as specified in the legends.Results from the 4FGL-DR3 are also overlaid.The left panel shows the SED, the middle panel shows the corresponding significance ( √ T S), and the right panel shows the count residuals.In addition to the nominal binning, the first energy bin is divided into two smaller bins, shown with open markers corresponding to the ROI size.This is done to make a more direct comparison to the 4FGL-DR3, as discussed in the text.In the middle panel, we indicate the predicted significance based on the additional exposure time for both the noise-limited case ( √ t) and the signal-limited case (t).The error bars on the data points are 1σ, and the ULs are plotted for bins with TS<4, and shown at the 95% C.L.

Table A1 .
NGC 4151 SED Data • × 10 • , and the model includes all sources within 15 • × 15 • .Although the ROI is smaller than the PSF at 20 MeV, it should be sufficient for the analysis, given the following consideration.The lowest-energy bin spans the energy range ROI, and include all sources within a 40 • × 40 • region.For another test, we use a 26 • × 26 • ROI and include all sources within a 60 • × 20 •