The Awakening of a Blazar at Redshift 2.7 Temporally Coincident with the Arrival of Cospatial Neutrino Event IceCube-201221A

We report on multiwavelength studies of a blazar NVSS J171822+423948, which is identified as the low-energy counterpart of 4FGL J1718.5+4237, the unique γ-ray source known to be cospatial with the IceCube neutrino event IC-201221A. After a 12 yr long quiescent period undetected by Fermi-LAT, γ-ray activities with a tenfold flux increase emerge soon (a few tens of days) after the arrival of the neutrino. Associated optical flares in the Zwicky Transient Facility g, r, and i bands are observed together with elevated Wide-field Infrared Survey Explorer infrared fluxes. Synchronized variations suggest that both the γ-ray emission and the neutrino event are connected to the blazar. Furthermore, the optical spectrum reveals emission lines at a redshift z = 2.68 ± 0.01. Thus, it is the first candidate for a neutrino-emitting blazar at the redshift above 2. Discussions of theoretical constraints of neutrino production and comparisons with other candidates are presented.


Introduction
The IceCube Neutrino Observatory has detected high-energy astrophysical neutrinos (e.g., IceCube Collaboration 2013; Aartsen et al. 2015;IceCube Collaboration et al. 2023), ushering in a new era in neutrino astronomy.The neutrino signal is compatible with an isotropic distribution in the sky, suggesting that a substantial fraction of the observed neutrinos are of extragalactic origin.Neutrinos interact weakly with matter and travel undeflected in magnetic fields, providing key insights into extreme cosmic environments that are otherwise opaque (Aartsen et al. 2018).Interactions involving photohadronic (pγ) or hadronuclear (pp) processes are responsible for the generation of neutrinos, while accompanying γ-ray radiations are also anticipated (Halzen & Kheirandish 2022).Besides routine leptonic radiations, hadronic cascades also produce electromagnetic emissions at lower energies.Hence, potential connections between electromagnetic emissions and neutrino events are suggested.Various types of extragalactic γray emitters are proposed as potential neutrino sources, such as active galactic nucleus (AGN) jets (e.g., Böttcher et al. 2013;Becker Tjus et al. 2014), γ-ray bursts (e.g., Waxman 1995), tidal disruption events (Stein et al. 2021;Reusch et al. 2022;Jiang et al. 2023), starburst galaxies (e.g., Loeb & Waxman 2006), as well as galaxy clusters (e.g., Murase et al. 2008).Intriguingly, evidence for neutrino emission from a nearby Seyfert galaxy NGC 1068 has been revealed (IceCube Collaboration et al. 2022).However, because of the relatively large angular uncertainties of neutrinos, the origin of extragalactic astrophysical neutrinos remains unresolved.
In addition to spatial information, temporal information is also crucial for determining the origin of neutrinos.IceCube has detected a 0.3 PeV muon neutrino (IC-170922A) cospatial with a blazar TXS 0506+056 (z = 0.3365; Paiano et al. 2018), at the time when the blazar experienced unprecedented multiwavelength flares, yielding a 3σ significance for the correlation of observed multimessenger radiations (IceCube Collaboration et al. 2018a).Blazars are a rare subclass of AGNs that include flat-spectrum radio quasars (FSRQs) and BL Lacertae objects (BL Lacs), of which strong relativistic jets are well aligned with our light of sight (Blandford et al. 2019).Blazars, which are the dominant population in the extragalactic γ-ray sky, exhibit overwhelming and highly variable Dopplerboosted nonthermal jet emissions that cover a wide range of energy regimes from radio to very high-energy γ-rays (Madejski & Sikora 2016).There are a growing number of cases in which both spatial and temporal coincidences have been observed between flaring blazars and incoming neutrinos (e.g., Kadler et al. 2016;Garrappa et al. 2019;Franckowiak et al. 2020;Giommi et al. 2020;Liao et al. 2022;Sahakyan et al. 2023).However, such circumstances are still rather rare.Furthermore, from an observational perspective, the connections between electromagnetic emissions and neutrinos are rather intricate.Electromagnetic activities were not observed during neutrino flaring of TXS 0506+056 in 2014-2015(IceCube Collaboration et al. 2018b).However, one of the brightest γ-ray blazars, PKS 1502+106, lies within the localization uncertainty region of the neutrino event IC-190730A (e.g., Rodrigues et al. 2021).It was in a quiescent state from optical emission to γ-rays, but with high flux-level radio emission, when the neutrino arrived (Kun et al. 2021).
On 2020 December 21, a 0.174 PeV track-like neutrino event (IceCube-201221A, hereafter IC-201221A) with a 56.4% probability of being an astrophysical origin (i.e., Gold alert streams) with the arrival direction given as R.A. was reported (IceCube Collaboration 2020).Subsequent electromagnetic observations were triggered, but no corresponding activities were reported (e.g., Fernandez & HAWC Collaboration 2020;Garrappa et al. 2020).However, the fourth data release of the Fermi γ-ray source catalog (4FGL-DR4, Ballet et al. 2023) includes a listing for a γ-ray source, 4FGL J1718.5+4237, in this direction.However, no corresponding low-energy counterpart is provided.Given the absence of a known source in the 4FGL-DR3 catalog (Abdollahi et al. 2022), the detection of 4FGL J1718.5+4237 is likely attributed to its recent flaring activity, which suggests a possible temporal correlation with IC-201221A.In light of these observational findings, we thoroughly examine 4FGL J1718.5+4237 in this Letter, including the search for its low-energy counterpart, and explore its potential association with the neutrino event IC-201221A.Analyses of multiwavelength data, including an optical spectroscopic observation by Hale 5 m telescope, along with the results, are presented in Section 2; while discussions and a short summary are given in Section 3. We adopt a ΛCDM cosmology with Ω M = 0.32, Ω Λ = 0.68, and a Hubble constant of H 0 = 67 km −1 s −1 Mpc −1 (Planck Collaboration et al. 2014).

