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

We report on a long-lasting, elevated gamma-ray flux state from VER J0521+211 observed by VERITAS, MAGIC, and Fermi-LAT in 2013 and 2014. The peak integral flux above 200 GeV measured with the nightly-binned light curve is $(8.8 \pm 0.4) \times 10^{-7} \;\text{ph}\;\text{m}^{-2}\; \text{s}^{-1}$, or ~37% of the Crab Nebula flux. Multiwavelength observations from X-ray, UV, and optical instruments are also presented. A moderate correlation between the X-ray and TeV gamma-ray fluxes was observed, and the X-ray spectrum appeared harder when the flux was higher. Using the gamma-ray spectrum and four models of the extragalactic background light (EBL), a conservative 95% confidence upper limit on the redshift of the source was found to be z<=0.31. Unlike the gamma-ray and X-ray bands, the optical flux did not increase significantly during the studied period compared to the archival low-state flux. The spectral variability from optical to X-ray bands suggests that the synchrotron peak of the spectral energy distribution (SED) may become broader during flaring states, which can be adequately described with a one-zone synchrotron self-Compton model varying the high-energy end of the underlying particle spectrum. The synchrotron peak frequency of the SED, as well as the radio morphology of the jet from the MOJAVE program, are consistent with the source being an intermediate-frequency-peaked BL Lac object.

Interesting differences in radio morphology between subclasses of blazars, beyond the wellknown "blazar sequence" (Fossati et al. 1998), have been recently proposed (e.g., Kharb et al. 2008;Hervet et al. 2016).Particularly, some blazars, typically belonging to the IBL/LBL subclass, exhibit a combination of stationary knots and superluminal knots in their jets as observed with the Very Long Baseline Array (VLBA).Gamma-ray flares contemporaneous with the ejection of superluminal knots have been observed in a few TeV IBLs (e.g., Arlen et al. 2013;Abeysekara et al. 2018;MAGIC Collaboration et al. 2018), providing tantalizing evidence for recollimation shocks (e.g., Komissarov & Falle 1998;Narayan et al. 2009).
The discovery of VER J0521+211 in the TeV gamma-ray band was made from VERI-TAS observations of the radio and X-ray source RGB J0521.8+2112,motivated by a cluster of photons with energies >30 GeV in Fermi-LAT data (Archambault et al. 2013).The positions of the source determined with X-ray, GeV and TeV gamma-ray data are all within 0.1 • .The spatial association across the electromagnetic spectrum and correlated variability, most prominent between X-ray and gamma-ray bands, strongly support the association between VER J0521+211 and RGB J0521.8+2112.
As is the case with many other BL Lac objects, the lack of optical emission features leads to a redshift uncertainty of the source.A redshift of z = 0.108 was reported by Shaw et al. (2013), but later studies were unable to confirm it (Archambault et al. 2013;Paiano et al. 2017).
Instead, a lower limit of z > 0.18 (Paiano et al. 2017) and an upper limit of z < 0.34 (Archambault et al. 2013) were reported.
VER J0521+211 is classified as an IBL as the peak of its spectral energy distribution in a low state lies between 10 14 Hz and 10 15 Hz.However, the source exhibits a harder X-ray spectrum and HBL-like behavior during flares.(Archambault et al. 2013).In HBLs, it is particularly interesting to study the flux correlation between X-ray and TeV gamma-ray bands, where the SED peaks lie.Despite being model-dependent, the quantitative description of the correlation between X-rays and TeV gamma rays can be informative of the radiative processes of the relativistic particles (e.g.Katarzyński et al. 2005).For example, a particle injection in a one-zone synchrotron self-Compton (SSC) model in the Thomson regime would lead to a quadratic correlation between the gamma-ray and X-ray flux, but could lead to a linear correlation if the inverse-Compton (IC) scattering happens in the deep Klein-Nishina (KN) regime (e.g.Aleksić et al. 2015;Baloković et al. 2016).
In this work, we report on the results from multiwavelength (MWL) observations of the blazar VER J0521+211 between 2012 Nov and 2014 Feb, focusing on the TeV gamma-ray and X-ray behaviors.

VERITAS
The Very Energetic Radiation Imaging Telescope Array System (VERITAS) is an array of four imaging atmospheric Cherenkov 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 ∼0.1 • .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 hours, 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 Nov to 2014 Feb, with a total exposure of 23.6 h 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 analysis are compatible and here we use the analysis package described in Cogan (2008).

