Dissecting the γ -Ray Emissions of the Nearby Galaxies NGC 1068 and NGC 253

Intrigued by recent high-energy study results for nearby galaxies with γ -ray emission and in particular NGC 1068 that has been detected as a neutrino-emitting source by the IceCube Neutrino Observatory, we conduct a detailed analysis of the γ -ray data of the galaxies NGC 1068 and NGC 253, obtained with the Large Area Telescope on board the Fermi γ -ray Space Telescope. By checking their possible spectral features and then constructing light curves in the corresponding energy ranges, we identify spectral-change activity from NGC 1068 in the 2 GeV energy range and long-term, statistically signi ﬁ cant changes for NGC 253 in the 5 GeV energy range. In the former, the emission appears harder in two half-year periods than in the otherwise “ quiescent ” state. In the latter, an ∼ two-fold decrease in the detection signi ﬁ cance after MJD = 57023 is clearly revealed by the test-statistic maps we obtain. Considering the previous studies carried out and the various models proposed for the γ -ray emissions of the two sources, we discuss the implications of our ﬁ ndings. We suspect that a jet ( or out ﬂ ow ) in NGC 1068 might contribute to the γ -ray emission. The nature of the long-term statistically signi ﬁ cant changes for NGC 253 is not clear, but since the part of the GeV emission may be connected to the very-high-energy ( VHE ) emission from the center of the galaxy, it could be further probed with VHE observations.


INTRODUCTION
The launch of the Fermi Gamma-ray Space Telescope (Fermi) and the observations carried out with the large Area Telescope (LAT) onboard it have confirmed the theoretical expectations (e.g., Voelk et al. 1989;Paglione et al. 1996;Domingo-Santamaría & Torres 2005) that star-forming galaxies can emit significant γ-rays (Abdo et al. 2010a).Due to high density of cosmic rays, emanated from copious supernova remnants that are related to the strong star formation (e.g., Ackermann et al. 2012 and references therein), such γray emission is produced through the proton-proton collisions and/or leptonic processes (i.e., bremsstrahlung or inverse Compton scattering; e.g., Dermer 1986).Thus far, more than 10 star-forming galaxies, within the local group or nearby, have been detected at γ-rays (e.g., Ajello et al. 2020;Xi et al. 2020;Xing & Wang 2023).These galaxies lie along a correlation line of the γ-ray versus the infrared (or radio) luminosities (Abdo et al. 2010b;Ackermann et al. 2012;Ajello et al. 2020;Xi et al. 2020), which is considered to indicate the high cosmic-ray densities related to their star-formation property.
However along with this understanding, there are complications.Some of the galaxies contain an active nuclear or other related components (e.g., the jet or outflow), which can contribute γ-ray emissions.The extreme cases are NGC 3424 and UGC 11041, which made themselves detectable with the Fermi LAT by having flare-like events (Peng et al. 2019).Another interesting case is NGC 1068, as it has been shown that its γ-ray emission is more intense than that expected from modeling of the starburst and related cosmic-ray density (Lenain et al. 2010;Yoast-Hull et al. 2014;Eichmann & Becker Tjus 2016;McDaniel et al. 2023).Alternative sources that possibly power the γ-ray emission have been proposed, such as the jet (Lenain et al. 2010) or the outflow (Lamastra et al. 2016) in it.Recently, very-high-energy (VHE) TeV neutrino emission has been detected from this galaxy by the IceCube Neutrino Observatory (Ice-Cube; IceCube Collaboration et al. 2022), making it the second known extra-galactic neutrino-emitting source.By taking the neutrino emission into consideration, Eichmann et al. (2022) have considered a combination of two zones, the starburst region plus the nuclear corona, and Inoue et al. (2022) have suggested the disk winds as the power source.
To fully understand the physical properties of NGC 1068 and other galaxies of the same group (e.g., Ackermann et al. 2012), we have conducted detailed studies of their γ-ray emissions, in which we focused on the variability analysis.In this paper, we report our analysis results for NGC 1068 and another nearby starburst galaxy NGC 253.Spectral and/or flux variations were found, which add more features for helping our understanding of possible physical processes occurring in them.Below we first briefly summarize the properties of the two galaxies related to this work in Section 1.1.The analysis and results are presented in Section 2. In Section 3, we discuss the results.

