Multiwavelength Observations of the Blazar VER J0521+211 during an Elevated TeV Gamma-Ray State

The Very Energetic Radiation Imaging Telescope Array System (VERITAS) is an array of four imaging atmospheric Cerenkov telescopes (IACTs) located in southern Arizona (30°40'N, 110°57'W, 1.3 km a.s.l.; Park & The VERITAS Collaboration 2015). It is sensitive to gamma rays in the energy range from 85 GeV to >30 TeV with an energy resolution of ∼15% (at 1 TeV). The angular resolution (68% containment at 1 TeV) is ∼01. VERITAS is capable of making a detection at a statistical significance of 5 standard deviations (5σ) of a point source of 1% Crab Nebula flux in ∼25 hr, with an energy threshold of 240 GeV when a set of a priori data selection cuts optimized on sources with a moderate power-law index (from −2.5 to −3) is used.

After its initial discovery in the TeV band (Archambault et al. 2013), VER J0521+211 was observed again by VERITAS from 2012 November to 2014 February, with a total exposure of 23.6 hr after data quality selection and dead-time correction. We analyzed the data using two independent analysis packages (Cogan 2008; Maier & Holder 2017), with a priori cuts optimized for lower-energy showers (see, e.g., Archambault et al. 2014). The results of both analyses are compatible, and here we use the analysis package described in Cogan (2008).

MAGIC is an array of two IACTs located on the Canary Island of La Palma (28°45'43''N, 17° 53'24''W, 2.2 km a.s.l.; Aleksić et al. 2016). MAGIC is capable of detecting a point source of ∼0.7% of the Crab Nebula flux above 220 GeV in ∼50 hr with a statistical significance of 5σ (Aleksić et al. 2016). Above 220 GeV, its energy resolution is ∼16% and its angular resolution (68% containment) is ≲007.

VER J0521+211 was observed by MAGIC on four nights between 2013 October and 2013 December, triggered by elevated optical and GeV gamma-ray fluxes (Buson 2013). After data quality selection and dead-time correction, a total effective time of ∼4.5 hr was taken on 2013 October 15, October 16, November 29, and December 2. Details on the stereo analysis of the MAGIC data can be found in Aleksić et al. (2011).

The Large Area Telescope (LAT) on board the Fermi satellite is a pair-conversion gamma-ray telescope sensitive to energies from ∼20 MeV to >300 GeV (Atwood et al. 2009).

We performed an unbinned likelihood analysis covering the period of the TeV gamma-ray observations using the LAT ScienceTools v11r5p3 and Pass-8 P8R2_SOURCE_V6_v06 instrument response functions (Atwood et al. 2013). We selected SOURCE class events with an energy between 100 MeV and 300 GeV within 10° from the direction of VER J0521+211, with a maximum zenith-angle cut of 90°. We included all point sources within 20° from the source direction in our model, while only leaving the normalization parameters free for those within 5° from the source, and fixed the parameters to the Fermi Large Area Telescope Third Source Catalog (3FGL; Acero et al. 2015) values for those farther away. We also included the 3FGL Galactic (gll_iem_v06) and isotropic (iso_P8R2_SOURCE_V6_v06) diffuse emission in the model.

The X-Ray Telescope (XRT) on the Neil Gehrels Swift Observatory is a grazing-incidence focusing X-ray telescope, sensitive to energies from 0.2 to 10 keV (Gehrels et al. 2004; Burrows et al. 2005). There were 32 Swift-XRT observations taken on VER J0521+211 in the period of interest, 28 of which were taken in photon-counting (PC) mode between 2013 October and 2014 February, and four in window-timing (WT) mode in 2012 November. We processed these XRT data using xrtpipeline. ⁷⁷ We selected photons with an energy between 0.3 and 10 keV and within a circular source region of radius 20 pixels (∼ 472) to estimate the flux of the source, while those within an annular background region with inner and outer radii of 70 and 120 pixels ($\sim 2\buildrel{\,\prime}\over{.} 75-4\buildrel{\,\prime}\over{.} 72$) were selected to estimate the background contribution. For those observations with a count rate in the source region above 0.5 counts s⁻¹, we checked that the point-spread function is well described by the King function (Moretti et al. 2005), and therefore it is not necessary to correct for the pileup effect.

We used an absorbed log-parabola model to fit the X-ray spectrum:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={e}^{-{N}_{{\rm{H}}}\sigma (E)}K{\left(\displaystyle \frac{E}{1\,\mathrm{keV}}\right)}^{-\alpha -\beta \mathrm{log}(E/1\,\mathrm{keV})},\end{eqnarray} \tag{ 1 }$$

where in the absorption component N_H represents the column density of neutral hydrogen and σ(E) is the photoelectric cross section, while in the power-law component K is the normalization, α is the photon index at 1 keV, and β describes the spectral curvature, i.e., the energy dependence of the photon index. We fixed the value of N_H to the total neutral hydrogen density, N_H = 4.38 × 10²¹ cm⁻², derived from the Leiden/Argentine/Bonn (LAB) survey of Galactic Hi (Kalberla et al. 2005), following Willingale et al. (2013).