MAGIC
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. 2016a).MAGIC is capable of detecting a point source of ∼ 0.7% of the Crab Nebula flux above 220 GeV in ∼ 50 hours with a statistical significance of 5 σ (Aleksić et al. 2016b).Above 220 GeV, its energy resolution is ∼ 16% and its angular resolution (68% containment) is 0.07 • .VER J0521+211 was observed by MAGIC on four nights between 2013 Oct and 2013 Dec, triggered by elevated optical and GeV gammaray fluxes (Buson 2013).After data quality selection and dead time correction, a total effective time of ∼ 4.5 hours was taken on 2013 Oct 15, Oct 16, Nov 29, and Dec 2. Details on the stereo analysis of the MAGIC data can be found in Aleksić et al. (2011).

Fermi-LAT
The Large Area Telescope (LAT) onboard 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 gammaray 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 further away.We also included the 3FGL Galactic (gll_iem_v06) and isotropic (iso_P8R2_SOURCE_V6_v06) diffuse emission in the model.

Swift-XRT
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 keV 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 Oct and 2014 Feb, and four in window-timing (WT) mode in 2012 Nov. We processed these XRT data using xrtpipeline. 1We selected photons with an energy between 0.3 keV and 10 keV and within a circular source region of radius 20 pixels (∼ 47.2 ) to estimate the flux of the source; while those within an annular background region with inner and outer radii of 70 and 120 pixels (∼ 2.75 −4.72 ) were selected to estimate the background contribution.For those obser-vations with a count rate in the source region above 0.5 cts s −1 , 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 pile-up effect.
We used an absorbed log-parabola model to fit the X-ray spectrum: 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 21 cm −2 , 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 of 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, 2 and later converted to flux using the zero-point for each of the UVOT filters (the central wavelengths of which are: 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 -8.7 • , the uncertainty on the dust extinction is large.A reddening value of E(B − V )=0.703 was found for the direction towards 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 261496 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.

Steward Observatory
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 Nov 25 and 2014 Apr 1.We use the V -band polarimetry results from that are publicly available3 in this work.

Gamma-Ray Light Curves
From the VERITAS observations between 2012 Nov and 2013 Feb, which yielded a qualityselected live time of 4.6 h, 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 −7 photon m −2 s −1 .From the 16.4 hours of observations between 2013 Oct and 2014 Feb, the source was detected at 62σ and (3.9 ± 0.1) × 10 −7 photon m −2 s −1 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 −7 photon m −2 s −1 .
The TeV gamma-ray light curve (> 200GeV) 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 Oct and 2014 Feb (see the top panel of Figure 2), as a best-fit constant-flux model gives a χ 2 value of 442.4 with 19 degrees of freedom (DOF), resulting in a negligible pvalue.However, a fit of the TeV gamma-ray flux on the five nights with VERITAS observations between 2012 Nov and 2013 Feb to a constantflux model yields a χ 2 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.
The TeV gamma-ray flux above 200 GeV reached a peak of (8.8 ± 0.4) × 10 −7 photon m −2 s −1 (∼37% of the Crab Nebula flux), as measured by VERITAS on the night of 2013 Dec 3. No intra-day variability was found in the VERITAS data obtained on the night of the peak flux.On 2013 Dec 2, the night immediately before the peak TeV flux, MAGIC measured a lower flux above 200 GeV of (5.1 ± 0.9) × 10 −7 photon m −2 s −1 .
The monthly-binned GeV gamma-ray light curve measured by Fermi-LAT between 2012 Oct and 2014 May is shown in the second panel from top of Figures 1 and 2. A fit of the GeV gamma-ray light curve to a constant-flux model yields a χ 2 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 one-year period between 2013 Jan 29 and 2014 Jan 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 seven times the 3FGL catalog value, indicating a long-lasting elevated GeV gamma-ray flux state.During the period between 2013 Oct and 2014 Feb, 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.

Bayesian Blocks from the Gamma-Ray Light Curves
A Bayesian blocks analysis (Scargle et al. 2013) of the TeV gamma-ray light curve was performed using the Python package Astropy (Astropy Collaboration et al. 2013, 2018) with an empirical prior equivalent to a false alarm probability p 0 = 3 × 10 −7 .Two change points (i.e., the edges of the Bayesian blocks), 2013 Dec 2 and 2013 Dec 6, were obtained from the VERITAS and MAGIC flux measurements, as illustrated by the light blue vertical dashed lines in Figure 1 and Figure 2.These two change points define three Bayesian blocks, designated as BB1, BB2, and BB3, within each of which the flux is consistent with being constant and no variability is detected above 5-σ statistical significance.The date ranges and the average fluxes of the three Bayesian blocks are shown in Table 1.Note that there was only one night with 2.3 h of VERITAS observations on 2013 Dec 3 in the BB2 interval, with an averaged flux above 200 GeV of (8.8 ± 0.4) × 10 −7 photon m −2 s −1 (∼37% of the Crab Nebula flux), the highest of all three Bayesian blocks.The flux in BB3 is comparable to that in 2012, (2.4 ± 0.2) × 10 −7 photon m −2 s −1 .MWL SEDs were constructed for the three time intervals from the Bayesian blocks analysis of the TeV gamma-ray light curve (see Section 3.2).
The same Bayesian block analysis was performed for the monthly binned Fermi -LAT light  2013-09-25 2013-10-15 2013-11-04 2013-11-24 2013-12-14 2014-01-03 2014-01-23 2014-02-12 Figure 2. MWL light curves of VER J0521+211 focusing on the data between late 2013 and early 2014.The vertical dashed lines and the hatch fills illustrate the three Bayesian blocks (hatched filled in blue, red, and green, sequentially).The date ranges of the Bayesian blocks are listed in Table 1.The vertical red dotted line illustrates the one night with the highest TeV gamma-ray flux, during which activity at other wavelengths can also be seen.Only 1-σ statistical uncertainties are shown.
Bayesian blocks calculated from the TeV gamma-ray light curve shown in Figure 1.
The GeV gamma-ray flux remained roughly constant in the year 2013, before and after which the flux tapered off.

X-ray and UV/optical Light Curves and Fractional Variability
The light curves measured by the Swift-XRT, the Swift-UVOT, and the Steward Observatory are shown in the lower panels in Figures 1 and  2. The flux of the source at all of the studied wavelengths, as well as the optical polarization, both the electric vector position angle (EVPA) and the polarization fraction, 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.03among 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.
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 0.14 ± 0.01 Steward Observatory (Pol.Frac.) 0.178 ± 0.001 Steward Observatory (EVPA) 0.244 ± 0.001 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.

Optical Polarization Variability
Significant variability in both the electric vector position angle (EVPA) and the polarization fraction were observed.The V -band EVPA roughly ranges between 25 • and 50 • during the time TeV gamma-ray observations were taken.The jet position angle from radio observations was measured in late 2012 and early 2014 to range from around 287 • for the innermost feature at ∼5 mas (from the core) down to about 240 • for the outermost feature at ∼50 mas (Lister et al. 2019).The polarization fraction varied between around 4% and almost 12%.0.178 ± 0.001 and 0.244 ± 0.001, respectively.Changes in polarization fraction were also observed about a month later, with possibly associated variability in X-ray and optical/UV fluxes around MJD 56660, but without significant flux elevation in the VHE band.Therefore, it is inconclusive whether or not the EVPA rotation and polarization fraction changes are associated with the gamma-ray variability.
It is worth noting that the behavior of the optical polarization of this source around the time of high TeV gamma-ray flux is quite different from the dramatic swing of the EVPA during some major gamma-ray flares observed for highpower blazars (see e.g., Abdo et al. 2010).In those cases, a fast swing of the EVPA by 90 • or 180 • accompanied by a decrease (almost to 0) of the polarization fraction was observed during a major GeV gamma-ray flare.However, the relation between EVPA rotations, polarization fractions and gammma-ray flares in BL Lac objects appears to differ from case to case (see e.g., Arlen et al. 2013;MAGIC Collaboration et al. 2019).

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 GeV) [−α−β log 10 (E/300 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.
The MAGIC spectra measured on three of the four nights, 2013 Oct 15, Oct 16, and Nov 29, are shown in Figure 5, together with a logparabola fit.The best-fit spectral parameters and reduced χ 2 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.

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 (α + β 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 "harderwhen-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 −7 ).A linear fit yields a slope of −0.2 ± 0.04.
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  2011), as we used in this work, the frequency of the synchrotron peak would be below 10 15 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 Xray band.Because the extinction correction is the same for all periods, the relative change is robust against the chosen reddening value.

X-Ray/Gamma-Ray Correlation
A moderate correlation between TeV gammaray and X-ray fluxes observed within a oneday 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 makes it hard to interpret the Pearson r correlation.A linear fit, a quadratic fit, and a power-law fit (F TeV ∝ F a X ) were performed.All three models yielded rather large reduced χ 2 values (between 8.4 and 9.6), likely due to the intrinsic scatter of the gammaray/X-ray correlation.The best-fit power-law index is a = 1.55 ± 0.12, indicating the correlation between TeV gamma-ray and X-ray fluxes is likely between quadratic (a = 2) and linear (a = 1).

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 U V W 2 filter), reaching the highest value of 0.76 (a pvalue of 0.0015) between the UV flux observed with the highest-frequency U V W 2 filter and the TeV gamma rays (as shown in the lower right panel of Figure 8).A power-law fit (F TeV ∝ F b UVOT ) indicates a more-than-quadratic relation with indices b > 2 between these bands, although the goodness-of-fit is poor with reduced χ 2 values ranging between 6.4 (between U V W 2 and VHE) and 16.7 (between U V M 2 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 U V W 2 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.

Broadband SED
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 were considered for 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.4 with and without the extragalactic background light (EBL) absorption (Franceschini et al. 2008).The blue crosses and filled diamonds, the red squares and the green filled circles represent the first (BB1), second (BB2, flaring), and the third (BB3) Bayesian blocks, respectively.Data are from Swift-UVOT, Swift-XRT, Fermi-LAT, and MAGIC (only for BB1, shown as blue filled diamonds) and VERITAS, going from the low to high frequencies.As a comparison, the magenta downward triangles and gray upward triangles show the flaring-state and low-state SEDs in 2009 from Archambault et al. (2013).
flux above 100 MeV stayed above five times the value from 3FGL for over five months between 2013 Sep and 2014 Feb, 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 gammaray flux, before increasing again to the high-est 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.

Cross-band Correlations
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 a 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 15 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).
To help constrain the redshift of the source, upper limits were derived from the gammaray 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 VER-ITAS 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 powerlaw 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 gammaray 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 2in Aleksić et al. (2011) where F int (E) is the extrapolated GeV flux, representing the maximally allowed intrinsic flux, at energy E, F obs is the observed TeV gammaray 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 upper 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 (Domínguez et al. 2011;Finke et al. 2010; Franceschini   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.

An SSC Model for the Broadband SED
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.
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 14 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 triangles and gray upward 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 that 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 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 Xray 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-whenbrighter" trend in both bands could be caused by the synchrotron SED peak shifting towards 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 −1 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 −1 first appeared at around 5.5 mas (∼11 pc) in 2012 December (Lister et al. 2018), 4 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 intra-day 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 4 (B/1 G) −1/2 [δ/(1 + z)] −1/2 (ν syn /10 16 Hz) 1/2 ≈ 2 × 10 4 (for BB1/BB3), where B is the strength of the magnetic field, δ is the Doppler factor (δ = [Γ(1−βcosθ)] −1 = 26, where β = (1−1/Γ 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 8 s × (B/1 G) −2 × γ −1 ≈ 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 4 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 inverse-Compton scattering, and the total observed cooling time can be as short as two weeks.Since the dynamic timescale is approximately R/c ∼1.5 day 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 ef-fects, which become substantial when the photon energy hν > m e c 2 in the electron rest frame.Considering the IC scattering between the photons at the synchrotron peak frequency in the comoving frame E peak ≈ 0.1 eV (primed quantities refer to those in the comoving frame) and the electrons that produce the peak synchrotron photons, we have γ syn E peak ≈ 2 keV < m e c 2 .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δγ 2 syn E peak ≈ 4 GeV, below the high-energy SED peak.
From the SED modeling, the SSC component peaks at a frequency ν ssc ≈ 4.4 × 10 24 Hz for the first two blocks, and moves to a lower frequency ν ssc ≈ 2.6 × 10 24 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 (Eq. ( 17) in Tavecchio et al. 1998) , where the spectral indices below and above the IC peak are taken from Fermi-LAT and VERITAS observations as α 1 ≈ 0.9 and α 2 ≈ 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 inverse Compton 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 α 1 .Moreover, the KN limit Doppler factor, defined as Eq. ( 17) in Tavecchio et al. (1998), depends sensitively on α 1 .If we use a value of α 1 ≈ 0.8 as an example, the KN limit Doppler factor becomes about 6, placing the gamma-ray SED peak in the Thomson regime.
Table 4. Parameters used for the one-zone SSC models shown in Figure 9.
Γ 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; log E min , log E max , and log E break are the logarithms of the minimum, maximum, and break electron energies, respectively; p 1 and p 2 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)  The maximum electron energy E e,max ranges between 10 11.9 eV and 10 12.2 eV (γ between about 1.6 × 10 6 and 3.1 × 10 6 ), and the observed maximum gamma-ray energy is E γ ∼ 1 TeV, corresponding to an energy of E γ = E γ (1 + z)/δ ∼ 40 GeV in the comoving frame.Therefore, we can estimate the energy range of the low-energy photons that participate in inverse-Compton scattering and produce the observed gamma rays to be roughly between E γ m 2 e c 4 /E 2 e,max ≈ 0.01 eV and m 2 e c 4 /E γ ∼ 6.8 eV in the comoving frame.These correspond to the observed photon energies between 0.25 and 175 eV, or approximately 6 × 10 13 and 4 × 10 16 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 δm e c 2 /γ max ∼4 eV were up-scattered by the electrons at the maximum energy, the inverse-Compton 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.4eV (corresponding to the frequency of the U V W 2 filter) are probably responsible for the observed VHE flux.
The electron break energy E brk ranges between 10 11.15 eV and 10 11.25 eV (γ between 2.8 × 10 5 and 3.5 × 10 5 ).These electrons emit synchrotron radiation at around 2 × 10 17 Hz (roughly 1 keV), and have a synchrotron cooling time of a few days in the observers' frame.Assuming the mean free path (λ) of these electrons is comparable to their gyroradius of between 3.1 × 10 10 cm and 3.9 × 10 10 cm, we can estimate the electron diffusion coefficient D ≈ λ 2 /(2λ/c) to be a few times 10 20 cm 2 s −1 .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 2 /(4D) to be on the order of 10 5 yr in the comoving frame, much longer than the cooling time.
The light crossing time in the observer's frame is 2R (1 + z)/(cδ) ∼ 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 −11 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 (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 due 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.

SUMMARY
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 −7 photon m −2 s −1 (∼37% of the Crab Nebula flux), and the monthly GeV gamma-ray flux above 100 MeV remained higher than seven times the 3FGL catalog value for a period that lasted one year.
The fractional variability amplitude F var was observed to be the lowest for the lowestfrequency V -band flux, and increases with fre-quency to the X-ray band, and then decreases with frequency to the GeV energies, before increasing again to the highest at the highest TeV energies.Such a pattern of high variability amplitude F var being observed at frequencies above the two peaks in the SED is consistent with previous observations of gamma-ray blazars (see e.g., Baloković et al. 2016, and references therein).
A moderate more-than-linear correlation between the X-ray and TeV gamma-ray fluxes was observed, with a Pearson correlation coefficient of 0.7.The X-ray spectrum appeared harder when the flux was higher.Unlike the gamma-ray and X-ray band, the optical flux did not increase significantly during the studied period compared to the archival low-state flux.The combination of low optical variability amplitude, the higher X-ray variability amplitude, and the "harder-when-brighter" trend in the X-ray spectrum indicates that the synchrotron peak of the SED may become broader during flaring states of this source, and the synchrotron peak frequency can shift to higher frequencies.Such a behavior of the SED can be explained by activity from particles with the highest energies.Of particular interest, the Xray spectrum in 2013 and 2014 is much softer than that observed in a flaring state in 2009 (Archambault et al. 2013), while the normalization values are comparable, suggesting strong X-ray spectral variability on timesales of years.
The correlation between the optical and VHE fluxes is weaker than that between X-ray and VHE fluxes, while that between UV and VHE fluxes is the strongest and follows a more-thanquadratic relation, possibly from the combination of a varying particle density and an expansion of the emitting region (Katarzyński et al. 2005).The weaker correlation between the optical and VHE fluxes, as well as the low optical variability observed from 2009 to 2014, could also suggest that the optical emission from the source originates from a different region (see e.g., Rajput et al. 2019).
Upper limits at the 95% confidence level on the redshift of the source were derived following Georganopoulos et al. (2010) and Aleksić et al. (2011) to be z ≤ 0.3, z ≤ 0.278, z ≤ 0.308, and z ≤ 0.281 using four EBL models from Domínguez et al. (2011), Finke et al. (2010), Franceschini et al. (2008), and Gilmore et al. (2012), respectively.These upper limits are more constraining than the previously reported value of z < 0.34 (Archambault et al. 2013).A precise redshift measurement of this source in the future will be important to studies of the particle distribution as well as the EBL intensity.
The broadband SED can be adequately described with a one-zone synchrotron self-Compton model.The inverse-Compton scattering in the source happens in the Thomson regime below the high-energy SED peak and likely transitions to the KN regime at or above the high-energy SED peak.In contrast to the hard X-ray spectrum and HBL-like SED properties of VER J0521+211 observed during the flaring state in 2009, the soft X-ray spectrum, the SED synchrotron peak frequency, the "harderwhen-brighter" trend in both UV and X-ray band, and the radio morphology from the MO-JAVE program are all consistent with the source being an IBL object during the flaring state in 2013.
As a rare IBL that can be detected at TeV gamma-ray energies on timescales of months, more simultaneous MWL observations of VER J0521+211 can provide valuable insights into this type of source.Of particular interest are the gamma-ray band and X-ray band, as they have exhibited more variability compared to other energies.Continued radio interferometry monitoring is also interesting for IBLs, as there could be a potential association between the ejection of superluminal knots and gamma-ray flares.

Figure 3 .
Figure 3. Fractional variability amplitude (F var ) as a function of frequency computed from the light curves shown in Figure 1.The filled circle is calculated from the nightly binned VHE light curve, the filled square is from the monthly binned Fermi-LAT light curve, the filled diamond is from the Swift-XRT light curve binned by observations, the filled crosses are from the Swift-UVOT light curves with the six filters binned by observations, and the filled plus sign and downward triangle are from the Vband polarization fraction and EVPA measured by the Steward Observatory, respectively.

Figure 4 .
Figure 4.The VERITAS gamma-ray spectra of VER J0521+211 during three different flux states.The dashed lines show the best-fit log-parabola models for the spectra.The 1-σ uncertainties from the model parameters are shown in shaded colors.These flux states were chosen based on a Bayesian block analysis of the VERITAS TeV gamma-ray light curve.

Figure 5 .
Figure 5.The MAGIC gamma-ray spectra of VER J0521+211 from observations on three nights in 2013.The dashed lines show the best-fit logparabola models.

Figure 6 .
Figure6.Scatter plot to investigate the correlation between the best-fit X-ray photon index at 2 keV and the X-ray energy flux.The blue dashed line shows the best-fit linear model.

Figure 7 .
Figure 7. Scatter plot to investigate the correlation between the TeV gamma-ray flux and the X-ray energy flux.The flux measurements in the two bands were performed on the same nights.The gray dashed and dotted lines show the best linear and quadratic fits, respectively.The solid blue line shows the best power-law fit.

Figure 9 .
Figure 8.The correlations between the TeV gamma-ray flux and the optical/UV energy fluxes measured with the six Swift-UVOT filters.The frequency of the Swift-UVOT observations increases from the upper left panel to the lower right panel.Similar to Figure 7, the flux measurements in the two bands in each panel were performed on the same nights.The blue dashed lines show the bestfit power law to guide the eye.

Figure 10 .
Figure 10.(Upper) The optical depth τ γγ (E) of the gamma rays due to γγ pair absorption induced by EBL photons.The blue filled crosses, the red filled circles, and the green filled diamonds show the 95% upper limits derived from VERI-TAS observations during periods BB1, BB2, and BB3.The optical depth τ γγ (z, E) from four EBL models, Domínguez et al. (2011) (blue solid line), Finke et al. (2010) (orange dashed line), Franceschini et al. (2008) (copper dotted line), and Gilmore et al. (2012) (pink dash dotted line), evaluated at the maximum redshift allowed by the 95% upper limits are also shown.(Lower) The minimum difference between the 95% upper limits on the EBLinduced pair absorption optical depth τ γγ (E) calculated from VERITAS data and the optical depth τ γγ (z, E) calculated from EBL models as a function of redshift.
This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation and the Smithsonian Institution, by NSERC in Canada, and by the Helmholtz Association in Germany.This research used resources provided by the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy's Office of Science, and resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument.
Figure 1.MWL light curves of VER J0521+211 of the entire period of interest in this work, from late 2012 to early 2014.The vertical dashed lines and the colored hatch fills illustrate the three Bayesian blocks between late 2013 and early 2014, based on which we selected the three periods for SED modeling.The binning of the light curves is nightly for VERITAS and MAGIC, monthly for Fermi-LAT, and by observations for Swift and Steward Observatory.Only 1-σ statistical uncertainties are shown.

Table 2 .
Fractional variability calculated from the MWL light curves shown in Figure1.

Table 3 .
VHE gamma-ray spectral fit parameters using power-law and log-parabola models.
for the low-state SED in 2009 are also included for a comparison.