The strange case of the transient HBL blazar 4FGL J1544.3-0649

Blazars are a special type of active galactic nuclei (AGN) hosting a relativistic jet that happens to be pointing very close to the observer’s line of sight (e.g. Blandford & Rees, 1978; Padovani et al., 2017). The emission in these sources is dominated by non-thermal radiation across the electromagnetic spectrum, often exhibiting large variations on timescales ranging from years to a few hundred seconds. Blazars are usually grouped into BL Lacertae (BL Lac) objects and Flat Spectrum Radio Quasars (FSRQs), based on their properties in the optical band: the optical spectrum of FSRQs shows broad emission lines while that of BL Lacs is featureless, shows narrow lines or features of the host galaxy (Urry & Padovani, 1995).
Blazars are the most luminous and energetic persistent emitters in the known Universe, dominating in the extragalactic $\gamma$-ray sky. Among the total $\sim 5,800$ sources in the Fermi Large Area Telescope (Fermi-LAT) fourth source catalog of $\gamma$-ray sources (4FGL-DR2 ; Abdollahi et al., 2020; Ajello et al., 2020a; Ballet et al., 2020; Lott et al., 2020) $\sim 3,100$ are blazars, 47 are radio galaxies, and 19 are other types of AGNs. Interestingly, $\gamma$-ray emission has been detected also from blazars at very high redshifts (e.g., Ackermann et al., 2017; Sahakyan et al., 2020) which enables the study of the evolution of the most luminous relativistic jets over cosmic time.
The broadband spectral energy distribution (SED) of blazars has a characteristic double hump shape. The low energy peak (radio to UV/X-ray) is usually explained by synchrotron emission of relativistic electrons propagating along the jet, although recent models developed to explain the emission of high-energy neutrinos in blazars assume that this component may be due to proton synchrotron radiation, at least during flares (Mastichiadis & Petropoulou, 2020). The second SED component extends to High Energies (HEs; $>100$ MeV), its origin being debatable. In the leptonic scenarios, this component is interpreted as inverse Compton upscattering of synchrotron photons (Synchrotron self Compton [SSC] Maraschi et al. (1992); Bloom & Marscher (1996); Ghisellini et al. (1985)) or of external photons (external inverse Compton) (Sikora et al., 1994; Dermer & Schlickeiser, 1994). On the other hand, in hadronic models the HE component is mainly due to proton synchrotron emission (Aharonian, 2000; Mücke & Protheroe, 2001) or pion decay (Mannheim, 1993; Mannheim & Biermann, 1989; Mücke & Protheroe, 2001; Mücke et al., 2003; Böttcher et al., 2013). In the latter case, blazars are also sources of Very High Energy (VHE; $>100$ GeV) neutrinos (IceCube Collaboration et al., 2018a, b; Ansoldi et al., 2018; Keivani et al., 2018; Murase et al., 2018; Padovani et al., 2018; Sahakyan, 2018; Righi et al., 2019; Cerruti et al., 2019; Gao et al., 2019).
Depending on the peak of the low energy component ($\nu_{\rm peak-s}$) blazars can be further grouped into three subcategories, namely low-energy peaked blazars (LBLs; $\nu_{\rm peak-s}<10^{14}$ Hz), intermediate peaked blazars (IBLs; $10^{14}\>{\rm Hz}<\nu_{\rm peak-s}<10^{15}$ Hz) and high-energy peaked blazars (HBLs; $\nu_{\rm peak-s}>10^{15}$ Hz) (Padovani & Giommi, 1995; Abdo et al., 2010). Observations in the X-ray band reveal that sometimes the synchrotron peak of HBLs can reach energies of $\sim$1 keV, ($\sim 2\times 10^{17}$ Hz) or even larger, showing an extreme behaviour (extreme HBLs, e.g. Giommi et al., 1999; Costamante et al., 2001; Biteau et al., 2020). For example, during the flares of Mkn 501 the synchrotron peak reached $\sim 100$ keV (Pian et al., 1998) or during the flares of 1ES 1218+304 the X-ray spectral index hardened to $\Gamma\leq 1.80$ shifting the peak towards higher energies (Sahakyan, 2020). The third catalog of high energy peaked blazars (3HSP, Chang et al., 2019) includes several objects in this category. Recent observations in the HE and VHE $\gamma$-ray bands have revealed an additional class of BL Lac objects VHE $\gamma$-ray spectrum of which is characterized by a hard intrinsic photon index of up to 1 TeV after correction for the extragalactic background light (EBL) absorption effects (BL Lacs extreme in $\gamma$-rays; Bonnoli et al. (2015); Tavecchio et al. (2011)).
The BL Lacs having featureless optical spectra have been traditionally discovered in radio or X-ray surveys (e.g., Giommi et al., 2005). Currently, the $\gamma$-ray data from Fermi-LAT observations combined with the data at the lower energy bands (e.g., infrared, IR, Massaro et al., 2012) are a powerful additional method; nearly 38% of detected blazars are BL Lacs in the fourth catalog of AGNs detected by Fermi-LAT (Chang et al., 2019; Ajello et al., 2020b). Nearly all these $\gamma$-ray emitting BL Lacs have detected radio counterparts, which makes radio observations crucial also for the understanding of their physics. In this regard of a particular interest is the identification of radio-weak BL Lac objects, i.e., BL Lacs from which radio emission is detected in follow-up observations but which do not have counterparts in the major radio surveys. A handful observations of radio-weak BL Lac objects will have an impact on our understanding of AGN unification as well as on BL Lac associations based on the radio surveys (Massaro et al., 2017; Chang et al., 2019).
In May 2017 Fermi-LAT detected a new $\gamma$-ray source (4FGL J1544.3-0649) not associated with any previously known $\gamma$-ray object; it was bright in the $\gamma$-ray band in two consecutive weeks after the discovery (Ciprini et al., 2017). The X-ray follow-up observations by Neil Gehrels Swift Observatory (Gehrels et al., 2004), (hereafter Swift) found a new bright X-ray source at a position corresponding to the optical transient ASASSN-17gs detected at $V=17.3$ mag on May 25 (Ciprini et al., 2017). The MAXI team (Kawase et al., 2017) also reported a detection of this source on May 12-13 at a flux level similar to that observed by Swift, with a flux increase of a factor $\sim$2 in the period May 21 to May 25. The spectroscopic redshift of the host galaxy, 2MASX J15441967-0649156, was estimated to be $z=0.171$ using the MDM 2.4m Hiltner telescope (Chornock & Margutti, 2017). The source has a counterpart in the radio band with flux densities of 46.6 mJy at $1.4$ GHz (NRAO VLA Sky Survey; Condon et al., 1998), 35.8 mJy at 3 GHz (Very Large Array Sky Survey, epoch 1; Gordon et al., 2020) and 67 mJy at 150 MHz (TIFR GMRT Sky Survey; Intema et al., 2017). In addition, Bruni et al. (2018) monitored 4FGL J1544.3-0649 with the Effelsberg-100 m single dish radio telescope for a four-month follow-up after the brightening in the $\gamma$-ray band. These data imply a flat radio spectrum very similar to that of most blazars, and a radio luminosity of $\sim 5\times 10^{40}$ erg s^-1, at 1.4 GHz, more than a factor 10 larger than that of Mrk 501 and Mrk 421. Considering the post-burst data at 5 GHz, the radio to X-ray flux ratio is log (F${}_{\rm 1.4\>GHz}$/F${}_{\rm 1\>kev}$) = -4.3 typical of HBL Blazars. Follow up XMM-Newtown (XMM) observations showed that the source exhibits strong variability both in the X-ray flux and spectral shape on timescales of $\sim 10$ ks and the X-ray spectrum is described by a broken powerlaw (Ursini et al., 2019).
The transient nature of 4FGL J1544.3-0649 makes it a unique object and an interesting target for broadband studies. Understanding of the changes in the jet of 4FGL J1544.3-0649 which led to its brightening in the X-ray and $\gamma$-ray bands can shed light on the physics at work in the BL Lac jets. In this regard, the many monitoring observations of 4FGL J1544.3-0649 with Swift provide unprecedented data both in the optical/UV and X-ray bands while that of Fermi-LAT in the 100 MeV-500 GeV band allowing to shape the low and high energy components in the SED of 4FGL J1544.3-0649 and investigate its evolution in time. Also, this broadband data can be used for theoretical modeling allowing to derive the main parameters characterizing the jet of 4FGL J1544.3-0649.
The purpose of this paper is to investigate the origin of the broadband emission from 4FGL J1544.3-0649 and the astrophysical multi-messenger implications of the dual behaviour of this source. The paper is organized as follows. The Fermi-LAT and Swift data extraction and analyses are described in Section 2. In Section 3 the origin of the multiwavelength emission is investigated, and the discussion and conclusion are given in Section 4.