The Swift-Ultraviolet/Optical telescope (UVOT; Roming et al. 2005) was used in some of the observations of VER J0521+211 during the time period of interest. The counts from the source were extracted from an aperture with a radius of 5.0'' centered on the position of the source, and the background counts were estimated using the counts from four neighboring regions without any bright source, each having the same radius. The magnitude of the source was then computed using uvotsource, ⁷⁸ and it was later converted to flux using the zero-point for each of the UVOT filters (the central wavelengths of which are V 5468 Å, B 4392 Å, U 3465 Å, UVW1 2600 Å, UVM2 2246 Å, and UVW2 1928 Å) from Poole et al. (2008).

An extinction correction was applied following the procedure described in Roming et al. (2009), using the value of E(B − V) = 0.605 (Schlafly & Finkbeiner 2011). Since VER J0521+211 is located close to the Galactic plane, with a Galactic latitude of −87, the uncertainty on the dust extinction is large. A reddening value of E(B − V) = 0.703 was found for the direction toward VER J0521+211 using 472 elliptical galaxies to calibrate infrared dust maps (Schlegel et al. 1998). The value E(B − V) = 0.605 used in this work was derived more recently using 261,496 stars from the Sloan Digital Sky Survey (Schlafly & Finkbeiner 2011), which is 14% lower than the value from Schlegel et al. (1998). Schlafly & Finkbeiner (2011) pointed out that if substantial dust exists beyond a distance of 6 kpc, the reddening and extinction would be underestimated.

VER J0521+211 was monitored in the V band by the Steward Observatory (Smith et al. 2009). There were 76 observations on 37 nights between 2013 November 25 and 2014 April 1. In this work we use the V-band polarimetry results that are publicly available. ⁷⁹

3.1.1. Gamma-Ray Light Curves

From the VERITAS observations between 2012 November and 2013 February, which yielded a quality-selected live time of 4.6 hr, the source was detected at a statistical significance of 18σ at a time-averaged integral flux above 200 GeV of (2.4 ± 0.2) × 10⁻⁷ photons m⁻² s⁻¹. From the 16.4 hr of observations between 2013 October and 2014 February, the source was detected at 62σ and (3.9 ± 0.1) ×10⁻⁷ photons m⁻² s⁻¹ above 200 GeV.

An overall statistical significance of 30.5σ was obtained from the analysis of the MAGIC observations on the four nights during the period of interest. The energy threshold of these observations was estimated to be ∼65 GeV. The average flux above 200 GeV of these four MAGIC observations is (5.8 ± 1.0) × 10⁻⁷ photons m⁻² s⁻¹.

The TeV gamma-ray light curve (>200 GeV) of VER J0521+211 measured by VERITAS and MAGIC is shown in the top panel of Figure 1. The TeV gamma-ray flux of the source was variable between 2013 October and 2014 February (see the top panel of Figure 2), as a best-fit constant-flux model gives a χ² value of 442.4 with 19 degrees of freedom (dof), resulting in a negligible p-value. However, a fit of the TeV gamma-ray flux on the five nights with VERITAS observations between 2012 November and 2013 February to a constant-flux model yields a χ² value of 9.1 with 4 dof, equivalent to a p-value of 0.59. Therefore, we do not rule out constant flux from the source in late 2012 through early 2013.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The TeV gamma-ray flux above 200 GeV reached a peak of (8.8 ± 0.4) × 10⁻⁷ photons m⁻² s⁻¹ (∼37% of the Crab Nebula flux), as measured by VERITAS on the night of 2013 December 3. No intraday variability was found in the VERITAS data obtained on the night of the peak flux. On 2013 December 2, the night immediately before the peak TeV flux, MAGIC measured a lower flux above 200 GeV of (5.1 ± 0.9) × 10⁻⁷ photons m⁻² s⁻¹.

The monthly binned GeV gamma-ray light curve measured by Fermi-LAT between 2012 October and 2014 May is shown in the second panel of Figures 1 and 2. A fit of the GeV gamma-ray light curve to a constant-flux model yields a χ² value of 102.5 with dof = 22, equivalent to a negligible p-value, showing that the source was variable during the period of interest. The average flux value during the 1 yr period between 2013 January 29 and 2014 January 24 was about 11 times that from the 3FGL catalog (illustrated by the horizontal dashed lines in Figures 1 and 2), and the monthly flux stayed above 7 times the 3FGL catalog value, indicating a long-lasting elevated GeV gamma-ray flux state. During the period between 2013 October and 2014 February, when most of the observations at TeV gamma-ray and X-ray energies were performed, the GeV gamma-ray flux of the source did not show significant variability on monthly timescales.

3.1.2. Bayesian Blocks from the Gamma-Ray Light Curves

3.1.3. X-ray and UV/Optical Light Curves and Fractional Variability

The light curves measured by Swift-XRT, Swift-UVOT, and the Steward Observatory are shown in the lower panels of Figures 1 and 2. The flux of the source at all of the studied wavelengths, as well as the electric vector position angle (EVPA) and the polarization fraction in the optical V band, exhibited variability during the period of interest.