Fermi-LAT Data
We analyzed the first 14.4 yr (i.e., from 2008 August 4 to 2023 January 1) Fermi-LAT Pass 8 data (evclass = 128 and evtype = 3) between energy 100 MeV and 500 GeV.We used Fermitools software version 2.0.8 and the accompanying Fermitools-data version 0.18.A region of interest (ROI), centered on 4FGL J1718.5+4237, with a 10°r adius was considered.In order to avoid significant contamination from the earth limb, a zenith angle cut (<90°) was adopted, so as the quality-filter cuts (i.e., DATA_QUAL==1 && LAT_CONFIG==1).Unbinned likelihood analysis was carried out using the gtlike task to extract the flux and spectrum of γ-rays.The significance of detecting a γ-ray source was determined by the test statistic (TS; TS = −2ln(L 0 /L); Mattox et al. 1996).L and L 0 denote the maximum likelihood value for the model with and without the target source, respectively.The parameters of the background sources within the ROI as well as normalizations of the two diffuse templates (i.e., gll_iem_v07 and iso_P8R3_SOURCE_V3_v01) were left free, while other parameters were fixed to values listed in 4FGL-DR4.We examined whether new γ-ray sources emerge (i.e., TS > 25) by analyzing subsequently generated residual TS maps.If so, we then reran the likelihood fitting process with the updated background model.In the process of temporal analysis, low-intensity background sources (TS < 10) were excluded from the model file.Additionally, the pyLikelihood UpperLimits tool was employed to compute a 95% confidence level (CL) upper limit instead of a flux, in cases of tentative detections.
First, a complete analysis of the entire data set was performed, from which the detection of a significant γ-ray source (TS = 60) in this direction was confirmed.The average flux of γ-rays is given as (4.05 ± 1.07) × 10 −9 ph cm −2 s −1 , along with a relatively soft spectrum (Γ = 2.57 ± 0.19, where Γ is the index of the power-law function, dN/dE ∝ E −Γ ).The optimized location was constrained as R.A. 259°.642,decl.
42°. 6215, with a 95% CL error radius of 5¢.These results are consistent with those listed in 4FGL-DR4.Within the radius of localization error, two radio sources, NVSS J171822+423948 (hereafter cand.I) and NVSS J171853+423658 (hereafter cand.II), are found.Therefore, the association relationship cannot be determined by spatial information alone.Investigations of two candidates were carried out.The former has been identified as a flat spectral radio source (Myers et al. 2003), while for the latter, data only from the NRAO Very Large Array sky survey (NVSS) are accessible (Condon et al. 1998).The Sloan Digital Sky Survey (SDSS) reveals a possible dropout in the u band for cand.I (Ahumada et al. 2020), suggesting that it may be a high redshift quasar.Based on the SDSS g-band magnitude and radio flux density at 5 GHz, radioloudness of cand.I is given as high as ∼1500 (Kellermann et al. 1989), indicative of a blazar candidate.
A half-year time-bin γ-ray light curve was extracted to investigate its temporal behavior, as shown in Figure 1.A distinct feature of the light curve is the significant flux variation on a long timescale.Adopting the "variability index" test (Nolan et al. 2012), the significance of the variability is given as >5σ.
Considering the relatively limited angular resolution of Fermi-LAT for sub-GeV photons, light curves from nearby background sources were also extracted.Since no similar variations were found to those of the target, the variation of the target is suggested to be intrinsic rather than artificially caused by the neighbors.Furthermore, the Bayesian block analysis (Scargle et al. 2013) marks a 2 yr epoch (i.e., the last four time bins) in a high-flux state between MJD 59219 and 59945.A specific analysis yields a significant γ-ray source (TS = 136).The corresponding photon flux is (1.4 ± 0.24) × 10 −8 ph cm −2 s −1 , with a slightly harder spectrum (Γ = 2.43 ± 0.11) than the value derived from the analysis of the entire data set.No significant spectral curvature is found by sophisticated spectral templates.By utilizing the Fermi gtsrcprob tool, the most energetic photon is at approximately 10 GeV, while the majority of highenergy γ-ray photons possess an energy of around 2 GeV.The localization analysis provides an optimized position of R.A. = 259°.64and decl.= 42°.6165,with a 95% CL error radius of approximately 4¢.Both candidates remain within the error radius.This source had been maintained in a low-flux state for 12 yr since the beginning of Fermi-LAT observation.
Combining the low-state data produces only a tentative detection (TS = 15) of a flux (1.89 ± 0.58) × 10 −9 ph cm −2 s −1 , which is approximately 7 times lower than the peak value.
We extracted a zoomed-in 2 month time-bin light curve, focusing on the epoch of the high-flux state (see Figure 2).The γ-ray activity consists of two subflares, separated by a 1 yr long quiescent state.The first one lasts about 4 months (TS = 54), between MJD 59214 and 59336.The duration of the other flare is also several months (TS = 113).Individual analyses of data for these two subflares reveal a similar flux level, ;2 × 10 −8 ph cm −2 s −1 , approximately 1 order of magnitude of the flux level extracted from the first 12 yr Fermi-LAT data.The two subflares also share a similar spectrum (i.e., Γ ; 2.3).In addition, localization analyses provide similar results so that the two radio sources fall into the γ-ray radii.Weekly time-bin γ-ray light curves were further extracted, but no sign of variability at a timescale of days is found.

