Revealing High-z Fermi-LAT BL Lacs Using Swift and SARA Data with Photometric Analysis

BL Lacertae (BL Lac) objects are a subclass of blazar, distinguished by their featureless optical spectrum. The featureless spectrum presents a challenge in measuring the redshift of the BL Lacs. In this paper, we measure the redshift of BL Lacs using the photometric dropout technique. The space-based telescope Swift and the ground-based SARA telescopes are employed to provide magnitudes in the uvw2,uvm2,uvw1,u,b,v,g′,r′,i′,andz′ filters. We observe 60 BL Lacs and report reliable redshift upper limits for 41 of them. We discover three new high-z BL Lacs (z > 1.3) at redshifts of 1.74−0.08+0.05 , 1.88−0.03+0.07 , and 2.10−0.04+0.03 , bringing the number of high-z BL Lacs found by this method up to 19. Discussions are made on the implications for the blazar sequence, the Fermi blazar divide, and the gamma-ray horizon based on an analysis of the 4LAC catalog and all high-z BL Lacs found with the photo-z technique.


INTRODUCTION
Active galactic nuclei (AGNs) are laboratories for extreme physical processes that generate high luminosities (Urry & Padovani 1995).They are supermassive black holes at the center of galaxies (Fabian 2008), powered by active mass accretion (Marconi et al. 2004).The classification of AGN is based on their orientation relative to our line of sight (Urry & Padovani 1995).Blazars are a subclass of jetted AGN oriented with the relativistic jet pointing along the line of sight with a viewing angle θ m < 1/Γ, where Γ is the bulk Lorentz factor (Blandford & Rees 1978;Marcotulli et al. 2017).Blazars produce high-energy gamma-ray emission and show high variability (Perlman 2013).There are two characteristic bumps in the spectral energy distribution (SED) of blazars: the bump at lower energies (infrared to X-ray) is theorized to originate from synchrotron emission produced by electrons, and the bump at higher energies (X-ray to γ-ray) is interpreted as being due to inverse Compton scattering off either the synchrotron photons (Maraschi et al. 1994) or a circumnuclear photon field (Dermer & Schlickeiser 1994).
Blazars can be classified into different groups by their spectral properties (Padovani & Giommi 1995).Flat spectrum radio quasars (FSRQs) are blazars with broad emission lines in their spectra (equivalent width > 5 Å), while BL Lacs exhibit no or weak emission lines (equivalent width ≤ 5 Å) (Maraschi et al. 1994;Ajello et al. 2013).The lack of prominent emission lines of BL Lacs makes it challenging to measure their redshifts.Indeed, around 40 % of BL Lacs in the Fourth LAT AGN Catalog Data Release 2 (4LAC-DR2) lack measured redshifts (Ajello et al. 2020;Lott et al. 2020b).Additionally, blazars can be classified according to the synchrotron peak frequency: low synchrotron peaked (LSP) blazars have a synchrotron bump at lower energies where ν S peak ≤ 10 energies where 10 14 ≤ ν S peak ≤ 10 15 Hz; high synchrotron peaked (HSP) blazars have a synchrotron bump at higher energies where ν S peak ≥ 10 15 Hz (Abdo et al. 2010).Measuring redshifts of distant sources is essential in studying their energetic jets as well as a means of measuring the extragalactic background light (EBL, Ackermann et al. 2012).The EBL is the cumulative radiation from star formation and super-massive black hole accretion integrated over all redshifts.(Domínguez & Ajello 2015;Desai et al. 2019a).Direct measurements of the EBL are difficult because of the strong zodiacal light and foreground emission from our Galaxy (Moralejo et al. 2017).However, blazars can be utilized to measure the EBL as an alternative method to the direct measurement.The annihilation between γ-ray photons from blazars and EBL photons leaves a distinct attenuation in the spectra of blazars.This attenuation can be used to constrain the EBL models and its evolution with redshift (Biteau & Williams 2015;Abdollahi et al. 2018).Moreover, the attenuation increases as the redshift increases, motivating us to search for high-redshift sources.
Photometric redshift measurement has been applied to powerlaw sources like GRBs (Tagliaferri et al. 2005;Krühler et al. 2011).Rau et al. (2012) unitized the photometric method to measure the redshifts of blazars by fitting SED templates to multi-band photometry obtained by simultaneous Swift/UVOT (Roming et al. 2005) and Gamma-Ray burst Optical/NearInfrared Detector (GROND, Greiner et al. 2007) optical imaging.Photons from blazars are absorbed by the neutral hydrogen along the line of sight with a rest-wavelength blue-wards of 912 Å (Lyman-limit) and 1216 Å (Lyman-α forest), leaving two characteristic breaks in the source spectra.The breaks are redshift dependent and, therefore, can be used to measure the photometric redshift.The photometric technique enabled Rau et al. (2012) to measure the redshifts for nine BL Lacs out of a sample of 103 from the Second LAT AGN Catalog (2LAC) (Ackermann et al. 2011).Six of them were newly discovered BL Lacs with high redshifts (z > 1.3).Later, using the same technique except simultaneous observations, Kaur et al. (2017) found photometric redshifts for five sources from a sample of forty from the Third Fermi Large Area Telescope Source Catalog (3FGL) (Acero et al. 2015), and Kaur et al. (2018) discovered two more high-z BL Lacs.Rajagopal et al. (2020) continued the work and found three more sources with high redshifts, bringing the number of high-z BL Lacs found by this method up to 16.
Here we continue this campaign, using the Neil Gehrels Swift Observatory (Gehrels et al. 2004), 1.0 m SARA-RM (Southeastern Association for Research in Astronomy at Roque de los Muchachos, La Palma, Spain), and 0.6 m SARA-CT (Cerro Tololo, Chile) telescopes (Keel et al. 2016) to obtain magnitudes in 10 filters from ultraviolet to near infrared (nIR) (uvw2, uvm2, uvw1, u, b, v, g , r , i , z ) of the 60 selected sources.The photometric redshifts are estimated by the method developed by Rau et al. (2012), but the UVOT and SARA observations are not always simultaneous.A flat cosmological model ΛCDM with H 0 = 73 km/s/Mpc, Ω m = 0.3, Ω Λ = 0.7 is used in the content of this work.The structure of this paper is organized as follows: Section 2 introduces the instruments and observational methods used.Section 3 describes the procedures and techniques in the data analysis.Section 4 describes the SED fitting technique.Section 5 reports the redshifts for the sources we analyzed, which is followed by Section 6 that discusses our findings and applications of high-z BL Lacs.Section 7 is the summary of this work.