To examine the relative variability amplitude at different wavelengths, the fractional variability amplitude F_var (Vaughan et al. 2003; Poutanen et al. 2008) as a function of frequency is shown in Figure 3. The measured F_var value calculated from each light curve is shown in Table 2. The VHE flux exhibits the highest fractional variability amplitude of F_var = 0.54 ± 0.03 among all of the observed frequencies, with the X-ray flux yielding the second-highest F_var = 0.45 ± 0.02. The F_var values in the optical and UV bands are the lowest.

Zoom In Zoom Out Reset image size

Table 2. Fractional Variability Calculated from the MWL Light Curves Shown in Figure 1

Instrument	Fvar
VERITAS and MAGIC	0.54 ± 0.03
Fermi-LAT	0.36 ± 0.06
Swift-XRT (PC)	0.38 ± 0.02
Swift-XRT (WT)	0.34 ± 0.02
Swift-XRT (combined)	0.45 ± 0.02
Swift-UVOT (UVW2)	0.25 ± 0.03
Swift-UVOT (UVM2)	0.18 ± 0.05
Swift-UVOT (UVW1)	0.18 ± 0.02
Swift-UVOT (U)	0.16 ± 0.01
Swift-UVOT (B)	0.16 ± 0.01
Swift-UVOT (V)	0.14 ± 0.01
Steward Observatory (Pol. Frac.)	0.178 ± 0.001
Steward Observatory (EVPA)	0.244 ± 0.001

Download table as:  ASCIITypeset image

It is worth noting that F_var depends on the timescales, both long (duration) and short (bin width), of the light curve. A wider sampling coverage of the same light curve with a smaller bin width and a longer duration will generally increase F_var. The durations of the light curves used to calculate F_var are similar; so are the bin widths, with the exception of the monthly binned Fermi-LAT light curve.

3.1.4. Optical Polarization Variability

3.2.1. VHE Gamma-Ray Spectrum

The VERITAS spectra during the three periods BB1—BB3, as defined in Section 3.1.2, are shown in Figure 4. No evidence of spectral variability was observed between these flux states, despite the integrated flux being different by a factor of four. The best-fit spectral parameters for a power-law model and a log-parabola model (${N}_{0}{(E/300\,\mathrm{GeV})}^{[-\alpha -\beta {\mathrm{log}}_{10}(E/300\,\mathrm{GeV})]}$) for the three Bayesian blocks are shown in Table 3. The curved log-parabola model describes the spectra during BB1 and BB2 better than a power-law model. During the highest flux, in BB2, the VERITAS spectrum of the source extended beyond 1 TeV without any evidence of a sharp cutoff.

Zoom In Zoom Out Reset image size

Table 3. VHE Gamma-Ray Spectral Fit Parameters Using Power-law and Log-parabola Models

	Power Law		Log Parabola
	Γ	χ2/dof	α	β	χ2/dof
VERITAS BB1	3.2 ± 0.1	4.3	3.2 ± 0.1	1.0 ± 0.3	0.8
VERITAS BB2	3.1 ± 0.1	3.2	3.1 ± 0.1	0.9 ± 0.2	0.5
VERITAS BB3	3.4 ± 0.1	1.6	3.1 ± 0.1	0.5 ± 0.5	1.7
MAGIC October 15	3.0 ± 0.1	5.9	3.7 ± 0.2	1.9 ± 0.4	0.4
MAGIC October 16	2.9 ± 0.2	0.8	3.1 ± 0.3	0.8 ± 0.7	0.6
MAGIC November 29	3.1 ± 0.1	8.5	3.7 ± 0.2	1.7 ± 0.3	1.1

Download table as:  ASCIITypeset image

The MAGIC spectra measured on three of the four nights, 2013 October 15, October 16, and November 29, are shown in Figure 5, together with a log-parabola fit. The best-fit spectral parameters and reduced χ² values for these three nights are shown in Table 3. Note that the MAGIC spectra extend down to a lower energy (70 GeV) than those from VERITAS, thanks to the larger areas of the reflectors, capturing slightly more spectral curvature at low energies.

Zoom In Zoom Out Reset image size

3.2.2. X-Ray and UV/Optical Spectral Variability

The X-ray spectrum from each Swift-XRT observation was extracted. The detailed information of these observations and the results of the spectral fits are shown in Table 5 in the Appendix. The relation between the best-fit photon index at 2 keV ($\alpha +\beta \mathrm{log}2$) and the energy flux in XRT data is illustrated by Figure 6. The log-parabola energy-dependent photon index at 2 keV is comparable to the power-law index for this source and therefore chosen to investigate the X-ray index−flux correlation. A "harder-when-brighter" trend, i.e., a smaller photon index when the flux is higher, is apparent. The Pearson correlation coefficient between the index and flux is 0.76 (p-value ∼6 × 10⁻⁷). A linear fit yields a slope of −0.2 ± 0.04.

Zoom In Zoom Out Reset image size

