A $\gamma$-ray emitting NLS1 galaxy SDSS J095909.51+460014.3

The first 13.8-yr (i.e. 2008 August 4 to 2022 June 4) Fermi-LAT Pass 8 data (evclass = 128 and evtype = 3) are collected, with energy range from 100 MeV to 500 GeV. The Fermitools software (version 2.0.8), along with Fermitools-data (version 0.18), are adopted. In order to filter the photon data, the zenith angle cut (i.e. $<90^{\circ}$) as well as the quality-filter cuts (i.e. DATA_QUAL>0 && LAT_CONFIG==1) are applied in the gtselect and gtmktime tasks. gtlike task with Unbinned likelihood approach is used to extract $\gamma$-ray flux and spectrum. Test statistic (TS = 2$\Delta\log\zeta$, Mattox et al. 1996), where $\zeta$ represents maximum likelihood values of different models with and without the target source, is used to qualify the significance of $\gamma$-ray detection. During the likelihood analysis, diffuse $\gamma$-ray emission templates (i.e. gll_iem_v07.fits and iso_P8R3_SOURCE_V3_v1.txt) together with 4FGL-DR3 sources (Abdollahi et al., 2022) within 15$\degr$ of 4FGL 0959.6+4606 are considered. Parameters of the inner region background sources, a 10$\degr$ radius region of interest (ROI), as well as normalizations of the two diffuse templates are left free, while others are fixed in the 4FGL-DR3 values. We check whether there are new $\gamma$-ray sources based on the subsequently generated residual TS maps. If so, the likelihood fitting is then re-performed adopting the updated background model. When extracting $\gamma$-ray light curve, weak background sources (i.e. TS $<$ 10) are removed from the analysis model.

There is one visit with an exposure time of 1.9 ks from the Neil Gehrels Swift Observatory (Gehrels et al., 2004) on candidate B at 8th May 2017. The data are analyzed by the FTOOLS software (version 6.28). For the XRT photon-counting mode data, firstly, event cleaning xrtpipeline procedure with standard quality cuts is performed. Only 12 net photons of the target have been detected. Nevertheless, signal to noise ratio (i.e. S/N) of detecting the X-ray source is given as 3.4 by the ximage task. No significant excess beyond the background at the position of candidate A is found. The spectrum of source is extracted from a circular region with a radius of 12 pixels, that of background from a circular blank region with a larger radius (50 pixels). The xrtmkarf task is adopted to create the ancillary response files with the most recent calibration database. Considering the limited statistic, the data is not binned and we freeze the absorption column density to the Galactic value (i.e. $\rm 3.6\times 10^{20}$ $\rm cm^{-2}$, HI4PI Collaboration et al. 2016). Meanwhile, the photon index is set as a routine value (i.e. $\Gamma_{x}=1.5$, , where $\Gamma_{x}$ is photon index of the power-law function). An unabsorbed 0.3-10.0 keV flux of $\rm 3.6^{+2.2}_{-1.6}\times 10^{-13}$ erg $\rm cm^{-2}$ $\rm s^{-1}$ ($\mathcal{C}$-Statistic/d.o.f, 9.6/11; Cash 1979) is yielded. In addition, there is a snapshot in UVOT uw2 band, from which the magnitude is extracted by the aperture photometry (i.e. the uvotsource task). During the extraction, a 5$\arcsec$ circular aperture is selected for the target while a larger source-free region is for the background.

We collect the multi-photometric monitoring data in W1 and W2 bands (centered at 3.4, 4.6 $\mu$m in the observational frame) from Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) and the Near-Earth Object WISE Reactivation mission (NEOWISE-R; Mainzer et al. 2014; WISE Team 2020). We remove the bad data points with poor image quality (”qi_fact”$<$1), a small separation to South Atlantic Anomaly (”SAA”$<$5) and in the moon mask areas (”moon mask”=1) (Sheng et al., 2017, 2020). We bin the data in each epoch (nearly half year) using weighted mean value to probe the long-term variability of the target, the corresponding uncertainty is calculated by the propagation of measurement errors.

We further cross-match the archive of the Very Large Array Sky Survey (VLASS; Lacy et al., 2020) and the Faint Images of the Radio Sky at Twenty-Centimeters (FIRST; Helfand et al., 2015) with optical positions of the two objects, candidate A and candidate B. The VLASS and FIRST frequencies, angular resolutions, radio integrated and peak flux densities, and background noise are listed in Table 1. The VLASS fluxes are obtained by modeling the source with a Gaussian fit on the image plane using the Common Astronomy Software Applications.