Multiwavelength observations of 4FGL J1544.3-0649  allow us to build the SEDs with simultaneous and archival data as well as to investigate the time evolution of SEDs. The $\gamma$-ray data for each adaptively binned interval has been combined with the Swift optical/UV and X-ray data, with the Effelsberg-100 m observations of 4FGL J1544.3-0649  (in the frequencies from 5 to 15 GHz Bruni et al., 2018) and with archival data from the VOU-Blazar (Chang et al., 2020) and the SSDC³³3https://tools.ssdc.asi.it/SED SED tools, to form the SED/light curve animation that is available here youtu.be/9feaNWi1RDs. Moving over each $\gamma$-ray light curve interval, the spectrum observed in different bands evolves showing the time evolution of the SED. The months before the large $\gamma$-ray flare on MJD 58000 (May 2017) the source is barely detected in the $\gamma$-ray band, i.e., the SED on MJD 57274.73 contain upper limits and large uncertainties in the flux estimation. Then the $\gamma$-ray emission substantially increases which is sometimes accompanied by increases in the other bands. The $\gamma$-ray emission from the source is still detectable with a high significance until MJD 58865.87, but in the last interval (MJD $58947.23\pm 81.36$) the detection significance drops to $TS=17$. To investigate the nature of the multiwavelength emission from 4FGL J1544.3-0649, data from the following periods are considered:

• Period 1 (P1): averaged $\gamma$-ray emission state on MJD $57938.35\pm 12.80$ with contemporaneous Swift UVOT/XRT (Obsid: 10145008) and radio data.

• Period 2 (P2): from MJD 57984.00 to MJD 57994.33 corresponding to the period when the highest X-ray flux was observed (Obsid: 10145018).

• Period 3 (P3): 105.06-days interval centered on MJD $58571.12$ when the source was in a low state in the $\gamma$-ray, X-ray and optical/UV bands (Obsid: 10145057).

The broadband SED of 4FGL J1544.3-0649  shows the usual double humped shape that can be interpreted within a common one-zone synchrotron/SSC model (Fig. 3). It is assumed that the broadband emission is produced by relativistic electrons in a spherical blob of radius $R$ with a uniform magnetic field $B$. The emitting region is seen in the direction close to the line of sight which moves with the bulk Lorentz factor of $\Gamma\approx\delta$. It is assumed that the emitting electrons have commonly used powerlaw with an exponential cut-off energy distribution in the form of $N^{\prime}_{\rm e}(\gamma^{\prime}_{\rm e})=N^{\prime}_{0}\>\gamma^{\prime-\alpha}_{e}\>Exp[-\gamma^{\prime}_{\rm e}/\gamma^{\prime}_{\rm cut}]$ for $\gamma^{\prime}_{\rm e}\geq\gamma^{\prime}_{\rm min}$ where $\gamma^{\prime}_{\rm min}=E^{\prime}_{e}/m_{e}\>c^{2}$ is the minimum electron energy. These electrons interacting with the magnetic field radiate via the synchrotron mechanism which explains the low energy component in the SED while the HE hump is described by up-scattering of the synchrotron photons by the same electrons.
The model includes seven free parameters to describe the emitting electrons ($\alpha$, $N^{\prime}_{0}$, $\gamma^{\prime}_{\rm min}$ and $\gamma^{\prime}_{\rm cut}$) and specify the emitting region ($\delta$, $B$ and $R$). In order to reduce the number of free parameters, we assume the emission is produced in a compact region with $R=5\times 10^{16}\>{\rm cm}$ which corresponds to flux changes on days time scales, as it has been found in the $\gamma$-ray band (adaptive bins in Fig. 1). In order to constrain $\delta$ and $\gamma^{\prime}_{\rm min}$ initially the SED on MJD $57938.35\pm 12.8$ (P1) with (quasi) simultaneous data ranging from radio to HE $\gamma$-ray bands is modeled, resulting in $\delta=25.6$ and $\gamma^{\prime}_{\rm min}=481.1$ which are typical values usually obtained for BL Lacs. These were fixed parameters when modeling the SEDs observed on MJD 57989.16 (P2) and MJD 58571.12 (P3) for which the low energy data are lacking. In order to account for the contribution of the host galaxy (see the excess in the IR/optical part of the SED on MJD $57938.35\pm 12.8$), in addition to synchrotron/SSC component also we adopted a template of giant elliptical galaxy. The fitting is performed with the open source package JetSet (Massaro et al., 2006; Tramacere et al., 2011, 2009).
The SED modeling results are shown in Fig. 3 and the corresponding parameters presented in Table 1. In the left panel of Fig. 3, the Swift UVOT and ASAS-SN data (downloaded from the ASAS-SN Sky Patrol web site Kochanek et al., 2017)⁴⁴4https://asas-sn.osu.edu/) allow to estimate or at least constrain the contribution of the host galaxy; the peak of the host galaxy component is at $6.13\times 10^{-12}\>{\rm erg\>cm^{-2}\>s^{-1}\>}$. When the SED in the quiescent state is considered (P3), the data can be reproduced when $\alpha=1.96$ and $\gamma^{\prime}_{\rm cut}=1.13\times 10^{5}$ and the magnetic field is $B=1.37\times 10^{-2}\>{\rm G}$. In this case low energy component peak is at $\sim 2.42\times 10^{16}\>{\rm Hz}$ which, however, shifts to higher frequencies when the active states of the source are considered; i.e., during MJD $57938.35\pm 12.8$ the X-ray data defines the peak to be at $\sim 10^{17}\>{\rm Hz}$ whereas during the hardest X-ray state the peak is shifted well beyond $\sim 10^{18}\>{\rm Hz}$. This shift changes also the maximum energy up to which the particles are effectively accelerated; $\gamma^{\prime}_{\rm cut}$ increases to $1.16\times 10^{6}$ and $6.82\times 10^{6}$ when the SEDs on P1 and P2 are considered, respectively. The modeling resulted in a similar powerlaw index and magnetic field for these active X-ray periods, namely $\alpha=2.25$ and $7.82\times 10^{-3}\>{\rm G}$ for P1 and $\alpha=2.19$ and $8.54\times 10^{-3}\>{\rm G}$ for P2. In all the considered periods the jet is particle dominated with $U_{\rm e}/U_{\rm B}>3000$ for P1 and P2 and $U_{\rm e}/U_{\rm B}\simeq 350$ for P3. However, $U_{\rm e}$ strongly depends on $\gamma^{\prime}_{\rm min}$ which is not well constrained, and when allowing $\gamma^{\prime}_{\rm min}$ to vary $U_{\rm e}/U_{\rm B}$ ratio can change. The jet power in the form of magnetic field and electron kinetic energy are calculated by $L_{B}=\pi cR^{2}\Gamma^{2}U_{B}$ and $L_{e}=\pi cR^{2}\Gamma^{2}U_{e}$, respectively, and are given in Table 1.