Hale P200 Spectroscopic Observation
The spectroscopic observation of cand.I was performed at MJD 60223, using the double spectrograph, with dichroic D55, mounted on the Hale 200 inch telescope at Palomar Observatory (P200; Oke & Gunn 1982).The total exposure time for the observation is 2400 s.To match the seeing conditions, we chose a slit width of 1 5. Data were reduced with Pypeit (Prochaska et al. 2020).We conducted flux calibration using a standard star on the same night.Galactic extinctions were corrected (Schlafly & Finkbeiner 2011).A simple power-law model was used to fit the continuum component after excluding the wavelength windows of the emission lines, the telluric region, and the blueward region of the peak of Lyα.
The spectrum shows strong Lyα, moderately weak C IV, and faint C III emission lines, which give a redshift of 2.68 ± 0.01, as shown in Figure 3.The C IV emission line was fitted with a single Gaussian.FWHM C IV is given as ∼4397 ± 1714 km s −1 , of which the uncertainty is derived by Monte Carlo simulations, also listed in Figure 3. Based on the estimated λ L 1350 ∼ 1.12 × 10 46 erg s −1 and the bolometric correction factor BC 1350 = 3.81 (Shen et al. 2011), the bolometric luminosity reaches to L bol,disk ∼ 4.28 × 10 46 erg s −1 .Then, the mass of the SMBH can be inferred as M M log 9 BH ( ) ☉ (Vestergaard & Peterson 2006).Therefore, the Eddington ratio r Edd ∼ 0.3 is suggested.Although the peaking frequency of the synchrotron emission for FSRQs is typically at the submillimeter/infrared domain, temporal jet contribution in optical bands may not be negligible when fast and largeamplitude optical variations associated with jet activities are exhibited.In order to investigate the influence of the jet emission in our optical spectrum, based on the relationship between L C IV and λ L 1350 of SDSS quasars (Shen et al. 2011), λ L 1350 is deduced from the C IV line luminosity, which is consistent with that directly obtained from the continuum spectral component.Therefore, the jet is likely at a quiescent state at the time of the spectroscopic observation and its contribution in optical bands is not significant then.In addition, a description of archival SDSS measurements (uband data excluded) that correspond to a low optical flux level (see Figure 4), by a multitemperature radial profile (Shakura & Sunyaev 1973), has been performed.A similar value of L bol,disk ∼ 3 × 10 46 erg s −1 is yielded, compared with the result from the optical spectrum.