Properties of NGC 1068 and NGC 253
At a distance of 10.1±1.8Mpc (Tully et al. 2008), NGC 1068 is a nearby galaxy that has been extensively studied at multi-wavelengths (e.g., Pasetto et al. 2019;Inoue et al. 2020 and references therein).It contains an active nucleus in the center, viewed by us with an angle of ∼70 • (Lopez-Rodriguez et al. 2018) and thus classified as type 2. Its GeV γ-ray emission was detected in early observations with the Fermi LAT (Lenain et al. 2010;Ackermann et al. 2012), and no >200 GeV very-high-energy (VHE) emission has been detected (Acciari et al. 2019).Given it being an active galactic nucleus (AGN), its GeV γ-ray emission was searched for variability in studies (Lenain et al. 2010;Ackermann et al. 2012;Ajello et al. 2020Ajello et al. , 2023)), but no significant variability was found.
NGC 253 is one of the nearest galaxies to us (at a distance of 3.5±0.2Mpc; Rekola et al. 2005) and has also been extensively studied at multi-wavelengths.It is nearly edge-on with an inclination angle of i ∼ 72 • (e.g., Puche et al. 1991) and appears with an approximate angular size of 0. • 48×0.• 10 in the optical (Paturel et al. 2003).One particular feature of it is its starburst in the nuclear region (e.g., Rieke et al. 1988;Ott et al. 2005).Related to the starburst, VHE and GeV γ-ray emission from it has been detected with the ground-based High Energy Stereoscopic System (HESS; Acero et al. 2009) and Fermi LAT (Abdo et al. 2010a), respectively.Whether this galaxy contains an AGN or not is under investigation (e.g., Weaver et al. 2002,Müller-Sánchez et al. 2010, and Gutiérrez et al. 2020 and references therein).We selected 0.1-500 GeV LAT events (evclass=128 and evtype=3) from the updated Fermi Pass 8 database in a time range of from 2008-08-04 15:43:36 (UTC) to 2023-04-06 00:00:00 (UTC), approximately 14.7 yr.The region of interest (RoI) for each target was set to be 20 • × 20 • centered at the center of the galaxy.For NGC 1068, R.A. = 40.• 6696 and Decl.= −0.• 01329 (equinox J2000.0), and for NGC 253, R.A. = 11.• 8881 and Decl.= −25.• 2888 (equinox J2000.0).We excluded the events with zenith angles greater than 90 • to reduce the contamination from the Earth limb and included those with good time intervals (selected with the expression DATA QUAL > 0 && LAT CONFIG = 1).The instrumental response function P8R3 SOURCE V3 and the software package of Fermitools-2.2.0 were used in the analysis.