A combined analysis was performed using observations taken during each of the three Bayesian blocks. The best-fit power-law model suggests that the X-ray spectrum during the last Bayesian block was marginally softer (with an index at 2 keV of 2.73 ± 0.11) than that for the first Bayesian block (with an index at 2 keV of 2.57 ± 0.06).

The UV/optical spectrum of the source can be obtained by combining the extinction-corrected fluxes at the six wavelengths from the Swift-UVOT filtered observations. However, the uncertainty in the measured reddening affects the extinction correction, most notably in the UV band, and subsequently affects the interpretation of the frequency of the synchrotron peak in the broadband SED. If the true reddening is closer to the value E(B − V) = 0.703 from Schlegel et al. (1998), the frequency of the synchrotron peak of the SED would be above the UV band, making this source an HBL during the flaring state. However, if the true reddening is the value E(B − V) = 0.605 from Schlafly & Finkbeiner (2011), as we used in this work, the frequency of the synchrotron peak would be below 10¹⁵ Hz, still in the IBL regime.

Because of the uncertainty in the measured reddening, we did not include data taken with the two Swift-UVOT filters with the highest frequencies in the UV/optical spectra. The UVOT spectra at the lower four frequencies during the three periods were fitted with a power-law model to estimate the spectral variability. The UV spectrum is marginally harder during BB2, with a best-fit index of 2.23 ± 0.09, whereas those during BB1 and BB3 are 2.41 ± 0.06 and 2.37 ± 0.05, respectively. This is consistent with the "harder-when-brighter" behavior in the X-ray band. Because the extinction correction is the same for all periods, the relative change is robust against the chosen reddening value.

3.2.3. X-Ray/Gamma-Ray Correlation

A moderate correlation between TeV gamma-ray and X-ray fluxes observed within a 1-day window is apparent, as shown in Figure 7, with a Pearson correlation coefficient r between the two bands of 0.74 (p-value 0.0014). Note that the small number of data points and the potential non-Gaussian distribution of the observed fluxes of the source make it hard to interpret the Pearson's r correlation. A linear fit, a quadratic fit, and a power-law fit (${F}_{\mathrm{TeV}}\propto {F}_{{\rm{X}}}^{a}$) were performed. All three models yielded rather large reduced χ² values (between 8.4 and 9.6), likely due to the intrinsic scatter of the gamma-ray/X-ray correlation. The best-fit power-law index is a = 1.55 ± 0.12, indicating that the correlation between TeV gamma-ray and X-ray fluxes is likely between quadratic (a = 2) and linear (a = 1).

Zoom In Zoom Out Reset image size

3.2.4. Optical–UV/Gamma-Ray Correlation

The correlations between the TeV gamma-ray flux and the optical/UV energy fluxes measured with the six Swift-UVOT filters are shown in Figure 8. The Pearson correlation coefficients range between 0.47 and 0.76 (the p-values range between 0.01 and 0.08 except for the UVW2 filter), reaching the highest value of 0.76 (a p-value of 0.0015) between the UV flux observed with the highest-frequency UVW2 filter and the TeV gamma rays (as shown in the bottom right panel of Figure 8). A power-law fit (${F}_{\mathrm{TeV}}\propto {F}_{\mathrm{UVOT}}^{b}$) indicates a more-than-quadratic relation with indices b > 2 between these bands, although the goodness of fit is poor with reduced χ² values ranging between 6.4 (between UVW2 and VHE) and 16.7 (between UVM2 and VHE). However, the fractional variability amplitudes F_var in the optical/UV bands are much lower than in the VHE and X-ray bands, increasing with frequency from 0.14 ± 0.01 for the lowest-frequency V filter to 0.25 ± 0.03 for the highest-frequency UVW2 filter, as shown in Figure 3. The low F_var values for optical/UV fluxes would generally lead to a steeper optical/UV/VHE correlation with a large index b. We note that the F_var values are not affected by extinction correction.

Zoom In Zoom Out Reset image size

The broadband SEDs were constructed for the three Bayesian block intervals, BB1, BB2, and BB3, and are shown in Figure 9. The EBL absorption of the VHE gamma-ray spectra was considered following Franceschini et al. (2008), assuming a redshift of z = 0.18 (Paiano et al. 2017). The X-ray spectra were corrected for the neutral hydrogen absorption (see Section 2.4). The UV/optical spectra were corrected for extinction (see Section 2.5). The SEDs are modeled with a one-zone SSC model as described below in Section 4.4.

Zoom In Zoom Out Reset image size

The blazar VER J0521+211 was observed at an elevated gamma-ray state between 2013 and 2014 with VERITAS, MAGIC, and Fermi-LAT. The TeV gamma-ray flux above 200 GeV reached a peak of (8.8 ± 0.4) × 10⁻⁷ photons m⁻² s⁻¹ (∼37% of the Crab Nebula flux) on the night of 2013 December 3. The MAGIC observations on 2013 December 2, 1 day before the night of the peak TeV gamma-ray flux, yielded a flux above 200 GeV of (5.1 ± 0.9) × 10⁻⁷ photons m⁻² s⁻¹, implying variability on a daily timescale. The GeV gamma-ray flux above 100 MeV stayed above 5 times the value from 3FGL for over 5 months between 2013 September and 2014 February, with an average flux of 10 times the 3FGL value.