Zwicky Transient Facility Light-curve Data
We collected the latest Zwicky Transient Facility (ZTF) data (Masci et al. 2019; PTF Team 2020). 5Based on coadded reference images in the g, r, and i bands, magnitudes within 5″ of cand.I were extracted to construct optical light curves, shown in Figure 2.Only frames satisfied with catflag=0 and chi < 4 were selected.Optical activities in multiepoch, for instance, around MJD 58600, 59300, and 59700, are exhibited.The variability amplitudes in all three bands are approximately 1 magnitude.For cand.II, no significant excesses against the background emission were found from the reference images in all three bands.

Wide-field Infrared Survey Explorer Data
We investigated the infrared counterparts of the two candidates.An infrared source, Wide-field Infrared Survey Explorer (WISE) J171822.76+423945.1, is clearly detected at the position of cand.I, while no sources were found within 10″ of the radio position of cand.II, according to the AllWISE Source Catalog (Wright et al. 2010;WISE Team 2019).The W1 (3.4 μm) and W2 (4.6 μm) magnitudes given by the catalog are 16.587 ± 0.059 and 16.152 ± 0.123, respectively.We first attempted to construct the long-term light curves after excluding bad exposures and binning the data within each epoch (within 1 day) following Jiang et al. (2021a).Although with larger errors (typically ;0.2 mag), there are two obvious flux maxima at MJD 59274 (i.e., W1 = 15.412 ± 0.294 mag and W2 = 14.439 ± 0.317 mag) and MJD 59804 (i.e., W1 = 15.297 ± 0.166 mag and W2 = 14.411 ± 0.26 mag).
To confirm this and obtain more accurate measurements, we also attempted to perform point-spread function photometry on the differential images of time-resolved WISE/NEOWISE coadds (Meisner et al. 2018) with the first epoch as a reference following Jiang et al. (2021b).Finally, we adopted the fluxes measured by subtraction of the image with the reference flux added (see Figure 1).