Analysis of the entire 13.8-yr Fermi-LAT data confirms that there is a significant (TS = 89) $\gamma$-ray source in this direction. The corresponding powerlaw spectrum index is given as $\Gamma_{\gamma}$ = 2.61 $\pm$ 0.14, consistent with the values listed in 4FGL-DR3 (Abdollahi et al., 2022). The optimized $\gamma$-ray location is at R.A. 149.866$\degr$ and DEC. 46.011$\degr$, with a 95% confidence level (C. L.) error radius of 4.6$\arcmin$. Because of the proximity (roughly 4$\arcmin$), both low-energy sources remain within the $\gamma$-ray localization error radius.

Since the spatial resolution of W1 and W2 bands of WISE is about 6$\arcsec$, the separated two infrared light curves of the counterparts provide crucial information of their temporal properties. As shown in Figure 1, variability amplitude of candidate A is about $\lesssim$ 0.5 mag, however, long-term variability with a brightening of $\sim$ 2.5 mag is seen in candidate B. Meanwhile, it is worth noting that the time epochs corresponding to the high flux states of the two sources are different. For candidate A, departed from the high flux state at the beginning of the WISE observations around MJD 55500, a long flux decline until MJD 57500 is detected. Since then, its infrared fluxes maintains to be quiescent. For the RLNLS1, at the start of the WISE observations (i.e. $<$ MJD 57500), no significant variations are detected. However, a distinct flux increase is then followed, peaking on MJD 58063. Other activities are also detected around MJD 59000.

A half-year bin $\gamma$-ray light curve is extracted to compare with those in infrared bands, see Figure 1. Although the statistic is limited, two bins with TS values larger than 10 stand out. There is no corresponding infrared activity of candidate A matching the $\gamma$-ray brightening. However, the epoch of $\gamma$-ray flux increase coincides with the time of high flux state of infrared emission of candidate B, which is supported by the public yearly bin $\gamma$-ray light curve provided in the 4FGL-DR3 catalog (Abdollahi et al., 2022). In addition, temporal properties of strong nearby $\gamma$-ray background sources are investigated, from which no similar behaviors to that of 4FGL 0959.6+4606 are found. Therefore, the brightening is likely intrinsic rather than artificially caused by the background sources.

A monthly bin $\gamma$-ray light curve, focusing on the data later than MJD 57500, is also extracted to search potential short-term variation as well as determine the time range of the high flux state. Two time bins with relatively large TS values appear, locating at MJD 58000 and 59388, respectively. Interestingly, there are a few time bins with TS values $\geq$ 5 clustered around the former one. Therefore, an individual analysis, selecting a time range of 15-month Fermi-LAT data, has been performed, see Figure 2. The analysis yields a significant (TS = 43) $\gamma$-ray source, of which an optimized $\gamma$-ray location at R.A. 149.875$\degr$ and DEC. 45.938$\degr$, along with a 95% C. L. error radius of 6.6$\arcmin$, are derived. A corresponding TS map is generated, where localization results of data from different time epochs are shown together, see Figure 3. In this case, the angular separation between the $\gamma$-ray location and it of candidate A is 7.8$\arcmin$, while the value for candidate B is 5.3$\arcmin$. Therefore, only the RLNLS1 falls into the $\gamma$-ray localization error radius. All the evidences suggest that candidate B is the low-energy counterpart of the $\gamma$-ray source. The photon flux of the $\gamma$-ray source is given as (1.40 $\pm$ 0.40) $\times$ $10^{-8}$ ph $\rm cm^{-2}$ $\rm s^{-1}$ and the $\Gamma_{\gamma}$ is constrained as 2.74 $\pm$ 0.21. Adopting a redshift of 0.399, a corresponding apparent luminosity of (3.1 $\pm$ 0.7) $\times$ $10^{45}$ erg $\rm s^{-1}$ is obtained. For the other monthly time bin at MJD 59388, further 10-day time bin light curve is extracted. In the time range between MJD 59369 and 59399, a $\gamma$-ray source (TS = 30) with a relatively hard spectrum ($\Gamma_{\gamma}$ = 1.90 $\pm$ 0.26) has been found, indicative of possible variability on timescale of a dozens of days. However, due to the limited statistics, in this case, both low-energy sources are embraced within the $\gamma$-ray localization uncertainty area. Moreover, no simultaneous infrared observations at this epoch are available.