Fermi-LAT Data and Source Model
Source models were generated based on the Fermi LAT 12-year source catalog (4FGL-DR3; Abdollahi et al. 2022) by running the script make4FGLxml.py.They each included all sources in 4FGL-DR3 within a radius of 25 • centered at each of the targets.The positions and the spectral parameters of the sources are provided in the catalog.In our analysis, the spectral parameters of the sources within 5 • from a target were set free and the other parameters were fixed at their catalog values.The spectral model gll iem v07.fit was used for the Galactic diffuse emission, and the spectral file iso P8R3 SOURCE V3 v1.txt for the extragalactic diffuse emission.The normalizations of these two diffuse components were always set free in our following analyses.
For each of the two targets, we set baseline spectral models following those given in the catalog.The γray counterpart to NGC 1068 was modeled as a point source with a power law (PL) spectrum, dN/dE = N 0 (E/E 0 ) −Γ , where E 0 was fixed at the catalog value 1018.59MeV, and the γ-ray counterpart to NGC 253 was modeled as a point source with a log-parabola (LP) spectral form, dN/dE , where E b was fixed at the catalog value 1048.36MeV.Setting a PL for NGC 1068, we performed the standard binned likelihood analysis to the data in 0.1-500 GeV.The obtained results, given in Table 1, are consistent with the catalog values.
We then extracted the γ-ray spectrum of NGC 1068 by performing maximum likelihood analysis to the LAT data in 10 evenly divided energy bins in logarithm from 0.1 to 500 GeV.In the extraction, only the spectral normalizations of the sources within 5 • from NGC 1068 were set free, and all other spectral parameters of the sources in the source model were fixed at the best-fit values obtained from the likelihood analysis.For the obtained spectral data points, we kept those with TS≥4 and derived the 95% flux upper limits otherwise.The spectrum and the PL model fit are shown in Figure 1.For the spectrum, the systematic uncertainties given in Section 3.5 of Abdollahi et al. (2020), which are 10% in 100-300 MeV and 5% in the higher energy ranges, were added in quadrature to the statistical ones.We noted that while the spectrum can be relatively well described with the PL model, the spectral data points with energies approximately greater than 2 GeV appear flat, which provided a hint for our following analysis.

Light curve analysis
We extracted the source's 0.1-500 GeV light curve by setting 180-day a time bin and performing the maximum likelihood analysis to each set of the time-bin data.In the extraction, only the normalization parameters of the sources within 5 • from NGC 1068 were set free and the other parameters were fixed at the best-fit values obtained from the above maximum likelihood analysis.The extracted light curve, for which a systematic uncertainty of 2% (see, e.g., Ajello et al. 2020) was added, is shown in the top panels of Figure 2. The light curve shows possible variations, but when we used the variability index TS var (Nolan et al. 2012) to evaluate, no significant variations were found as TS var ≃ 37.5, lower than the threshold value 49.6 (at a 99% confidence level, for 29 degrees of freedom).
Given the spectrum of the source and our suspicion about the spectral flatness above 2 GeV, we also extracted a 2-500 GeV light curve (bottom panels of Figure 2) with the same setup as the above and the same systematic uncertainty added.In this high energy range, there are two data points (especially the corresponding TS ones), at the beginning ∼MJD 54770 and the middle ∼MJD 57650, possibly indicating two flaring events.However by calculating TS var , we noted that the light curve did not show significant variations, as TS var ≃ 29.5.
In any case, we refined the time periods of the two possible events by obtaining a smooth light curve, for which a time bin of 180 day but with a moving step of 30 day was used.The obtained smooth light curve (with a 2% systematic uncertainty added) and TS curve are shown in Figure 3.The two time bins with the highest TS values were thus found to be MJD 54712.7-54892.7 and MJD 57502.7-57682.7.We defined the first and second time bins as P 1 and P 2 , respectively, and the rest of the time periods as the 'quiescent' P q .
2.2.3.Likelihood and Spectral Analysis for Pq, P1, and P2 We performed likelihood analysis to the data in P q , P 1 , and P 2 , while a PL spectral model was still used for NGC 1068.In this analysis, the spectral parameters of the sources within 5 • from the target were set free and all the other parameters (except the normalizations of the two background components) were fixed at their catalog values.The obtained results are given in Table 1.As can be seen, the emissions in P 1 and P 2 (Γ ≃ 2.02 and 1.85 respectively) were harder that that in P q (Γ ≃ 2.43).
The significances for the spectral differences were ≃2.2σ and 3.2σ respectively.We extracted spectra from the data during the three time periods to check the likelihood results.The same setup and procedures as that in Section 2.2.1 were used.We kept the spectral data points with TS ≥ 4 and showed the derived 95% flux upper limits otherwise.The same systematic uncertainties mentioned in Section 2.2.1 were added.The obtained spectra are shown in Figure 1, with the spectra of P 1 and P 2 compared to that of P q respectively.As can be seen, the spectral models from the likelihood analysis can generally describe the spectra, while only a ∼10 GeV spectral data point appears slightly deviating away from the model fit in P 1 .
For the P 1 and P 2 two possible spectral change events, we carefully checked for other possibilities, such as if they could be induced by the proximity to the Sun, the solar flares1 or gamma-ray bursts (von Kienlin et al. 2020) in the source field, or flares of background quasars during the time periods.The first two possibilities were excluded as no such cases occurred in the RoI during the time periods, which was confirmed by non-detection of short flaring events in such as 1-day-binned light curves we tested to construct.For the third one, there are many quasars, given in the SIMBAD database, that are located ∼2 ′ -30 ′ to NGC 1068, for example [VV2000] J024207.2+000038(Veron-Cetty & Veron 2000; see Figure 4).However, there are no reports of this source being a blazar.
In addition, there are 9 sources listed in 4FGL-DR3 within 5 • of NGC 1068, among which the closest one 4FGL J0242.9+0045(Figure 4) is 0. • 77 away.Among them, six were within 4 • and fainter than our target.We checked if any of them could be variable during P 1 and P 2 by performing binned likelihood analysis to the data for each of them.No significant variations were seen, as the spectral parameters of them during either P 1 or P 2 were consistent within uncertainties with those obtained from the likelihood analysis to the whole data.We also determined the positions of NGC 1068 during P 1 , P 2 , and P q in ≥ 2 GeV energy by running gtfindsrc in Fermitools, and they are shown in Figure 4, a TS map calculated from the ≥ 2 GeV data during the P 1 time period.As can seen, the positions are compatible, excluding the possibility that the apparent hardening of NGC 1068 was instead caused by spectral changes or flares of some nearby sources.