Implications of Multimessenger Observations
Blazars are well-known for significant broadband variability, on timescales from minutes to years (Wagner & Witzel 1995;Ulrich et al. 1997).A remarkable observational feature of FSRQs is that correlated optical/infrared and γ-ray flares are frequently detected (e.g., Abdo et al. 2010).Since the angular resolution of Fermi-LAT (typically ;0°.1) is limited compared to those of the optical and radio bands (3″), the observed correlated emissions are essential to determine the association relationship between γ-ray source and its low-energy counterpart.For the two radio candidates falling into the localization uncertainty area of 4FGL J1718.5+4237,optical and infrared temporal behaviors have been investigated here.Significant flux variations of cand.I have been detected, while no such evidence is found for cand.II.Furthermore, contemporaneous brightenings between optical/infrared emissions of cand.I and γ-ray emissions of 4FGL J1718.5+4237 are revealed.As shown in Figure 2, corresponding to a γ-ray activity around MJD 59300, distinct optical flares in all three independent bands are presented.Meanwhile, its infrared emission is at a high flux state.Similar behaviors are found for another γ-ray activity MJD 59800, though there are no available i-band observations then.Based on these observed features, in conclusion, cand.I (i.e., NVSS J171822+423948) is suggested as the low-energy counterpart of 4FGL J1718.5+4237.In the latest incremental version of the fourth catalog of AGNs detected by Fermi-LAT (Ajello et al. 2022), only 29 known γray emitting blazars beyond the redshift of 2.6 are embraced, among more than 3000 sources.Our optical spectroscopic observation reveals the redshift of cand.I as z = 2.68 ± 0.01, and hence it is noticed as one of the rare high redshift γ-ray emitting blazars, which are crucial for studying relativistic jets and mass growth of the central SMBHs at the early epoch of the evolution of the Universe (e.g., Volonteri 2010;Ackermann et al. 2017;Liao et al. 2018).
The contemporaneous brightening of cand.I suggests that its optical and γ-ray emission is likely from the same radiation region.Although no sign of fast variability in γ-ray is found, optical flares peaking at MJD 59319 put a constraint on the radiation region.For instance, the i-band flux increases 1.1 mag within 10 days, giving a doubling time, in cosmic rest frame,  t F F z ln 2 ln 1 1.9 doub 1 2 t = D ´+ ( ) ( ) days.A similar value is obtained from the flux variations in the r band, but the timescale in the g band is relatively longer.Hence, adopting a routine Doppler factor (δ) value of 10 (Liodakis et al. 2017), the radius of the radiation region is limited to r  τ doub cδ ; 0.02 pc.Assuming a conical jet geometry, the location of the radiation region is constrained to R loc ∼ δr ∼ 0.2 pc, while the radius of the broad-line region (BLR) is also suggested as ∼0.2 pc based on the detected λ L 1350 (Grier et al. 2019).Therefore, the radiation region is likely embedded within the radiation of BLR.Note that the observational cadence of ZTF is several days, so potential intraday optical variability could be missed.In fact, fast γ-ray variability for blazars of z  3 with a timescale down to a few hours has been reported (Li et al. 2018).Meanwhile, intraday infrared variability has also been found in one of these sources (Liao et al. 2019).All of these pieces of evidence indicate the presence of a compact dissipation region.
Benefiting from multimessenger observations, the connection of the flaring blazar and the neutrino event is investigated.First, cand.I is the only known γ-ray source spatially coincident with IC-201221A.However, due to its relatively large localization uncertainty, relying solely on spatial coincidence is not sufficient to claim a physical connection.Monte Carlo simulations have been performed to calculate a probability of chance.The central location of IC-201221A is randomized in the sky beyond the Galactic plane with a varying  (Ballet et al. 2023) that fall into the simulated uncertainty box.Obviously, additional temporal information is needed.As shown in Figure 1, the incoming of the neutrino seems to switch on γ-ray emission of the blazar.It maintains to be in a quiescent state lasting for more than 12 yr beneath the detection limit of Fermi-LAT.However, it becomes detectable soon after the arrival of IC-201221A, roughly several tens of days.The amplitude of variability between different flux levels is 1 order of magnitude.Since the yearly time bin γ-ray light curve for each 4FGL-DR4 source is available,6 the simulations are reperformed including temporal information.At this time, only sources exhibiting a clear flux increase in the 13th time bin (i.e., 3 times of the 14 yr averaged flux) are considered.The probability of chance reduces to ∼0.03 from ;1.A similar case in which the emergence of a cospatial transient γ-ray source coincides temporally with neutrino detection is reported for GB6 J2113+1121/IC-191001A (Liao et al. 2022).Apart from a long period in a quiet state (i.e., over 8 yr), γ-ray flux enhancement about 5 times of TXS 0506+056 in long-term are found, corresponding to the incoming of IC-170922A (IceCube Collaboration et al. 2018a).Furthermore, in a multiwavelength perspective, strong optical flares with a peaking time 115 days after the arrival time of the neutrino (i.e., MJD 59204) are detected, although prior optical activities are exhibited, see Figure 2. In addition, accompanying infrared flux increases at MJD 59274 are also presented.Conclusively, cand.I is suggested as a suggestive neutrino emitter.