To further investigate the validness of the $\gamma$-ray brightening, a 15-month time bin light curve is also extracted, see Figure 4, in which the 7th bin corresponds to the high flux state mentioned in Figure 2 and 3. TS values are lower than 10 (i.e. $<2\sigma$) for the first six time bins in the light curve. However, TS value of the 7th bin reaches to 43 (i.e. 5.9$\sigma$). After the trial factor correction, the detection significance holds at 5.5$\sigma$, suggesting that the brightening is unlikely due to fluctuation of the background emissions. The significance of the flux enhancement is also looked into. Since the detection sensitivity of Fermi-LAT is strongly dependent on length of the exposure time and the target maintains at the low flux state at the beginning, analysis of the entire 7.5-yr data, corresponding to the first six time bins in the light curve, has been carried out. The analysis yields a marginal detection (i.e. TS =24), with a photon flux of (3.84 $\pm$ 0.64) $\times$ $10^{-9}$ ph $\rm cm^{-2}$ $\rm s^{-1}$. A quantity defined as, $\sigma_{var}=\Delta flux/\sqrt{(fluxerr^{2}_{1}+fluxerr^{2}_{2})}$, is used to qualify the significance of variation between two different flux states. Hence the significance is calculated as 2.5$\sigma$, indicating that probability of the $\gamma$-ray flux enhancement is roughly 99%. On the other hand, amplitude of infrared variability is large ($\Delta mag>$ 2.5 mag), and the corresponding measurement uncertainties are relatively small (typically less than 0.1 mag), compared with that in $\gamma$ rays. More importantly, the two infrared light curves exhibit the similar variability trends, suggesting that infrared brightening is not caused by random fluctuation.

Beside the long-term variability, the single WISE snapshots, typically a dozen of exposures within two days, allow us to investigate intraday variability. Assuming a Gaussian noise mode of WISE measurements (Wright et al., 2010), a $\chi^{2}$-test is adopted, of which the null hypothesis corresponds to a constant flux level (Mao & Yi, 2021). The value of the constant flux is optimized for each epochs to find a minimum reduced $\chi^{2}$ value. If such a value is still too large that the null hypothesis is rejected, the variability is decided to be significant. In addition to the estimated PSF photometry errors, systematic errors that are derived from the dispersion of the detections of the standard stars (i.e. $\sim 0.03$ mag, Jiang et al. 2012) have been also considered. No signs of significant infrared short-term variability (i.e. $<3\sigma$) of candidate A are found. For candidate B, significant short-term variability (i.e. ¿ 5$\sigma$) both in W1 and W2 bands around MJD 57334 and 58063 has been detected. Note that the latter is the time of maximum infrared flux value. Variability timescales in the source frame then are estimated as $\tau_{source}=\Delta t\times{\rm ln}2/{\rm ln}(F_{1}/F_{2})/(1+z)$, which suggest doubling timescales of $\simeq$ 9-hr and 17-hr, respectively.

The optical spectrum of candidate B observed on MJD 54525 is derived from the SDSS archive, and its median S/N is about 8. The analysis is with the help of PyQSOFit software¹¹1https://github.com/legolason/PyQSOFit (Guo et al., 2018). After shifting the spectrum to its rest frame and considering the Galactic extinction (Schlafly & Finkbeiner, 2011), firstly, it is simultaneously fitted with a global AGN power-law continuum and stellar contribution of host galaxy while emission lines excluding the Fe ii multiplets (Boroson & Green, 1992; Vestergaard & Wilkes, 2001) are masked. The galaxy emission is described by summation of independent component analysis templates (Lu et al., 2006). Then a simultaneous fitting of the H$\beta$ emission line region (i.e. 4385–5500Å) is carried out, in which the local continuum, the Fe ii multiplets together with the H$\beta$ and the [O iii] doublet emission lines are considered, see Figure 5. However, strength of the Fe ii multiplets can not be well constrained due to the relatively low S/N. He ii $\lambda$4687Å line also falls into this region but is not distinct from the continuum emission. Monte-Carlo (MC) simulations are performed to estimate the measurement uncertainties. Each data point is treated as a gaussian probability distribution, of which the mean and standard deviation are set as the measured value of the data point and its error. Then 1000 times MC samplings for each data point bring mock spectra. By repeatedly fitting the mock spectra, distributions of the parameters are obtained. The measurement uncertainties are calculated based on from the 68% range (centered on the median) of the distributions (Shen et al., 2011). A single Lorentzian profile gives an acceptable description of the H$\beta$ line, of which the full width at half maximum (FWHM) is obtained as (1424 $\pm$ 180) km/s. A two Gaussian model does not bring significant improvement of the fitting. In addition, flux ratio of total [O iii] to total H$\beta$ is estimated as 0.7. Therefore, candidate B is confirmed as a NLS1, consistent with the results of Rakshit et al. (2017). Using the empirical virial BH mass relation (Shen et al., 2008) and the estimated continuum luminosity $\rm\lambda\rm L_{\lambda}$(5100 Å) = (6.45 $\pm$ 2.76)$\times$ $10^{43}$ erg $\rm s^{-1}$, we get a constraint of $\rm log(M_{BH}/M_{\sun})\gtrsim$ 6.8, in case of a disk-like geometry of the broad line region (BLR). Moreover, the Eddington ratio $L_{bol}/L_{edd}\lesssim 0.6$ can be inferred, using the bolometric luminosity correction $L_{bol}\sim(8.1\pm 0.4)L_{5100\AA}$ (Runnoe et al., 2012).