MWL variability from the optical to VHE bands was observed in 2013 and 2014, with the fractional variability amplitude F_var increasing with frequency from the lowest 0.14 ± 0.01 for V-band flux to 0.45 ± 0.02 for X-ray flux, and then decreasing to 0.36 ± 0.06 for GeV gamma-ray flux, before increasing again to the highest value 0.54 ± 0.03 for TeV gamma-ray flux. Such a pattern in F_var is consistent with previous observations of gamma-ray blazars (see, e.g., Baloković et al. 2016, and references therein), with the highest F_var observed at the frequencies above the two peaks in the SED.

A moderate correlation between the TeV gamma-ray and X-ray fluxes from the source was observed. The best-fit power-law model suggests that the correlation is more than linear, with an index of 1.47 ± 0.13, consistent with the expectation of a varying electron density and adiabatic cooling in an SSC model (Krawczynski et al. 2002; Katarzyński et al. 2005). The TeV gamma-ray and X-ray correlation could also give insights into whether the IC scattering occurs in the Thomson regime or KN regime, as the correlation may become closer to linear in the KN regime with the diminished scattering cross section (e.g., Aleksić et al. 2015; Baloković et al. 2016). More data are needed to quantify the correlation between X-rays and TeV gamma rays with a higher precision.

The correlation between the optical and TeV gamma-ray fluxes was weaker than that between X-ray and TeV gamma-ray fluxes, while the correlation between UV and TeV gamma-ray fluxes was comparable to that between the X-rays and TeV gamma rays, with a slightly higher correlation coefficient observed between the UVOT measurements with the UVW2 filter at the highest central frequency of ∼1.55 × 10¹⁵ Hz. The UV/TeV gamma-ray correlation is more than quadratic, which could occur if the emitting region expands as the particle density increases (Katarzyński et al. 2005).

Despite multiple measurement attempts (Shaw et al. 2013; Paiano et al. 2017), the redshift of VER J0521+211 remains uncertain. To facilitate the discussion, we adopt the value z = 0.18, the lower limit reported by Paiano et al. (2017). At this redshift, the luminosity distance is 878 Mpc (using Ω_m = 0.29, Ω_Λ = 0.71, and H₀ = 69.6 km s⁻¹ Mpc⁻¹; Bennett et al. 2014), and the angular scale is roughly 3 pc mas⁻¹.

To help constrain the redshift of the source, upper limits were derived from the gamma-ray spectra using two methods based on Mazin & Goebel (2007) and Georganopoulos et al. (2010).

The first method (Mazin & Goebel 2007) assumes that the intrinsic gamma-ray spectrum follows a power law with an index of Γ_int > 1.5 (or more conservatively, Γ_int > 0.7). This method was used to derive the 95% upper limit on the redshift of z < 0.34 (and z < 0.44 for the conservative assumption of Γ_int > 0.7) for VER J0521+211 from VERITAS observations in 2009 and 2010 (Archambault et al. 2013). Using the same method, we find that the VERITAS spectrum on the night of the flare (BB2) offers a stronger constraint than other time periods, yielding a 95% upper limit on the redshift of z < 0.38 (and z < 0.53 for the conservative assumption of Γ_int > 0.7). We note that the deabsorbed VERITAS spectra at z = 0.38 are convex, with the measured fluxes at high energies significantly exceeding the best-fit power-law model, indicating that the redshift of the source is likely to be well below these limits.

The second method (Georganopoulos et al. 2010) assumes that the intrinsic TeV gamma-ray flux does not exceed the extrapolation of the GeV gamma-ray spectrum, which is minimally affected by the EBL for TeV blazars. A 95% upper limit on the optical depth from the EBL absorption can be calculated following Equation (2) in Aleksić et al. (2011):

$$\begin{eqnarray*}&&{\tau }_{\gamma \gamma }(E)\leqslant \mathrm{log}\left[\displaystyle \frac{{F}_{\mathrm{int}}(E)}{{F}_{\mathrm{obs}}-1.64{\rm{\Delta }}{F}_{\mathrm{obs}}}\right],\end{eqnarray*}$$

where F_int(E) is the extrapolated GeV flux, representing the maximally allowed intrinsic flux, at energy E; F_obs is the observed TeV gamma-ray flux; and ΔF_obs is the uncertainty of F_obs. The 95% upper limits on the optical depth τ_{γ γ} (E) derived from Fermi-LAT and VERITAS spectra for periods BB1, BB2, and BB3 are shown in the top panel of Figure 10. The spectrum from BB2 offers the strongest constraints on τ_{γ γ} (E). The optical depths τ_{γ γ} (z, E) for sources at a range of redshift of 0.25 ≤ z ≤ 0.35 from four EBL models (Franceschini et al. 2008; Finke et al. 2010; Domínguez et al. 2011; Gilmore et al. 2012) were computed and compared to the constraints from VERITAS observations. The minimum difference between the observational constraints and the model-predicted optical depth, ${\rm{\Delta }}{\tau }_{\gamma \gamma ,\min }(z,E)$, as a function of redshift is shown in the bottom panel of Figure 10. The 95% upper limits on the redshift of the source (corresponding to ${\rm{\Delta }}{\tau }_{\gamma \gamma ,\min }(z,E)=0$) were found to be z ≤ 0.3, z ≤ 0.278, z ≤ 0.308, and z ≤ 0.281 for the models from Domínguez et al. (2011), Finke et al. (2010), Franceschini et al. (2008), and Gilmore et al. (2012), respectively. These upper limits on the redshift are more constraining than the first method and those reported by Archambault et al. (2013). A precise measurement of the spectroscopic redshift of this source in the future will be important, as it helps constrain both the particle distribution and the EBL intensity.