Discussions and Summary
It is interesting to put some general theoretical constraints on the neutrino production of cand.I based on the multimessenger data.The number of muon (antimuon) neutrinos detected by IceCube during a time interval ΔT at a decl.δ is given as ,max PeV represent the limits of the energy range that one expects 90% of neutrinos in the gamma-ray follow-Up (GFU) channel (Oikonomou et al. 2021), while f n m denotes the differential flux of muon neutrinos.For IC-201221A, the GFU_Gold effective area of IceCube is approximately A eff ≈ 6 m 2 (Abbasi et al. 2023).Neutrinos are assumed to follow a power-law spectrum (ò − γ with γ = 2).Corresponding to the duration of the high γ-ray flux state obtained by the Bayesian Blocks approach, ΔT is set to 2 yr.Thus, the integrated muon neutrino energy flux is derived as ∼1.6 × 10 −10 erg cm −2 s −1 , corresponding to the detection of one neutrino.Taking into account the high redshift of the blazar, the integrated muon neutrino luminosity is as high as  n m ∼ 10 49 erg s −1 , corresponding to an averaged muon neutrino luminosity   m /In(8 PeV/80 TeV) ∼ 2 × 10 48 erg s −1 .Alternatively, ΔT is set to 14 yr, representative of the time length of the entire Fermi-LAT data set.In this case, the integrated muon neutrino energy flux is reduced to 2.2 × 10 −11 erg cm −2 s −1 , and so is the averaged muon neutrino luminosity, ∼3 × 10 47 erg s −1 .Photopion (p + γ ⟶ p + π) processes are adopted as the responsible neutrino production mechanism, considering the luminous accretion disk emission of the blazar.The energy of the detected neutrino ò ν ; 0.6 PeV in the cosmic rest frame, suggests that the energy of the emitting proton is ò p ≈ 20ò ν ≈ 12 PeV.Quantities with subscripts "obs" refer to the observer frame, while those with primes denote the frame comoving with the jet, whereas other quantities are measured in the cosmic rest frame unless specified otherwise.Following estimates in Murase et al. (2018), ò ν ≈ 0.6 PeV (0.1 keV/ t ¢ )(δ/ 10).Considering that for FSRQs external photon fields are likely responsible for the target photons of photopion processes, the corresponding energy in the cosmic rest frame is inferred, ò t =  t d ¢ ∼10 eV, where  t ¢ is the energy in the jet comoving frame.The UV emissions could be radiations from BLR.During the interaction, protons lose three-eighths of their energy to neutrinos, resulting in an all-flavor neutrino luminosity given by   , where the optical depth to pπ processes ( f pπ ) is used to describe the efficiency of neutrino production.The remaining five-eighths of the lost proton energy contributes to the production of electrons and pionic γ-rays.Subsequent electromagnetic cascades start and synchrotron cooling is dominated because of Klein-Nishina suppression in inverse Compton processes.The connection between neutrino radiations and the cascade is given as ¢ Gauss)(ò p /12 PeV) 2 (10/δ)(3.68/(1+ z)), in which B 4 ¢ = Gauss is taken.Note that the magnetic field strength of a few Gauss is typical for FSRQs, as suggested by spectral energy distribution (SED) modeling studies (e.g., Ghisellini et al. 2010;Rodrigues et al. 2024).Thus, a muon neutrino luminosity (antimuon) of 2.1 × 10 46 erg s −1 is anticipated.Compared to the luminosity of muon neutrinos inferred by the detection of IC-201221A, the Poisson probability of detecting such a neutrino is ≈0.01.In the case of ΔT = 14 yr, a   L g g , at ò syn,obs of 2 GeV, is ∼1.4 × 10 46 erg s −1 (Ballet et al. 2023), which gives an expected muon (antimuon) neutrino luminosity of 5.6 × 10 45 erg s −1 .The corresponding neutrino detection probability is suggested as ≈0.02.An upward fluctuation (∼2σ) is needed if the neutrino is radiated from the blazar, which is consistent with the value found from TXS 0506+056/ IC-170922A (e.g., Cerruti et al. 2019;Gao et al. 2019;Petropoulou et al. 2020).
The proton luminosity of the jet, derived by neutrino luminosity, is compared with the Eddington luminosity of the accretion disk, 1.3 × 10 47 erg s −1 for a black hole mass of 10 9 M ☉ .In the case that BLR photons (ò t ∼ 10 eV) are target photons, the optical depth to pπ processes is suggested as s » ṕ -cm 2 is the effective cross section (Murase et al. 2016).Meanwhile, the accretion-disk luminosity of NVSS J171822+423948 is obtained as L AD = 4.28 × 10 46 erg s −1 .Thus, f pπ ≈ 6.3 × 10 −2 is given.For ΔT = 2 yr,     L L p p = n n /(3f pπ /8) ≈ 2.7 × 10 48 erg s −1 .Assuming a powerlaw spectrum for protons ( )≈ 3.8 × 10 47 erg s −1 (Sahakyan et al. 2023), with a set of Γ = 10.Such a value is approximately 3 times higher than the Eddington luminosity.The required proton power exceeds the Eddington luminosity, which is consistent with the recent studies on the neutrino candidates (e.g., Cerruti et al. 2019;Gao et al. 2019;Rodrigues et al. 2021).Alternatively, if ΔT = 14 yr, the corresponding proton luminosity of the jet  p,jet ˆis ∼9.9 × 10 46 erg s −1 , which is comparable to the Eddington luminosity.
Besides TXS 0506+056, an intermediate synchrotron peak BL Lac (ISP; Padovani & Giommi 1995;Abdo et al. 2010), several other BL Lacs have been reported as neutrino-emitting blazar candidates, including GB6 J1040+0617 (low-synchrotron peaked blazars, LSP; Garrappa et al. 2019), MG3 J225517 +2409, and PKS 0735+178 (intermediate-synchrotron peaked blazars, ISP; Franckowiak et al. 2020;Sahakyan et al. 2023), as well as 3HSP J095507.9+355101(high-synchrotron peaked blazars, HSP; Giommi et al. 2020).Note that an accompanying X-ray flare of 3HSP J095507.9+355101,not in γ-rays, is detected as responsible for the arrival of a neutrino (Giommi et al. 2020).FSRQs are also proposed as potential neutrino emitters, such as PKS B1424-418 and GB6 J2113+1121 (Kadler et al. 2016;Liao et al. 2022).Interestingly, a minor γray flare of 1H 0323+342, an HSP radio-loud narrow-line Seyfert galaxy, is detected as temporally coincident with a cospatial neutrino (Franckowiak et al. 2020).Because of the strong optical emission lines, cand.I is classified as an FSRQ.Meanwhile, the shape of its broadband SED suggests that it is an LSP, see Figure 4. Accretions of the central SMBHs of FSRQs and RLNLS1s are efficient, leading to intense radiation fields external to the jet.If such photons act as target photons in photopion processes, rather than self-radiation from the jet alone, the production of neutrinos could be significantly enhanced.On the other hand, dense soft photons also likely play an important role in γγ absorption processes, and hence severe absorption creates a soft γ-ray spectrum and weakens the observational connection between γ-ray emissions and neutrinos (Kun et al. 2021).For LSP/ISP BL Lacs, modest emissions from the accretion system are anticipated, and serious dilution from the jet emission is an obstacle to detecting potential broad-line emissions (e.g., Padovani et al. 2019).
In summary, a thorough investigation of 4FGL J1718.5 +4237 has been performed, which is the unique known γ-ray source within the angular uncertainty region of IC-201221A.First, the low-energy counterpart of the γ-ray source is determined.Although two radio sources fall into its localization radius, contemporaneous brightening of cand.I (i.e., NVSS J171822+423948) between γ-ray and optical/infrared emissions are detected.Meanwhile, no optical and infrared sources are found toward the other radio source.Therefore, cand.I is suggested to be associated with 4FGL J1718.5+4237.Moreover, optical spectroscopic observation suggests that cand.I is a rare high redshift (z = 2.68 ± 0.01) γ-ray emitting blazar.The multiwavelength temporal information is also important to establish a potential physical connection between the neutrino and the blazar.It remains in a quiescent state for 12 yr undetected by Fermi-LAT, but soon (i.e., a few tens of days) after the arrival of the neutrino, a tenfold γ-ray flux increase emerges.Furthermore, 115 days after the arrival time of the neutrino, strong optical flares in ZTF g, r, and i bands are exhibited; meanwhile, infrared fluxes are at a high level then.Monte Carlo simulations, with both spatial and temporal information embraced, give a by-chance probability of ∼0.03, suggesting that the blazar is a possible neutrino emitter and the first candidate above the redshift of 2. Theoretical constraints of neutrino production and comparisons with other candidates have been discussed.Future multiwavelength campaigns, especially X-ray observations in which contributions of cascades of Bethe-Heitler pairs may be nonnegligible, together with detailed theoretical SED modeling, would further investigate its jet properties and connection to the neutrino.