We have presented the results of the analysis of optical/UV, X-ray and $\gamma$-ray observations of 4FGL J1544.3-0649, a remarkable blazar showing transient-like high-energy multiwavelength emission, a behaviour that was never observed before in AGNs. This source, in fact, remained below the sensitivity limits of all the high-energy surveys and observations until May 2017, when it brightened to such a level to be detectable in short exposures by Fermi-LAT and by the MAXI X-ray sky monitor (Kawase et al., 2017).
During the brightening, the 2-10 keV X-ray flux strongly varies increasing from $(1.0-5.0)\times 10^{-12}\>{\rm erg\>cm^{-2}\>s^{-1}}$ to $(1.28\pm 0.05)\times 10^{-10}\>{\rm erg\>cm^{-2}\>s^{-1}}$. This variability, which is even larger when the flux is estimated in each Swift orbit, is clearly associated to a harder-when-brighter behavior, as also noted by Ursini et al. (2019) based on XMM observations. Such flux and photon index correlation is well known in HBL blazars (e.g. Giommi et al., 1990), and is usually interpreted by the relation between the acceleration and cooling processes. As shown in Kirk et al. (1998), when the injection decreases and particles are loosing energy by synchrotron cooling or any other cooling process that is faster at higher energies, the spectrum becomes steeper.
During the $\gamma$-ray brightening of the source, the X-ray spectrum modifies as well. The available optical/UV and X-ray data defines the peak of the low energy component to be within $\nu_{\rm peak-s}=10^{15}-10^{17}$ Hz, implying 4FGL J1544.3-0649 is a classical HBL. However, occasionally (e.g., on MJD 57898.70, 57989.16, 57998.92, 58009.40 and 58021.40, see the SED/Light curve animation) the source shows a synchrotron component peak above $10^{18}$ Hz which is more typical for extreme synchrotron BL Lacs. Yet, in all these periods the $\gamma$-ray photon index hardens as well, e.g., $1.74\pm 0.20$ on MJD 57898.70, $1.57\pm 0.18$ on MJD 57989.16 or $1.61\pm 0.19$ on MJD 58021.40, shifting the HE peak to VHEs. The peaks in the SED are determined by $\gamma^{\prime}_{\rm cut}\simeq\sqrt{\frac{\nu_{\rm c,\>peak}}{\nu_{\rm s,\>peak}}}$ implying in P1 and P2 when $\nu_{\rm peak}\geq 10^{18}$ Hz the second peak is around $2\times 10^{26}-2\times 10^{27}$ Hz ($1-10$ TeV) (solid line in Fig. 3). However, this peak and the flux is suppressed by Klein Nishina effects and EBL absorption, which for the distance of 4FGL J1544.3-0649 starts at $\sim 5\times 10^{25}$ Hz ($\sim 200$ GeV; dot dashed line in Fig. 3). Despite this reduction, during the bright $\gamma$-ray and X-ray states the source flux at 1 TeV is expected to be $\simeq 10^{-11}\>{\rm erg\>cm^{-2}\>s^{-1}}$, well above the sensitivity of the current ground-based Cherenkov telescopes. In addition, the hard $\gamma$-ray photon index observed in the MeV/GeV band and the applied SSC model predict also a very hard photon index in TeV band after EBL correction. So, 4FGL J1544.3-0649 mimics properties similar to those of extreme BL Lacs in $\gamma$-rays (Bonnoli et al., 2015; Tavecchio et al., 2011), which are an interesting sub-class not only because of their possible association with VHE neutrinos and cosmic rays (Padovani et al., 2016; Resconi et al., 2017) but also because they pose major difficulties for the understanding of particle acceleration and emission processes in blazar jets.
The multiwavelength light curve of Fig. 1, shows that the optical/UV flux increase is delayed compared to that in the X-ray and $\gamma$-ray bands. In general, lags in different energy bands are expected either from changes in the global parameters (magnetic field, Lorentz factor, etc.) or from the temporal evolution of emitting particles. For example, Bednarek & Wagner (2008) proposed that if the low energy component is produced when the blob is closer to the base of the jet and the HE photons are emitted at larger distances from the jet base, at which the blob has accelerated to a higher Lorentz factor, this could result in delays between the low energy and HE flares. Or, in complex jet geometry if the flare is caused by the collision between a relativistic shock wave and a stationary feature in the jet, it results in a complex variation of fluxes in different bands (Sokolov et al., 2004). However, a simple consideration of particle acceleration and cooling can naturally explain the observed delay in the multiwavelength light curve of 4FGL J1544.3-0649. The SED in different epochs is well reproduced by a one-zone synchrotron/SSC model with physically reasonable parameters. Within this model, the flare and modification of the spectra can be interpreted by an increase in the luminosity of electrons and $\gamma^{\prime}_{\rm cut}$; during the bright X-ray period the luminosity increased $\sim 4.9$ times whereas $\gamma^{\prime}_{\rm cut}$ $\sim 62$ times. This implies that the X-ray and $\gamma$-ray flare was most likely caused by injection of new freshly accelerated electrons. The cooling time of these electrons is