Zoom In Zoom Out Reset image size

Three periods were identified with a Bayesian block analysis during which the source exhibited different TeV gamma-ray flux states (BB1, BB2, and BB3 as mentioned above). Broadband SEDs from UV to TeV gamma rays were constructed for each period. These SEDs are adequately described by a one-zone SSC model calculated following Krawczynski et al. (2002). The sets of parameters used to generate the SEDs shown in Figure 9 are listed in Table 4. These parameters differ from those in the SSC model for the low-state SED from VERITAS data in 2009 (Archambault et al. 2013), where an even lower magnetic field strength of B = 0.0025 G results in a long synchrotron cooling time and a particle-dominated jet away from equipartition.

Table 4. Parameters Used for the One-zone SSC Models Shown in Figure 9

State	Γ	θ	B	R	we	$\mathrm{log}\tfrac{{E}_{\min }}{\mathrm{eV}}$	$\mathrm{log}\tfrac{{E}_{\max }}{\mathrm{eV}}$	$\mathrm{log}\tfrac{{E}_{\mathrm{break}}}{\mathrm{eV}}$	p1	p2
		(deg)	(10−2 G)	(1017 cm)	(10−3 erg cm−3)
BB1	25	2.2	1.5	1.01	1.32	9.7	12.0	11.25	3.2	4.2
BB2	25	2.2	1.5	1.05	1.45	9.7	12.2	11.25	3.15	4.15
BB3	25	2.2	1.5	1.1	1.0	9.7	11.9	11.15	3.25	4.25
2009 a	30	⋯	0.25	4.0	⋯	10.25	12.0	⋯	3.0	⋯

Note. Γ and θ are the bulk Lorentz factor and the angle with respect to the line of sight of the relativistic emitting region, respectively; R is the radius of this region; B is the strength of the magnetic field; w_e is the energy density of the emitting electrons; $\mathrm{log}{E}_{\min }$, $\mathrm{log}{E}_{\max }$, and $\mathrm{log}{E}_{\mathrm{break}}$ are the logarithms of the minimum, maximum, and break electron energies, respectively; p₁ and p₂ are the spectral indices of the electrons below and above the break energy, respectively. The corresponding parameters listed in Table 3 in Archambault et al. (2013) for the low-state SED in 2009 are also included for a comparison.

^a Archambault et al. (2013).

Download table as:  ASCIITypeset image

The simple SSC modeling provides several constraints on the properties of the source that are discussed below.

The peak frequency of the synchrotron radiation from the SSC model is roughly 5 × 10¹⁴ Hz, in the optical band and the IBL regime. However, as mentioned in Section 2.5, the synchrotron peak frequency closely depends on the intrinsic optical/UV spectrum and is heavily affected by extinction correction. If more substantial extinction from dust exists, the synchrotron peak would be in the UV band and the source in the HBL regime.

We note that it has been observed that the synchrotron peak frequency can shift between different flux states for blazars (e.g., Acciari et al. 2011; Valverde et al. 2020). Particularly, VER J0521+211 has been reported to exhibit HBL-like properties during a previous flare on 2009 November 27 (Archambault et al. 2013). The "flare" and "low-state" SEDs from Archambault et al. (2013) are shown, as a comparison to this work, as the magenta downward-pointing triangles and gray upward-pointing triangles in Figure 9. The TeV gamma-ray spectrum, including normalization and index, during the entire period studied in this work was comparable to that during the "flaring" state in 2009 November. However, while the X-ray spectra in this work were similar in normalization, they were much softer than those during the "flaring" state in 2009 November (with a photon index value of 2.0 ± 0.1; Archambault et al. 2013). This suggests strong synchrotron X-ray spectral variability on timescales of years for VER J0521+211. Moreover, the optical fluxes in this work were comparable to the low-state flux in Archambault et al. (2013), with a rather low variability amplitude. As the optical flux is only mildly affected by the extinction correction, this low variability amplitude in the optical band is robust.