Figure 1 .
Figure 1.Upper panels: smoothed γ-ray TS maps (6°× 6°scale with 0°. 1 per pixel, 4FGL J1718.5+4237not included in the background model).The left one is based on the Fermi-LAT data between MJD 54683 and MJD 59204 (prior to the arrival of the neutrino).While the right one corresponds to the Fermi-LAT data for the rest of the time.The green X-shaped symbol and rectangle represent the optimized position and positional uncertainty of the neutrino, respectively.The blue circle is the 95% CL γ-ray localization error radius of 4FGL J1718.5+4237, and the violet cross marks the radio location of cand.I. Bottom panel: γ-ray light curves of 4FGL J1718.5+4237 and infrared ones of cand.I.In the former, blue circles and red triangles correspond to flux estimations and upper limits, along with TS values displayed as red bars.The gray shaded area represents an epoch of high γ-ray flux state obtained by Bayesian blocks.The black vertical dashed line across the light curves marks the arrival time of IC-201221A.

Figure 2 .
Figure 2. Zoomed-in multiwavelength light curves.Upper panel: 2 month time-bin light curve γ-ray light curves of 4FGL J1718.5+4237.Middle panel: ZTF light curves of cand.I, solid markers represent the magnitudes of the target, while hollow ones correspond to the average values of the comparison stars in the same field.Bottom panel: WISE light curves of cand.I.A black vertical dashed line across all panels marks the arrival time of the neutrino.