Summary of the results
By checking the ≥ 2 GeV emission of NGC 1068, we found two half-year periods, during which the TS values appeared larger than those of the rest.The corresponding spectra were found to become harder from 2 GeV than the spectrum in quiescence, i.e., having more higher energy photons, and in one time period, the deviation of the spectrum from the latter reached 3.2σ.However, there were trials to be considered in our analysis, such as the used time bin and energy bin in light curves and spectra respectively.As we did not particularly test to find the optimal values for them, where they were approximately chosen based on the TS value of the source, it is hard to evaluate a trial number for our analysis.Also an energy range of ≥2 GeV was chosen based on the visual inspection of the source spectrum.If we take a trial number of ≥100 considering all the numbers of time bins and energy bins, the global significance is largely reduced to be ∼0.Thus we conclude that only Using the source model described in Section 2.1, we performed the standard binned likelihood analysis to the data in 0.1-500 GeV.The obtained results, given in Table 2, are consistent with the values given in 4FGL-DR3, while it can be noted that β is quite small, suggesting a PL model could likely provide an equally good fit.
We then extracted the γ-ray spectrum of NGC 253 by performing maximum likelihood analysis to the LAT data in 10 evenly divided energy bins in logarithm from 0.1 to 500 GeV.The same setup and procedure as those in Section 2.2.1 were used.For the obtained spectral data points, we kept those with TS≥4 and uncertainty/flux ≤ 50%, and derived the 95% flux upper limits otherwise.The spectrum (with the systematic uncertainties added) and the LP model fit are shown in Figure 5.We noted that while the spectrum can be relatively well described with the LP model, there is a possible flux drop-off at 5 GeV, which provided a hint for our following analyses.

Light curve analysis
We extracted the source's 0.1-500 GeV light curve using the same setup and procedure described in Section 2.2.2.The extracted light curve is shown in Figure 6, for which the 2% systematic uncertainty was added.The variability index TS var (Nolan et al. 2012) was determined to be TS var ≃ 24.9, lower than the threshold value 49.6 (at a 99% confidence level for 29 degrees of freedom).Thus no variability was found in the whole data.
Given the spectrum of the source and its possible spectral drop-off above 5 GeV, we also extracted a 5-500 GeV light curve (bottom panels of Figure 6) with the same setup as the above and the same systematic uncertainty added.In this high energy range, TS values of the data points seem to have decreased after ∼MJD 57000.However TS var ≃ 30.7 was obtained for the light curve, not indicating significant variations.
In any case, we defined the two time periods based on the light curve, phase 1 (P1) in MJD 54682.7-57022.7 and phase 2 (P2) in MJD 57022.7-60040.0.