where $u_{\rm syn}$ is the synchrotron photon energy density in the jet frame. Considering the peak frequency of synchrotron emission is $\nu=1.2\times 10^{6}\>B\>\delta\>\gamma^{2}$, the cooling time of the electron emitting photons with frequency $\nu_{15}=\nu/10^{15}\>{\rm Hz}$ is $t_{\rm cooling}=2.7\times 10^{4}\>(\delta/1+z)^{0.5}\>B^{-1.5}\>\nu_{15}^{-0.5}$ s. Considering the parameters obtained for the P1, the delay between X-ray emission ($10^{18}$ Hz) and optical/UV ($5\times 10^{14}$ Hz) is of the order of $\sim 114$ days which is very similar to that observed. However, in order to investigate in details the time lags and thus to infer the underlying physics, much denser monitoring of 4FGL J1544.3-0649 in optical/UV band is needed.
Figure 4 plots the high-energy light-curves of 4FGL J1544.3-0649  spanning the entire operational period of the Fermi mission. The upper panel, reporting the $\gamma$-ray data, illustrate the transient nature of this source, which was not detected by the LAT instrument until late 2015, peaked in 2017-2018 and then faded. The lower panel reconstructs the long-term light-curve of 4FGL J1544.3-0649  in the X-ray band, combining the Swift-XRT measurements (blue points), with the XMM data from Ursini et al. (2019) (red points) and the expected X-ray flux level estimated from the flux in the $\gamma$-ray band, based on the average x-ray/$\gamma$-ray flux ratio (open grey circles). An upper limit derived from the non-detection of the source in the RASS X-ray sky survey, which was carried out in 1990-1991 (Voges et al., 1999), is also shown in the leftmost part of the plot. The reconstructed X-ray light-curve implies a rise time of approximately two years, and a decay of approximately three years.