We measure the spectral index between 1.4 and 3.0 GHz by fitting the spectrum with a power-law using the definition

where S1 is the flux density at 1.4 GHz from FIRST and S2 is the average flux density at 3.0 GHz of the two epochs VLASS observations. Candidate A has a steep spectral slope with $\alpha=-0.72$ and candidate B has a flat spectral slope with $\alpha=0.15$. The two sources show a point source morphology in the VLASS and FIRST surveys. The ratio of peak/integrated flux density of candidate A is smaller than that of candidate B, which indicates that the former one is more compact.

Broadband SED of candidate B has been drawn based on the multiwavelength observations, especially focusing on the simultaneous brightening of the $\gamma$-ray and infrared emissions. The high flux state SED includes the 15-month $\gamma$-ray spectrum between MJD 57651 and 58104, Swift observations at MJD 57881, as well as WISE detections at MJD 58063, see Figure 6. Un-simultaneous data including archival GALEX (Seibert et al., 2012), SDSS (Ahumada et al., 2020) and ALLWISE (WISE Team, 2019) photometric data, as well as the VLASS (Lacy et al., 2020) and FIRST (Helfand et al., 2015) radio data, are also plotted. Considering it shares a similar redshift with PKS 1502+036, data of the latter from Abdo et al. (2009b) are drawn for comparison. A remarkable feature of the SED of candidate B is violent variation in the infrared band, which makes the jet emission is overwhelming even in the UV domain then. Shape of the SED in optical-UV bands of candidate B is similar to that of PKS 1502+036, suggesting that the peak frequencies of their synchrotron bump are likely $\lesssim$ $\rm 10^{14}$ Hz. Therefore, candidate B is classified as a low-synchrotron-peaked source (LSP; Padovani & Giommi 1995; Abdo et al. 2010a). For the high energy bump, both of the sources exhibit soft $\gamma$-ray spectra, though PKS 1502+036 is more luminous than candidate B.

The intraday infrared variability allow us to put a constraint of Doppler factor (i.e. $\delta$) of the jet. It should be large enough to avoid severe attenuation on $\gamma$-rays from soft photons via the $\gamma\gamma$ process. The corresponding opacity (Dondi & Ghisellini, 1995) is given as:

where $\sigma_{\rm T}$ is the Thomson scattering cross-section, $R^{\prime}$ is the absorption length, $x^{\prime}$ = $h\nu^{\prime}/m_{e}c^{2}$ and $x^{\prime}_{\rm t}$ correspond to the dimensionless energy of the $\gamma$ rays and target soft photon in the comoving frame, $n^{\prime}(x^{\prime}_{\rm t})$ is the differential number density per energy of the latter. The target photons from the jet itself could be responsible for the absorption. Hence radius of the jet dissipation region, constrained as $R_{j}^{\prime}\leq$ c$\tau_{source}\delta$, where $\tau_{source}$ is set as 17-hr, is equal to the absorption length. The most energetic photon associated with the source is about 4 GeV. Adopting the observed luminosity of soft photons at a few keV of 2 $\rm\times 10^{44}$ erg $\rm s^{-1}$, a constraint of $\delta\gtrsim$ 3 is given. Alternatively, lack of information about the external absorption photons at several tens of eV prevents us from setting a reliable constraint.

The high flux state SED is described by the classic single-zone homogeneous leptonic model including the synchrotron and IC processes, in which the synchrotron self-absorption process and the Klein$-$Nishina effect in the IC scattering are considered. A relativistic compact blob, embedded in the magnetic field and the external photon field, is responsible for the jet emission. Relativistic transformations of luminosity and frequency between the observational frame and the jet frame are $\nu L_{\nu}=\delta^{4}\nu^{\prime}L^{\prime}_{\nu^{\prime}}$ and $\nu=\delta\nu^{\prime}/(1+z)$, respectively. The emitting electrons are assumed to follow a broken power-law distribution,