OBSERVATIONS
A sample of 60 BL Lacs with no redshift is selected from the 3FGL (Acero et al. 2015) and Third Fermi-LAT Catalog of High-Energy Sources (3FHL, Ajello et al. 2017).Of the 60 selected BL Lacs, 12 were observed during cycle 131 , 162 , and 173 of the Neil Gehrels Swift Observatory Guest Investigator program.The rest were approved as Targets of Opportunity.The selected BL Lacs are observed in the UV-optical band with six filters (uvw2, uvm2, uvw1, u, b, v) on the UVOT instrument of Swift.The total exposure time of the UVOT instrument is ∼ 2000s for each source.The time assigned to each filter is weighted as uvw2 : uvm2 : uvw1 : u : b : v = 4 : 3 : 2 : 1 : 1 : 1 The selected targets are observed in the optical-nIR band from the ground using four SDSS filters (g , r , i , z ) installed on SARA-RM and SARA-CT.The exposure in every filter ranges from 15 to 40 minutes to ensure a good signal to noise ratio of the images.The time separation between the first and last filter is around 50 minutes to 1 hour and 30 minutes depending on the brightness of the sources.Table 1 shows the details of the observation of our sample.

SARA data analysis
We calibrate the images with bias, dark and flat frames using IRAF v2.61 (Image Reduction and Analysis Facility, Tody 1986): the bias frames remove the electronic background, the dark frames remove the current caused by thermal effect, and the flat frames correct pixel-to-pixel variations in the sensitivity of the CCD.The same IRAF software is also used to perform photometry.We use standard stars (Smith et al. 2006;Landolt 2009;Albareti et al. 2017) observed on the same night to calibrate the magnitudes of the sources.We measure the full width half maximum (FWHM) of the point-spread function for each night.A schematic representation of the source and background regions chosen for analysis are shown in Figure 1.
Figure 1.A circular region with radius of 3 times the FWHM is used as the source region.An annulus with inner radius of 5 times the FWHM and width of 3 times the FWHM is used as the background region.This setting makes the area of the background annulus ∼4 times the area of the aperture, which ensures good statistics on the background estimation.