Likelihood and spectral analysis for P1 and P2
In order to check possible property differences between the two time periods, we first performed likelihood analysis to their data.For the target, we tested both LP and PL spectral models.The results are summarized in Ta- ble 2. We found that in P1 and P2, the PL was preferred given the same likelihood values obtained from the both models.The obtained Γ values are consistent with each others within uncertainties.Spectra of P1 and P2 were extracted from the data during the two time periods.The setup and procedure were the same as in Section 2.2.1, and the systematic uncertainties were included in the results.With the spectra shown in Figure 5, it can be seen that no significant flux changes could be determined due to the large uncertainties.
We thus calculated the TS maps centered at NGC 253 in energy range of >5 GeV during , where the second one was set such that its length is equal to that of P1 for direct comparison.The TS maps are shown in Figure 7.As can be seen, the TS value is ∼114 in P1 while ∼26 in P2e.We performed the likelihood analysis to the data in 5-500 GeV in P1 and P2 (P2e), and obtained the fluxes of 11.8±2.4× 10 −11 and 5.6±1.5 × 10 −11 photon cm −2 s −1 (4.1±1.5 × 10 −11 photon cm −2 s −1 ).The flux values are consistent with the TS values as the flux in P2 or P2e is approximately half of that in P1.However, the flux difference between P1 and P2 is not significant due to large uncertainties (similarly indicated by the spectral data points in Figure 5).

Additional analyses
For discussing the γ-ray emission of NGC 253, we conducted additional analyses and provide the results in this section.
First, we checked if the γ-ray emission could be extended or not.In this analysis, the whole 0.1-500 GeV data were used.We set a uniform disk model for the emission in the source model, with its radius ranging from 0. • 1 to 1. • 0 at a step of 0. • 1. Performing the likelihood analysis, the likelihood values from different radii were obtained and compared to that from a point source.The analysis and comparison indicated the emission was consistent with being from a point source.
Second, we determined the position of ≥5 GeV emission of NGC 253, to check if the high energy emission was consistent with being from the center of NGC 253.Running gtfindsrc to the data, we obtained R.A.=11.• 90, Decl.=−25.• 28 (equinox J2000.0), with a 2σ uncertainty of 0. • 04.The error circle contains the center's position of NGC 253 (Figure 7).Note that given the LAT's ∼0.• 25 containment angle at 5 GeV 2 , smaller than the major axis of the galaxy, the emission should arise from the central region.To further check this, we calculated a ≥ 10 GeV TS map (Figure 7), at which energy the containment angle is 0. • 2. The TS map shows that the emission appears smaller than the major axis of the galaxy.

Summary of the results
By checking the emission of NGC 253, we found that its ≥5 GeV part had a possible flux drop around MJD 57023, with the flux before the date being ∼2 times that after the date.This change is revealed by the TS maps calculated for the two time periods.However the change was not significant as the flux uncertainties, particularly that of the second part, are large.Thus we only found a detection significance decrease for NGC 253.In addition, the γ-ray emission in 0.1-500 GeV is consistent with being a point source, and the high-spatial resolution TS maps obtained at higher, ≥5 GeV or ≥10 GeV energies indicate that the emission likely arises from the central region of NGC 253.

DISCUSSION
For the purpose of fully studying the γ-ray emissions of the nearby galaxies NGC 1068 and NGC 253, which can potentially help improve our understanding of the physical processes in them, we carried out detailed analysis.In their γ-ray emissions, we found possible spectral or flux variations at ≥2 GeV or ≥5 GeV high-energy ranges.Below we discuss possible implications of the findings respectively for the two galaxies.