Figure 3 .
Figure3.P200 The spectrum is depicted in blue, the power-law continuum in green, and the Gaussian fitting for the C IV emission lines in orange.The red shadows represent the positions of the dichroic or the uncorrected telluric.Fitting results for the emission lines, including broad and narrow components of the Lyα line, as well as C IV and C III] lines, are listed.

Figure 4 .
Figure 4.The broadband SED of cand.I.The red, green, and purple points represent data in the 4 month flare state, while the blue points correspond to data in the entire 2 yr flare state.The orange parallelogram represents the Fermi-LAT data before IceCube-201221A, while the gray shadowed butterfly is for the 14 yr Fermi-LAT spectrum.The horizon lines correspond to energy-averaged neutrino energy flux for different epochs.Meanwhile, the X-shaped symbols are theoretical anticipations of neutrinos.
IC represents the Compton-Y parameter for pairs from the cascades (typically 1;Murase et al. 2018).Optimistically, cascade emission is capable of contributing significantly in the γ-ray domain.Since the majority of the observed high-energy γ-ray photons are around 2 GeV,   L g g ∼ 5.3 × 10 46 erg s −1 , at ò syn,obs of 2 GeV, is extracted from the observed γ-ray spectrum in the 2 yr high-flux state.On the other hand, theoretically,  p syn,obs p is set as the same energy.Assuming δ = 10 and UV emissions as target photons ( 0 = 18 eV, where m p c 2 is the proton rest energy), as well as δ = 10, eV/δm p c 2 ) ≈ 5 × 10 49 erg s −1 .Therefore, in the AGN frame, the proton luminosity of the jet  (Ghisellini & Tavecchio 2008),BLR photons in the AGN frame, r BLR stands for the BLR radius, f COV = 0.1 represents the covering factor(Ghisellini & Tavecchio 2008),