Swift UVOT data analysis
We analyze the UVOT data following the standard pipeline from Poole et al. (2008).The images are calibrated using the remote Swift Calibration Database (CALDB) v.1.0.2 to remove bad pixels, correct sensitivities, apply dark frames, etc.A circular region of 5 is used as the source region, and 20 to 30 is used as the background region.We use the command UVOTSOURCE provided in HEASoft v.6.284 to extract the AB magnitudes from the images.If the image file contains multiple extensions, UVOTIMSUM is used to align and stack the extensions.Both SARA and Swift magnitudes are corrected for the Galactic foreground extinction (Kataoka et al. 2008).

Variability correction and cross-calibration
Blazars are highly variable in the optical band on timescales from minutes to days (Racine 1970;Carini et al. 1991;Urry 1993;Otero-Santos et al. 2022).Since both Swift and SARA telescopes observe the sources in sequence, we must account for the variability caused by the non-simultaneous observation as a part of the systematic uncertainty in our photo-z technique.Following Rau et al. (2012) and Rajagopal et al. (2020), we include a systematic uncertainty ∆ = 0.1 mag for both UVOT and SARA filters based on the changes in the observed magnitudes in the uvw2 filter.
Cross-calibration between UVOT and SARA filters is also performed following Rau et al. (2012).We assume that the SED of the sources of interest remains unchanged.We approximate the SED by a power law.Color terms g − r and b − g are fitted using a quadratic equation (Rajagopal et al. 2020 (1) Using Equation 1, we calculate the offset in the b filter.This offset is applied to all UVOT filters.Cross-calibrated magnitudes are reported in Table 2.
The FORTRAN code LePhare v.2.2 (Arnouts et al. 1999;Ilbert et al. 2006) is used to perform SED fitting and find the photometric redshifts.LePhare utilizes a χ 2 fitting technique to determine a best-fit SED model for each source, based on the observational data in 10 filters.Three different libraries are included during this fitting process.The first library contains 60 power-law (F λ ∝ λ −β , β ∈ [0, 3]) templates.The second library consists of templates of normal galaxies and AGN/galaxy hybrids (Salvato et al. 2008(Salvato et al. , 2011)), which are included to model host-galaxy-dominated sources.The third library is composed of various stellar templates to avoid false associations (Bohlin et al. 1995;Pickles 1998;Chabrier et al. 2000).
Table 1.Swift-UVOT and SARA Observations along with visual extinction values.a Best-fit or upper limit of the photometric redshift.
b Photometric redshifts with 1 σ confidence level.c χ 2 value for ten degrees of freedom.
e Spectral slope for the power-law model of the form f The fittings are not reported because they are not good (χ 2 = 0) due to upper limits for six or more filters.
Following the criteria simulated in Rajagopal et al. (2020), sources with z > 1.3 and E(B −V ) ≤ 0.30 have reliable measured photometric redshifts within |∆z/ (1 + z sim )| < 0.15 accuracy, where z sim is the input redshift for the simulation.Another criterion is the integral of probability density function P z > 90%, which indicates that the measured redshift is within 0.1(1 + z phot ) of the best-fit value.Reliable upper limits from the power-law model are provided for sources with P z ≤ 90% and χ 2 ≤ 30 for ten degrees of freedom.Specifically, we only report the photometric magnitudes for sources with E(B −V ) > 0.3.

RESULTS
The best-fit photometric redshifts or upper limits for the BL Lacs, along with the χ 2 and P z values for the power-law and galaxy template models are reported in Table 3.Three high-z BL Lacs are found from the analysis of sixty BL Lacs, while forty-one are reported with reliable upper limits.The SEDs in UVOT and SARA filters for the four high-z sources are presented in Figure 2.Under the criteria mentioned above (P z ≥ 90%, E(B − V ) ≤ 0.3, and z > 1.3), the three high-z sources are 3FHL J0301.4−5618,3FGL J1032.5+6623, and 3FGL J1419.5+0449.The redshifts of the three sources are 1.74 +0.05 −0.08 , 1.88 +0.07 −0.03 , and 6.DISCUSSIONS

Cosmic Gamma-ray Horizon
Very high energy (VHE) photons (>100 GeV) do not travel unimpeded.The extragalactic background light (EBL) impedes the travel of VHE photons via photon-photon interactions in the ultraviolet, optical, and infrared bands.Due to these interactions there is a redshift dependent opacity for high-energy photons in the universe (Finke et al. 2010;Domínguez & Ajello 2015).This opacity to VHE photons is evaluated by the cosmic gamma-ray horizon (CGRH, Domínguez et al. 2013).The CGRH is defined as the energy at which the optical depth is one (τ = 1), and it is a function of redshift.
The gamma-ray SEDs of the four high-z BL Lacs are shown in Figure 4. We fit the Fermi-LAT Fourth Source Catalog (4FGL, Abdollahi et al. 2020;Lott et al. 2020a) and 3FHL data with a power law attenuated by EBL absorption.The EBL absorption (Domínguez et al. 2011) is applied with a factor e −τ(E,z) in the fits.The redshift uncertainty determined by the χ 2 fitting is considered in the fits, which is demonstrated by the gray shaded area in the plots.The fits show that redshift uncertainty has a negligible effect on the source flux.Moreover, the effect of redshift uncertainty decreases with higher redshift.In addition, the attenuation caused by the EBL can be seen in the joint fits of the 3FHL and 4FGL data.Fossati et al. (1998) introduced the blazar sequence to find a unified model for all types of blazars.The blazar sequence is based on three relations among the physical properties of blazars (Prandini & Ghisellini 2022).The first is an anti-correlation between the rest frame synchrotron peak frequency (ν sy pk, rest ) and the luminosity (L sy pk ) at this peak frequency.FSRQs occupy the region with lower ν sy pk, rest and higher L sy pk .On the contrary, BL Lacs show higher ν sy pk, rest and lower L sy pk .However, this phenomenon could also be ascribed to selection effects (Giommi et al. 2012).BL Lacs with high ν sy pk, rest and high L sy pk could exist among those that do not have a redshift measurement available.Thus, measuring the redshift of BL Lacs is important to constrain the blazar sequence.The second relation is the anti-correlation between ν sy pk, rest and Compton dominance.Compton dominance is defined as the ratio of the inverse-Compton peak luminosity (L IC pk ) to L sy pk .BL Lacs and FSRQs populate different regions of this plane: FSRQs have higher Compton dominance compared to BL Lacs.The last relation is also an anti-correlation, but it is between the spectral index (Γ γ ) and ν sy pk, rest .They follow similar behavior as the second one mentioned above: FSRQs have larger gamma-ray index and populate lower synchrotron frequencies, while BL Lacs have smaller gamma-ray index and populate higher synchrotron frequencies.  z) due to EBL absorption.The EBL model is taken from Domínguez et al. (2011).The uncertainty region (the grey shaded area) shown for this joint fit is exclusively coming from the redshift uncertainty, according to Domínguez et al. (2011), and the gray area is not easily visible due to the small effects from the redshift uncertainty.The joint fit shows where the energy cutoff is.Note that the energy cutoff is solely due to the EBL attenuation.