NGC 1068
From the beginning when γ-ray emission was detected from NGC 1068, it has been realized that the emis-sion is more intense than that estimated from typical starburst-induced γ-ray radiation models (Lenain et al. 2010;Yoast-Hull et al. 2014;Eichmann & Becker Tjus 2016).Thus alternative possibilities have been proposed.Lenain et al. (2010) used a jet model, as a kpc-scale radio jet is seen as a component of the galaxy (Wilson & Ulvestad 1983), although recently reported multi-epoch Very Long Baseline Array radio observations of NGC 1068 actually question the presence of (small-scale) radio jets (Fischer et al. 2023).Lamastra et al. (2016) considered the AGN-driven outflow (Krips et al. 2011) that induces shocks to produce high-energy particles.Indeed, AGN-driven outflows have been reported to be detected in γ-rays (Ajello et al. 2021).
Recently to explain the observed neutrino emission, high-energy particles are suggested to be produced in the corona around the central supermassive black hole (Murase et al. 2020;Inoue et al. 2020).While these particles emit γ-rays in the same processes that produce the neutrinos, the optical depth to the high-energy photons is high (e.g., Murase 2022) and thus the photons can not escape the central region until further cascading to MeV energies.By including the processes of the corona with those of the starburst, Eichmann et al. (2022) provided a fit to both γ-ray and neutrino spectra, while the γ-ray one was dominantly contributed by the starburst component.Alternatively, Inoue et al. (2022) proposed disk winds as the power source for the neutrino and γ-ray emissions, with the latter arising from the interaction of the outgoing wind with the surrounding torus.We note that Ajello et al. (2023) very recently conducted analysis of the γ-ray emission of NGC 1068, mainly by extending the spectrum to as low as 20 MeV.They argued that the low-energy, ≤500 MeV part of the emission comes from the corona and the other part from the starburst.
Our results also suggest that the γ-ray emission probably contains multiple components.As the typical starburst models generally work, as proved from studies of other starburst galaxies (e.g., Lenain et al. 2010;Yoast-Hull et al. 2014;Eichmann & Becker Tjus 2016), a component from the starburst should be present.Then at 2 GeV energies, our finding of spectral change events suggests there is another component.Given the variations, it is possibly from the jet (it is not clear whether the outflow scenario would produce variable emission), since jets' emissions are well known to be highly variable.In the model of Lenain et al. (2010), the γ-ray emission is due to the external inverse Compton process (EIC), from a blob of plasma that has a bulk Doppler factor δ b and a radius r b .Taking half a year as the variability time scale ∆t, we may estimate the size of the blob, ∼ c∆tδ b ≃ 0.15δ b pc, where c is the speed of the light.This size value is smaller than that (≃6.5 pc) estimated from the spectral-energy-distribution (SED) fitting in Lenain et al. (2010), but we should note that there are large uncertainties on parameters in the modeling and the size value (or the variability time scale) is in a reasonable range for this type of the emission regions.
For the two variation events, we have checked multiwavelength data to find any possible correlated variations, in particular at X-rays.Oh et al. (2018) reported a hard X-ray light curve in 14-195 keV obtained with the Burst Alert Telescope (BAT) onboard the Neil Gehrels Swift Observatory (Swift), and the light curve ends at ∼MJD 56500.The mission Monitor of All-sky X-ray Image (MAXI; Matsuoka et al. 2009) also detects the galaxy and provides a 2-20 keV light curve.However no correlated variations were found from these two light curves.We note that Marinucci et al. (2016) reported a transient flux excess in ≥ 20 keV hard X-rays however in 2014 Aug. (before the second variation event), and it was likely caused by a temporary decrease of the column density of the obscuring material (in the torus) towards the AGN.Future sensitive hard X-ray or MeV observations of NGC 1068 may help find hints for understanding the γ-ray properties, probing if there are related changes, as the emission in the bands is suggested to be related to the γ-ray emission from the galaxy's central region (Murase et al. 2020;Inoue et al. 2020).