The low optical variability implies that the SED at/below the synchrotron peak frequency does not change as much as that above the peak. The relatively flat spectrum from UV to X-ray bands indicates that the synchrotron peak is likely broad, making it difficult to measure the peak frequency and classify the blazar. On the other hand, the UV/X-ray spectrum of the source becomes harder with higher flux (as illustrated by the UV photon indices in Section 2.5 and the index-flux relation in Figure 6). The similar behavior of the observed UV and X-ray fluxes implies that the optical/UV bands are also likely located above the synchrotron SED peak, confirming the source as an IBL even in the higher flux state in 2013. The "harder-when-brighter" trend in both bands could be caused by the synchrotron SED peak shifting toward higher frequencies during flares. Such a behavior of the SED can be explained by changes to particles with the highest energies, e.g., from an evolving maximum or break particle energy or a change in the particle distribution shape above the break.

As the synchrotron peak frequency of VER J0521+211 is around the borderline between IBL and HBL, we review the radio features in the innermost jet region, which also provide insights into the blazar subclasses. Lister et al. (2019) identified six moving features in the inner jet from 15 GHz VLBA observations since 2009, including one that moves at a large apparent speed of ∼0.77 ± 0.05 mas yr⁻¹ at ∼5 mas from the core in mid-2014. On the other hand, three innermost features move at a very small apparent velocity, two of which even exhibited statistically significant inward motion. These slow features may be quasi-stationary. Such a combination of quasi-stationary and fast-moving radio knots is typically found in IBL/LBLs (e.g., Hervet et al. 2016). Therefore, the MOJAVE observations are consistent with the source being an IBL.

Interestingly, a 15 GHz moving feature at an apparent speed of ∼ 0.21 ± 0.02 mas yr⁻¹ first appeared at around 5.5 mas (∼11 pc) in 2012 December (Lister et al. 2018), ⁸⁰ coincident with the period when the source exhibited elevated gamma-ray flux. Similar coincidences have been observed before in IBLs, but between fast gamma-ray flares on intraday timescales and moving knots at distances much closer to the core observed with the VLBA at 43 GHz (e.g., from BL Lacertae; Arlen et al. 2013; Abeysekara et al. 2018). The interactions between the moving and stationary radio features in the jet could provide a possible particle acceleration mechanism that leads to enhanced gamma-ray emission.

Based on the SSC model parameters (see Table 4), we can investigate the emission process of the particles at the source. The Lorentz factor of those electrons emitting the peak synchrotron radiation γ_syn ≈ 5.4 × 10⁴(B/1 G)^−1/2 ${[\delta /{(1+z)]}^{-1/2}({\nu }_{\mathrm{syn}}/{10}^{16}\,\mathrm{Hz})}^{1/2}\approx 2\times {10}^{4}$ (for BB1/BB3), where B is the strength of the magnetic field, δ is the Doppler factor ($\delta ={[{\rm{\Gamma }}(1-\beta \cos \theta )]}^{-1}=26$, where $\beta ={(1-1/{{\rm{\Gamma }}}^{2})}^{1/2}$), z is the redshift, and ν_syn is the peak frequency of the synchrotron radiation. The synchrotron cooling time of these electrons emitting at the peak frequency is t_syn ≈ 7.74 × 10⁸ s × (B/1 G)⁻² × γ⁻¹ ≈ 2000 days, corresponding to ∼78 days in the observer's frame. The electrons responsible for the peak synchrotron emission during BB2 have a slightly higher Lorentz factor of 2.1 × 10⁴ and a cooling time of 73 days in the observer's frame. As the synchrotron radiation energy density in the model is a few times larger than the magnetic field energy density, the cooling is likely dominated by the IC scattering, and the total observed cooling time can be as short as 2 weeks. Since the dynamic timescale is R/c ∼ 1.5 days in the observers' frame, the adiabatic cooling could also dominate the decay of the flux of the source.

Using the Lorentz factor of the electrons emitting at the synchrotron peak frequency γ_syn, we checked for the importance of the KN effects, which become substantial when the photon energy h ν^⋆ > m_e c² in the electron rest frame. Considering the IC scattering between the photons at the synchrotron peak frequency in the comoving frame $E{{\prime} }_{\mathrm{peak}}\approx 0.1\,\mathrm{eV}$ (primed quantities refer to those in the comoving frame) and the electrons that produce the peak synchrotron photons, we have ${\gamma }_{\mathrm{syn}}E{{\prime} }_{\mathrm{peak}}\approx 2\,\,\mathrm{keV}$ < m_e c². Therefore, the IC scattering of the peak of the synchrotron radiation field by electrons with γ_syn occurs in the Thomson regime (Acciari et al. 2011) and produces the observed IC radiation up to $4\delta {\gamma }_{\mathrm{syn}}^{2}E{{\prime} }_{\mathrm{peak}}\approx 4\,\,\mathrm{GeV}$, below the high-energy SED peak.