The SED of 4FGL J1544.3-0649 (Fig. 5) assembled with Swift-XRT, Swift-UVOT, Fermi-LAT (green points) and archival data (grey points), shows a highly variable X-ray spectrum, changing its intensity by large factors and often peaking in the X-ray band. For comparison the average SED of the well known HBL blazar Mrk501, rescaled to match the radio emission of 4FGL J1544.3-0649, together with its X-ray spectra observed during its minimum and maximum intensity state, is plotted as light orange points in the same figure. The nearly identical shape of the radio to optical SEDs and the range of observed X-ray fluxes, confirms that 4FGL J1544.3-0649 during the active phase is a typical HBL blazar, with X-ray and $\gamma$-ray emission that are even stronger than that of Mrk 501. The upper limit from the RASS all sky survey (red arrow) sets the magnitude of the intensity variations in the soft X-ray band, which ranged from below $\approx 10^{-13}\>{\rm erg\>cm^{-2}\>s^{-1}}$ in 1990 to such a high value to place 4FGL J1544.3-0649 amongst the 20 brightest X-ray blazars.

The existence of transient (HBL) blazars was never considered in the literature before the discovery of 4FGL J1544.3-0649 and, consequently, never taken into account in population studies. However, if this dual behaviour is relatively common, it might have implications in blazar research, especially in high-energy and multi-messenger astrophysics, depending on how numerous transient blazars are and how frequently they enter a bright phase. We don’t know if 4FGL J1544.3-0649 represents the tip of the iceberg, or it is simply an isolated event. The discovery of 4FGL J1544.3-0649 is a clear selection effect. Similar transient sources would probably go unnoticed in case their top flux was only a factor of a few lower than that of 4FGL J1544.3-0649, or with a softer $\gamma$-ray spectrum, or even in cases where the $\gamma$-ray intensity remains low during strong X-ray flares associated to neutrino emission as in scenarios similar to that presented by Mastichiadis & Petropoulou (2020). The identification of similar objects would add crucial information to our knowledge of the population of extragalactic sources, especially in the high-energy domain. Areas of research that might be affected would be e.g. the determination of the LogN-LogS, luminosity function and cosmological evolution of HBL blazars, the identification of still unassociated $\gamma$-ray and VHE sources, and the contribution of discrete sources to the $\gamma$-ray and possibly high-energy neutrino cosmic backgrounds. A transient nature of some electromagnetic counterparts of astrophysical neutrinos would on one side increase the level of difficulty in the process of the association, but on the other it would strengthen the association in case X-ray or $\gamma$-ray flares are detected in association with the neutrino arrival. A recent study (Giommi et al., 2020b) reported a 3.23 $\sigma$ post trial excess of HBL Blazars in the error regions of a sample of 70 Ice Cube astrophysical neutrinos with good positional accuracy. If HBL blazars are indeed counterparts of high-energy neutrinos, some of them could be transient sources somewhat fainter than 4FGL J1544.3-0649 that would difficult to locate since they are not listed in catalogues and no full coverage of the still large Ice Cube neutrino positional uncertainty can be achieved with current X-ray imaging instruments. A full census of these so far elusive objects would be possible only when sensitive all-sky X-ray (or $\gamma$-ray) monitor detectors become available. Meanwhile a careful search for weak transients in the Fermi-LAT data on the positions of flat-spectrum radio sources with radio to infrared flux ratio similar to that of HBL objects could lead to the discovery of additional transient blazars and set limits to their duty cycle and space density.

We acknowledge the use of data, analysis tools and services from the Open Universe platform, the ASI Space Science Data Center (SSDC), the Astrophysics Science Archive Research Center (HEASARC), the Fermi Science Tools, the All-Sky Automated Survey for Supernovae (ASAS-SN), the Astrophysics Data System (ADS), and the National Extra-galactic Database (NED).
PG acknowledges the support of the Technische Universität München - Institute for Advanced Studies, funded by the German Excellence Initiative (and the European Union Seventh Framework Programme under grant agreement no. 291763) and support by the DFG Cluster of Excellence "Origin and Structure of the Universe".
NS acknowledges the support by the Science Committee of RA, in the frames of the research project No 20TTCG-1C015.
This work used resources from the EGI infrastructure with the dedicated support of CESGA (Spain).
We are grateful to M. Petropoulou, P. Padovani and R. Middei for useful comments.

All the Swift-XRT spectra, best fit values and other imaging analysis results will be made available as part of a publication presenting the analysis of XRT data of all the blazars frequently observed by Swift (Giommi et. al. 2020, in preparation), in the framework of the Open Universe for blazars program. The remaining data underlying this article will be shared on reasonable request to the corresponding author.