NGC 253
As the starburst models can explain the observed GeV and VHE γ-ray emissions from NGC 253 (e.g., Abdo et al. 2010a;Lacki et al. 2011;Abramowski et al. 2012;Eichmann & Becker Tjus 2016;H. E. S. S. Collaboration et al. 2018) and no flux variations would arise from this sort of starburst-induced emission, it is unexpected to find any variations probed in this work.Both the VHE detection (Acero et al. 2009;Abramowski et al. 2012) and our analysis results point the emission region at the center of the galaxy.We thus suspect any possible flux variations, as hinted by the TS value difference, are possibly related to the central region, where its potential AGN has been searched (e.g., Weaver et al. 2002;Müller-Sánchez et al. 2010).From the work reported by Lehmer et al. (2013), we can see there are four relatively bright X-ray sources in the center, with one of them (N-2003) suggested as an AGN candidate.We checked the hard X-ray light curve of NGC 253 reported from the Swift/BAT Hard X-ray Transient Monitor program (Krimm et al. 2013), but no significant variations could be seen in the light curve.We also checked the X-ray data from targeted observations with the Swift X-ray Telescope (XRT), and no variation patterns were seen in the flux measurements of the central X-ray source (which presumably consists of the four sources) before and after ∼MJD 57023 (2015 Jan. 1).There are also many sets of archival data from Chandra X-ray observations of NGC 253, but by checking the data for the resolved X-ray sources, no clear related variations were found.
Finally, we note that the VHE emission may be closely connected with the ≥5 GeV part.In Figure 5

Figure 2 .
Figure 2. 180-day binned light curves and TS curves of NGC 1068 in 0.1-500 GeV (top panels) and in 2-500 GeV (bottom panels).Fluxes with TS≥4 are kept in the light curves.

Figure 3 .
Figure 3. 180-day binned smooth light curve (top) and TS curve (middle).The time periods P1 and P2, during which the TS values reach to ∼40, are marked by green dash-dotted and yellow dashed lines respectively.

Figure 4 .
Figure 4. TS map of a 3 • ×3 • region centered at NGC 1068 in ≥ 2 GeV energy range during the P1 time period.The positional error circles (2σ values, which are 0. • 12, 0. • 3, and 0. • 04) determined in P1, P2, and Pq time periods, respectively, are marked as dashed green, cyan, and black circles.The pink plus marks the position of NGC 1068 and the green cross marks the position of [vv2000] J024207.2+000038, a nearby quasar example.

Figure 6 .
Figure 6.180-day binned light curves and TS curves of NGC 253 in 0.1-500 GeV (top panels) and 5-500 GeV (bottom panels).Fluxes with TS≥4 are kept in the light curves.In the bottom panels, a blue dotted line indicates MJD 57022.7,an approximate time for the flux changes.

Figure 7 .
Figure 7. TS maps of a 3 • × 3 • region centered at NGC 253 in energy ranges of 5-500 GeV during P1 (left) and P2e (middle) and of 10-500 GeV during the whole time period (right).The optical shape of the galaxy is shown as a green ellipse, with the center marked as a cross.The circle in the left and right panels marks the 2σ error circle determined for the 5-500 GeV emission of the galaxy.
, we show the HESS VHE data points of the galaxy (H.E. S. S. Collaboration et al. 2018), which were obtained from the data collected in years 2005 and 2007-2009.The ≥5 GeV spectral data points in P2 appear to be at the same flux level as the VHE data points.If they are closely connected, similar behavior might be seen in the VHE emission.Unfortunately HESS did not conducted any further observations of the galaxy after year 2009.Hopefully with other VHE facilities at work, possible changes of the VHE emission may be detected in the near future, which could provide hints for further understanding of the γ-ray properties of NGC 253.This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.This research has made use of the MAXI data provided by RIKEN, JAXA and the MAXI team.We acknowledge the use of the Fermi Solar Flare Observations facility funded by the Fermi GI program.

Table 1 .
Likelihood analysis results for NGC 1068

Table 2 .
Likelihood analysis results for NGC 253