From the SED modeling, the SSC component peaks at a frequency ν_ssc ≈ 4.4 × 10²⁴ Hz for the first two blocks and moves to a lower frequency ν_ssc ≈ 2.6 × 10²⁴ Hz for the last, low-flux block. The shift of the SSC peak to higher frequencies during flares has been previously observed in blazars (e.g., Acciari et al. 2011). We also calculate the limit of the KN Doppler factor (Equation (17) in Tavecchio et al. 1998) δ < δ_KN $={[4{\nu }_{\mathrm{syn}}{\nu }_{\mathrm{ssc}}/(3{({m}_{e}{c}^{2}/h)}^{2})]}^{1/2}$ $\exp \{1/({\alpha }_{1}-1)+1/[2({\alpha }_{2}-{\alpha }_{1})]\}{}^{-1/2}$, where the spectral indices below and above the IC peak are taken from Fermi-LAT and VERITAS observations as α₁ ≈ 0.9 and α₂ ≈ 2.1, respectively. In this case, the Doppler factor δ = 26 is below the KN limit δ_KN ranging between about 40 and 60. As a result, the IC scattering that produces the peak SSC emission occurs in the KN regime. Note, however, that the observed SSC peak lies around 15 GeV; therefore, the best-fit index from the Fermi-LAT spectrum may not be a good estimator of α₁. Moreover, the KN limit Doppler factor, defined as Equation (17) in Tavecchio et al. (1998), depends sensitively on α₁. If we use a value of α₁ ≈ 0.8 as an example, the KN limit Doppler factor becomes about 6, placing the gamma-ray SED peak in the Thomson regime.

The maximum electron energy ${E}_{e,\max }$ ranges between 10^11.9 and 10^12.2 eV (γ between about 1.6 × 10⁶ and 3.1 × 10⁶), and the observed maximum gamma-ray energy is E_γ ∼ 1 TeV, corresponding to an energy of ${E}_{\gamma }^{\prime} ={E}_{\gamma }(1+z)/\delta \sim 40\,\,\mathrm{GeV}$ in the comoving frame. Therefore, we can estimate the energy range of the low-energy photons that participate in IC scattering and produce the observed gamma rays to be roughly between ${E}_{\gamma }^{\prime} {m}_{e}^{2}{c}^{4}/{E}_{e,\max }^{2}\approx 0.01\,\mathrm{eV}$ and ${m}_{e}^{2}{c}^{4}/{E}_{\gamma }^{\prime} \sim 6.8\,\mathrm{eV}$ in the comoving frame. These correspond to the observed photon energies between 0.25 and 175 eV, or between approximately 6 × 10¹³ Hz and 4 × 10¹⁶ Hz, ranging from the IR/optical bands up to the UV and soft X-ray band. In this scenario, we expect the observed emission at some of these lower energies to be strongly correlated with the gamma-ray emission. If the target photons of the observed energies above $\delta {m}_{e}{c}^{2}/{\gamma }_{\max }\sim$ 4 eV were upscattered by the electrons at the maximum energy, the IC scattering would occur in the KN regime, leading to spectral breaks in the synchrotron spectrum that are not observed. However, the observed correlations between optical/UV and VHE fluxes (as shown in Figure 8) are not stronger than the observed X-ray/VHE correlation (as shown in Figure 7), implying that the lower-energy electrons and the target photons with observed energies above ∼6.4 eV (corresponding to the frequency of the UVW2 filter) are probably responsible for the observed VHE flux.

The electron break energy E_brk ranges between 10^11.15 and 10^11.25 eV (γ between 2.8 × 10⁵ and 3.5 × 10⁵). These electrons emit synchrotron radiation at around 2 × 10¹⁷ Hz (roughly 1 keV) and have a synchrotron cooling time of a few days in the observers' frame. Assuming that the mean free path (λ) of these electrons is comparable to their gyroradius of between 3.1 × 10¹⁰ cm and 3.9 × 10¹⁰ cm, we can estimate the electron diffusion coefficient D ≈ λ²/(2λ/c) to be a few times 10²⁰ cm² s⁻¹. Taking the radius of the emitting region used in the SED modeling, we can estimate the escape timescale of the electrons at the break energy τ_esc ≈ R²/(4D) to be on the order of 10⁵ yr in the comoving frame, much longer than the cooling time.

The light-crossing time in the observer's frame is $2R^{\prime} (1+z)/(c\delta )\sim 3$ days, shorter than the shortest Bayesian block interval (∼5 days), and the model obeys causality.

A leptonic model has been used to describe the observed SED from VER J0521+211, allowing information about the jet environment to be inferred. It is possible that protons are accelerated in this same jet environment. By requiring protons to be magnetically confined within the emitting region (Hillas 1984), a rough estimation of the maximum proton energy can be calculated as E_p < 3 × 10⁻¹¹ PeV(B/1 G)(R/1 m) ≈ 500 PeV. If photohadronic processes occur for protons at the maximum energy, the corresponding maximum energy of the produced neutrinos is then E_ν ≤ γ_p m_π c²/4 ≈ 0.04E_p ≈20 PeV (Rachen & Mészáros 1998), where γ_p is the Lorentz factor of the protons and m_π is the mass of a charged pion. Since no considerations were made of the acceleration, radiation, or escape processes of the protons, the maximum proton energy could be lower owing to the constraints from these processes. However, we do not rule out the possibility of neutrino production at energies up to 20 PeV in VER J0521+211 from the confinement argument alone. Given the relatively long duration of the elevated TeV gamma-ray flux, VER J0521+211 is potentially a good candidate for astrophysical neutrino searches.