Blazar Sequence
In order to explore these relations, we use the SSDC Sky Explorer5 to calculate the synchrotron peak frequency and flux of high-z sources if they are not provided in the 4LAC catalog.We also use the equation from Marcotulli et al. (2020) to calculate the inverse Compton peak flux for sources with power-law gamma-ray spectra: The break energy E b is estimated from the power-law photon index Γ using the E b − Γ relationship reported in Ajello et al. (2015).The normalization constant K is set to reproduce the 1-100 GeV flux reported in the 4FGL catalog.Figure 5(a) shows the existence of high-z luminous BL Lacs that are not detected easily via traditional optical spectroscopy but are discovered by our photometric redshift.On the other hand, in both Figures 5(b) and 5(c), our objects are more in line with the standard regions occupied by BL Lacs.After three months of Fermi-LAT observations (Abdo et al. 2009), Ghisellini et al. (2009) reported a clear divide between BL Lacs and FSRQs on the gamma-ray index (Γ γ ) and gamma-ray luminosity (L γ ) plane.BL Lacs are less luminous (L γ < 10 47 erg s −1 ) and have harder spectra (Γ γ <2.2) than FSRQs.Different accretion rates might account for this phenomenon (Ghisellini et al. 2009).Figure 5(d) reproduces the Fermi blazar divide plot using the 4LAC catalog (Ajello et al. 2020) and the high-z sources found by our photometric method.We calculate the gamma-ray luminosity from 0.1 GeV to 100 GeV using equation 1 in Ghisellini et al. (2009).All the high-z sources this paper find have a harder spectrum (Γ γ < 2.2).Although there is an overlap between BL Lacs and FSRQs, the separation of the two populations is still distinct.However, our BL Lac objects have properties in between FSRQs and BL Lacs (they have high luminosities and hard spectra).Giommi et al. (2013) proposed that these BL Lacs are "blue quasars" or "masquerading BL Lacs" because the overwhelming synchrotron emission outshines their broad emission lines.According to Sbarrato et al. (2012) the luminosity of the broad line region is connected to the gamma-ray luminosity by L BLR ∼ 4L 0.93 γ .On average we assume that 10% of the disk luminosity (L disk ) is processed by the broad line region: L BLR = 0.1L disk .Then the ratio of L disk to L Edd can be found.For the four high-z sources in this work, the ratios are 0.027, 0.023, and 0.007 for 3FGL J1032.5+6623,3FGL J1419.5+0449, and 3FHL J0301.4−5618respectively.Two of the three sources have the ratio between 0.02 to 0.05, suggesting efficient accretion as typical for FSRQs.The overall properties of the BL Lacs with high photometric redshifts (z > 1.3) allow us to characterize them as potential masquerading BL Lacs, thus increasing this interesting population.