where $p_{1,2}$ are indices of the broken power-law particle distribution, and $\rm\gamma_{br}$, $\rm\gamma_{min}$ and $\rm\gamma_{max}$ correspond to the break, the minimum and maximum energies of the electrons respectively. Since the accretion system of NLS1s is radiatively efficient, EC processes are proved to be crucial to describe the $\gamma$-ray emission of the RLNLS1s (Abdo et al., 2009a, b). Distance (i.e. $D_{diss}$) between the jet dissipation region and the central SMBH is needed to determine the origin of the external soft photon field (Madejski & Sikora, 2016). Assuming a simple conical jet geometry, the distance is in proportion to the radius of the jet blob $R_{j}^{\prime}$ (Tavecchio et al., 2010), which can be constrained based on the WISE short-term variability. Taking a routine $\delta$ value (i.e. 15, Liodakis et al. 2017), the distance is suggested as $D_{diss}\simeq\delta R_{j}^{\prime}\sim 0.1$ pc. Meanwhile, the radius of BLR can be constrained $\sim 0.02$ pc based on the measured $\rm L_{\lambda}(5100\AA)$ (Bentz et al., 2013). Therefore, in EC process, infrared emissions from the dust torus with temperature of 1200 K (Błażejowski et al., 2000) are considered as the soft photon field, of which the energy density in the rest frame is set as 3 $\rm\times 10^{-4}$ erg $\rm cm^{-3}$ (Ghisellini et al., 2012). In perspective of the broad emission line revealed by the SDSS spectrum, contribution of the accretion disk emission (Shakura & Sunyaev, 1973) has been considered. The accretion disk extends from 3 $\rm R_{s}$ to 2000 $\rm R_{s}$, where $\rm R_{s}$ is the Schwarzschild radius. It produces a total luminosity of $L_{d}=\eta\dot{M}c^{2}$, in which $\dot{M}$ is the accretion rate and $\eta$ is the accretion efficiency. The efficiency is set to 0.057 for a low angular momentum black hole (Frank et al., 2002). The accretion disk emission is characterized by a multi-temperature radial profile (Frank et al., 2002), with a local temperature T at a certain radius $R_{disk}$ is given by,

where $\rm\sigma_{SB}$ is the Stefan–Boltzmann constant.

The theoretical description provides a well explanation of the high flux SED, shown in Figure 6. Corresponding input parameters are listed in Table 2. Lack of sub-mm/far-infrared data as well as a reliable X-ray spectrum obstacle to a further constraint of the modeling. Meanwhile, the disk model can not be well constrained from the high flux SED that the contribution of jet is overwhelming in the UV domain then. Nevertheless, the disk model manages to provide an acceptable description of the archival GALEX data that the jet emission was in the quiescent state. Similar with the known $\gamma$-NLS1s (e.g., Abdo et al., 2009b), $\gamma$-ray emission of candidate B is mainly from the EC processes. Furthermore, the input parameters here are consistent with that used in those sources (e.g., D’Ammando et al., 2012, 2013; Yang et al., 2018; Rakshit et al., 2021; Gokus et al., 2021), such as B $\sim$ 1 Gauss, $\delta\sim 15$ and $\gamma_{br}\sim$ 500. Interestingly, the $\gamma_{br}$ values of $\gamma$-NLS1s, alike of FSRQs, is generally lower than those of BL Lacs (e.g., Chen, 2018). In fact, majority of $\gamma$-NLS1s, are LSPs (Abdollahi et al., 2022). In addition to the leptonic scenarios, hadronic scenarios could also play an important role in $\gamma$-NLS1s, and hence they are possible neutrino emitters (Ahlers & Halzen, 2015). Recently, a cospatial incoming IceCube neutrino event is temporally coincident with a minor $\gamma$-ray flare of 1H 0323+342 (Franckowiak et al., 2020). Hadronic processes have been taken into account to describe the SED of PKS 1502+036 (Wang et al., 2023) We have looked up into the IceCube alert data archive ³³3https://icecube.wisc.edu/science/real-time-alerts/, unfortunately, no known neutrino event in the direction of candidate B is found. Moreover, no evidence of hadronic processes in the electromagnetic data, for example, an orphan $\gamma$-ray flare (e.g., Böttcher et al., 2013), has been found either. Future multi-messenger studies are needed to probe its potential as a neutrino source.