SUMMARY AND CONCLUSIONS
In this work we use 10-filter photometry with Swift-UVOT and SARA to constrain the redshift of 60 BL Lac objects without a spectroscopic measurement.Our method, based on the method used in Rau et al. (2012), can produce a photometric redshift measurement for BL Lacs at z ≥ 1.3 or an upper limit for BL Lacs at lower redshift.Out of the 60 BL Lacs, we discovered 3 new high-z (z ≥1.3) objects and provided upper limits for 41.For the remaining 15 objects, our fits to the photometric data failed to provide an acceptable solution.This work, together with the previous three in this photometric campaign, have discovered a total of 19 BL Lacs with a z > 1.3 photometric redshift.Increasing the number of high-z BL Lacs is crucial since they are rare but essential in studying the blazar population and constraining the EBL.Our work shows that these objects have characteristics (high luminosity and high synchrotron peak frequencies) in between those of FSRQs and BL Lacs and that they could belong to the elusive class of 'masquerading BL Lacs'.Moreover, we find that the photons detected from these objects probe the high-z portion of the CGRH.We also find that the redshift uncertainty has negligible impact on the modeling of the EBL, allowing the use of these newly discovered high-z BL Lacs in future measurements of the EBL.Y.S. and M.A. acknowledge funding under NASA contracts 80NSSC21K1888 and 80NSSC21K1400.The authors acknowledge the Swift team for scheduling all the Swift-UVOT observations.They also acknowledge the Southeastern Association for Research in Astronomy for providing SARA-CT and SARA-RM observations.A.D. is thankful for the support of the Ramón y Cajal program from the Spanish MINECO.

Figure 2 .
Figure 2. The SEDs of the three high-redshift BL Lacs.The solid circles with black edges represent photometric magnitudes from Swift-UVOT, SARA-CT, and SARA-RM.The magnitudes are ordered in the following way: uw2, um2, uw1, uuu, ubb, g , uvv, r , i , z .The solid and dashed lines are power-law and galaxy models, respectively.

Figure 4 .
Figure4.The SED data points are from the 4FGL (green dots) and 3FHL (purple stars) catalogs.The green and purple areas are the power-law fits with 1 σ uncertainty provided in the 4FGL and 3FHL catalogs.The green and purple uncertainty bands come from the Fermi-LAT data reduction and include statistical uncertainties The black curves are the best joint power-law fits of 4FGL and 3FHL data attenuated by a factor e −τ(E,z) due to EBL absorption.The EBL model is taken fromDomínguez et al. (2011).The uncertainty region (the grey shaded area) shown for this joint fit is exclusively coming from the redshift uncertainty, according toDomínguez et al. (2011), and the gray area is not easily visible due to the small effects from the redshift uncertainty.The joint fit shows where the energy cutoff is.Note that the energy cutoff is solely due to the EBL attenuation.

Figure 5 .
Figure5.The colored dots are from the 4LAC catalog(Ajello et al. 2020).The triangles represent the sources from previous photo-z papers(Rau et al. 2012;Kaur et al. 2017;Kaur et al. 2018;Rajagopal et al. 2020), while the squares are found in this work.The vertical dotted lines are the divisions for LSP, ISP, and HSP blazars.The orange crosses are the masquerading BL Lacs fromPadovani et al. (2022) .Note that the indices are not corrected for EBL absorption.

Table 2 .
Swift-UVOT and SARA photometry (AB magnitudes corrected for extinction)