Millimeter-VLBI Detection and Imaging of the Gravitationally Lensed γ-Ray Blazar JVAS B0218+357

Between 2017 and 2018 January, we observed JVAS B0218+357 with KaVA for nine sessions as part of coordinated MWL campaigns on this source. Each session lasted two consecutive nights, where we had a 5–8 hr track at 22 GHz for the first day and another similar track at 43 GHz for the following day. By default seven stations (three from KVN and four from VERA) joined each session, but occasionally VERA-Mizusawa or VERA-Ishigaki was missing due to local issues (see Table 1). All the data were recorded at 1 Gbps (a total bandwidth of 256 MHz with eight 32 MHz subbands). Left-hand circular polarization was received and the data were correlated by the Daejeon hardware correlator (Lee et al. 2015). The initial data calibration (amplitude, phase, bandpass) was performed using the National Radio Astronomy Observatory Astronomical Image Processing System (AIPS; Greisen 2003) based on the standard KaVA/VLBI data reduction procedures (Niinuma et al. 2014; Hada et al. 2017; Park et al. 2019). The data were finally averaged in frequency, but individual subbands were kept separate to avoid bandwidth smearing. Similarly, the data were time-averaged to only 16 s to avoid time smearing. Following the initial calibration, imaging was performed in the Difmap software (Shepherd et al. 1994) using the standard CLEAN and self-calibration procedures. The visibility amplitude of KaVA data was self-calibrated down to a solution interval of 10–120 minutes.

In this paper, we present the KaVA results for four out of the nine sessions where we had simultaneous KVN 86 GHz data (see below). Our regular monitoring of JVAS B0218+357 with KaVA was continued until early 2019 and detailed analyses of the structural evolution and time-domain properties using the whole KaVA data sets will be presented in future papers.

Truly in parallel to five out of the nine KaVA 43 GHz sessions, we also observed the source with the KVN-only array with 43 GHz/86 GHz dual-frequency simultaneous recording mode. These five dates are 2017 January 15 (MJD 57768), 2017 October 17 (MJD 58043), 2017 November 12 (MJD 58069), 2017 November 25 (MJD 58082), and 2018 January 5 (MJD 58123). For the session in 2017 November 25, only two stations were available due to an antenna motor issue at the Ulsan station, so we excluded this epoch from our study. A wideband 4 Gbps mode was used where each frequency band was recorded at 2 Gbps (a bandwidth of 512 MHz for each band). We first calibrated the 43 GHz data in the standard manner in AIPS. The derived 43 GHz fringe phase solutions were then transferred to the 86 GHz data by using the frequency-phase transfer technique (e.g., Rioja et al. 2011; Algaba et al. 2015; Zhao et al. 2018). This greatly reduced the rapidly fluctuating phase components at 86 GHz and allowed us to perform fringe fitting with a much longer solution interval than the typical coherence time at this frequency. Thanks to this strategy we are able to detect 86 GHz fringes for the target at sufficient signal-to-noise ratios (S/Ns) for most of the scans. The KVN data were averaged in frequency and time in the same manner as for the KaVA data. Imaging was performed in Difmap. Since KVN is a three-element array, we performed phase-only self-calibration using a solution interval of 30 s.

In Figure 1 we show the representative uv-coverage of our KaVA/KVN observations at each frequency. In Table 1 we summarize the basic information of individual KaVA/KVN images used in this paper. Typical angular resolution of KaVA (a maximum baseline length D = 2300 km) is 1.2 mas (22 GHz) and 0.6 mas (43 GHz), while that of KVN (D = 560 km) is 1 mas at 86 GHz. We assume the typical flux calibration accuracy of 10% and 20% for KaVA 22/43 GHz and KVN 86 GHz data, respectively.

Zoom In Zoom Out Reset image size

In Figure 2 we show a wide-field image of the lensed system obtained with KaVA at 22 GHz. The lensed images A and B are clearly detected, while the diffuse Einstein ring emission is totally resolved out in our VLBI imaging due to the lack of short baselines sensitive to such a large-scale component. In Figure 3, we show zoom-in milliarcsecond-scale images toward A and B at 22, 43, and 86 GHz, resolving the magnified parsec-scale structures in both of the lensed images. These images are produced by stacking over the four sessions and better characterize the morphology at each frequency. At 22 GHz the lensed image A shows a core–jet structure toward the north–northeast, while B exhibits a core–jet profile toward the east. The lensed image A additionally shows some diffuse extension in the direction perpendicular to the jet axis. These characteristics observed in our KaVA 22 GHz images are in good agreement with those seen in the previous VLBA 22 GHz images (Spingola et al. 2016). Hereafter we label A1 (core in A), A2 (jet component in A), B1 (core in B), and B2 (jet component in B), respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the other hand, the overall source brightness is significantly (a factor of 2) lower than the past levels. In Figure 4, we show a long-term 15 GHz light curve of the source obtained by the Owens Valley Radio Observatory (OVRO) 40 m telescope (Richards et al. 2011). As can be seen, the historical radio flux of the target was gradually decreasing between MJD 55500 and MJD 58300, and the period of our KaVA observations was roughly coincident with a historical low state of the source, although some local fluctuations (∼5% in rms fractional amplitude) on timescales of months were seen.

Zoom In Zoom Out Reset image size

At 43 and 86 GHz, the results presented here are the first robust VLBI detection and imaging of the lensed system (but see also Porcas & Patnaik 1996; Lee et al. 2008). In the KaVA 43 GHz images the core–jet structure is more clearly resolved and each lensed image can be well characterized by a bright core and a secondary component separated by ∼1.5 mas (in A) and ∼1.3 mas (in B) from the core. At 86 GHz, VLBI observations are generally challenging, but the unique capability of the KVN wideband multifrequency receiving system allows us to detect and image the source at adequate S/N. For the lensed image A (typically S/N > 30), the core–jet morphology toward the northeast is consistently seen up to 86 GHz. The lensed image B is detected at S/N ∼ 10 in each epoch. Its morphology is dominated by a single component, but there is a clear elongation toward the east, which is consistent with the structure seen at 22 and 43 GHz.

In Figure 5, we show the observed integrated radio spectra of A and B up to 86 GHz. For comparison, the low-frequency spectra from the literature (Mittal et al. 2006; Spingola et al. 2016) are also plotted. One can see that the (integrated) spectra for both A and B tend to become steeper with increasing frequency. Between 22 and 43 GHz both A and B show very stable spectral shapes over the four epochs with an averaged spectral index of α_A,KQ = −0.48 ± 0.05 or α_B,KQ = −0.49 ± 0.05, respectively. At 86 GHz the emission is quite variable and the spectral slope significantly changes with time. For the lensed image A, the integrated spectral index between 43 and 86 GHz varies between −0.58 and −1.70 with a mean value of α_A,QW ∼ −1.28. For the image B, α_B,QW is generally steeper than α_B,KQ with a mean slope of α_B,QW ∼ −2.0. However, the image quality of B at 86 GHz is rather limited and some amount of low-level flux (from the extended emission associated with the jet) in B may not be recovered at 86 GHz.

Zoom In Zoom Out Reset image size

In Table 2, we also show the spatially resolved spectral indices for the individual features A1, A2, B1, and B2. As typical for a blazar, the jet (A2 and B2) has steeper spectra than the core (A1 and B1). Nevertheless, it is notable that the core spectra themselves are substantially steep especially between 43 and 86 GHz.

Table 2.  Spatially Resolved Spectral Indices

Feature	${\alpha }_{\mathrm{KQ}}$ (22–43 GHz)	αQW (43–86 GHz)
A (total)	−0.48 ± 0.05	−1.28 ± 0.42
A1 (core)	−0.45 ± 0.17	−0.76 ± 0.30
A2 (jet)	−0.53 ± 0.11	−2.51 ± 0.85
B (total)	−0.49 ± 0.05	−2.00 ± 0.25
B1 (core)	−0.32 ± 0.15	⋯
B2 (jet)	−0.53 ± 0.10	⋯

Note. Each spectral index is a mean value over the four epochs. The error of each spectral index is the standard deviation over the four epochs. α_QW in B was estimated only for the total flux due to the low image S/N of B at 86 GHz.

Download table as:  ASCIITypeset image

In Figure 6, we plot the observed magnification ratio A/B (i.e., (A1+A2)/(B1+B2)) as a function of frequency. The flux ratio can be determined much more accurately than the individual fluxes since any systematic flux calibration errors for A or B are canceled out. To obtain a wider frequency coverage of A/B ratio, we also include the previous results obtained at 22 GHz and lower frequencies (Mittal et al. 2006; Spingola et al. 2016). The previous observations between 1.7 and 22 GHz reported a gradual increase of the ratio with frequency, from ∼2 at 1.7 GHz to ∼3.4–4 at 22 GHz. Our magnification ratio measurements at 22 GHz result in ∼3.4, which is in good agreement with the previous results. At 43 GHz, the measured ratios are similar to the ones at 22 GHz with a slightly larger scatter. At 86 GHz, the scatter is still very large and the apparent large ratios are likely caused by the insufficient flux recovery of lensed image B, so the obtained ratios at 86 GHz are less constrained than at 22/43 GHz and should be considered as upper limits.

Zoom In Zoom Out Reset image size

Using the KaVA 43 GHz images where the core and jet are well resolved, we also checked the magnification ratio for the core (A1/B1) and jet (A2/B2) separately. We find that A1/B1 tends to be systematically larger than A2/B2, such that A1/B1 = 3.3–4.5 (with a mean value 3.7) while A2/B2 = 3–3.3 (with a mean value 3.2). We will discuss the possible reasons for the difference in Section 4.2.

The gravitationally lensed blazar JVAS B0218+357 provides a rare opportunity to study the physics of VHE γ-ray-emitting relativistic jets at z ∼ 1. High-resolution radio observations are especially useful to directly resolve the detailed structure of the innermost jet regions. Previous extensive VLBI imaging of this source  (e.g., Biggs et al. 2003; Mittal et al. 2006; Spingola et al. 2016) revealed a core–jet morphology (in both A and B) at parsec scales that is typical for powerful blazars. However, Mittal et al. (2007) found that the low-frequency (1–15 GHz) radio spectra of the lensed image A (and accordingly the A/B flux ratios) were significantly affected by the free–free absorption (FFA) in the lens galaxy, preventing us from imaging the pure GL morphology. In addition, (even though the external effects are irrelevant) the jet base could become opaque by the effect of synchrotron-self-absorption (SSA), which might make the access to the high-energy γ-ray emission site difficult at low frequencies (e.g., Acciari et al. 2009; Spingola et al. 2016).

The KaVA/KVN observations presented here have for the first time revealed the detailed milliarcsecond-scale structures of this lensed system at millimeter wavelengths. In particular, we have found that the simultaneous millimeter spectra of the radio core and jet become progressively steep with frequency. This indicates that the radio structure becomes substantially transparent to any absorption effects at these frequencies. In fact, if we adopt the best-fit FFA parameters derived by the low-frequency radio spectra toward A (e.g., ${T}_{{\rm{e}}}={10}^{4}\,{\rm{K}}$ and $\mathrm{EM}=1.8\times {10}^{7}\,{\mathrm{cm}}^{-6}\,\mathrm{pc}$ where ${T}_{{\rm{e}}}$ and EM are the temperature and emission measure of the presumed ${{\rm{H}}}_{\mathrm{II}}$ region in the lens galaxy; Mittal et al. 2007), the FFA optical depth at ≥43 GHz results in ${\tau }_{\mathrm{FFA}}\leqslant 2\times {10}^{-3}$, being essentially negligible. Therefore, we start to see the absorption-free radio structure of JVAS B0218+357 at these frequencies. Note that we do not rule out the possible presence of a highly SSA-thick substructure in the core (i.e., a spectral turnover ${\nu }_{\mathrm{SSA}}\,\gt 86$ GHz). However, the contribution from such an inverted spectral component should be tiny at the observed frequencies.

It is notable that the millimeter-core spectra between 43 and 86 GHz are quite variable, while those between 22 and 43 GHz are rather stable over a year (Figure 5). If we attribute the observed spectral steepening at the high frequencies to the radiative cooling of synchrotron electrons, we can estimate the magnetic field strength of the radio-emitting region through B ∼ $1.3\,{t}_{\mathrm{syn},\mathrm{obs}}^{-2/3}{\nu }_{{\rm{b}},\mathrm{obs}}^{-1/3}{\delta }^{-1/3}{(1+z)}^{1/3}$ (Gauss), where B, ${t}_{\mathrm{syn},\mathrm{obs}}$ and ${\nu }_{{\rm{b}},\mathrm{obs}}$, δ are the observed timescale of synchrotron cooling (in yr) and observed break frequency (in GHz) and Doppler factor, respectively (e.g., Pacholczyk 1970). The 86 GHz emission is clearly variable between 2017 January and October (0.75 yr) and possibly even on shorter timescales of ∼1 month (2017 October–November). This suggests that B of the 86 GHz emission site is at least B ∼ $0.33{\delta }^{-1/3}$ G (for ${t}_{\mathrm{syn},\mathrm{obs}}=0.75\,\mathrm{yr}$) and possibly as large as ∼$1.5{\delta }^{-1/3}$ G (for ${t}_{\mathrm{syn},\mathrm{obs}}=0.08\,\mathrm{yr}$). The Doppler factor δ of the radio emission site is still not well constrained since the jet morphology is largely stationary, but if we assume that δ ∼ 20 (adopted in a spectral energy distribution (SED) model by Ahnen et al. 2016) is a maximum value, we can obtain B ∼ 0.12–0.55 G.

Interestingly, the derived B-field strength seems to be somewhat larger than that implied for the γ-ray-emitting site in the previous SED modeling (0.03 G; Ahnen et al. 2016). In fact, they indicate that the broadband radio-to-γ-ray emission cannot be modeled by a simple single-zone model but can be better described by a two-zone model where the γ-ray emission region is offset from the low-energy emission region. It should be noted that the above-mentioned model also did not explain the emission observed from this source at $\lesssim 100\,\,\mathrm{GHz}$. Our B-field estimate therefore seems to be consistent with this picture, implying that the primary millimeter emission region may originate in a more magnetized part of the jet.

Another remarkable feature seen in the data presented here is that the core–jet morphology in each lensed image remains extremely stable at least over two decades since the first VLBI images of this system were obtained (Patnaik et al. 1995). As discussed in Section 4.2, the projected distance of the jet feature from the core is estimated to be ∼10 pc in the source plane, implying a deprojected (viewing-angle-corrected) distance from the core to be of the order of ∼100 pc. This feature is reminiscent of a "recollimation shock" claimed in other relativistic jets exhibiting active γ-ray flares (e.g., Cheung et al. 2007; Marscher et al. 2008; Hada et al. 2018). If this is the case, the distant component could also be relevant to the γ-ray production. To test this hypothesis, it would be useful to cross-correlate the γ-ray light curves with spatially decomposed (core and jet) radio light curves.

Having absorption-free VLBI images of JVAS B0218+357, now it would be instructive to apply for a simple lens model. This helps to discuss the intrinsic source geometry and also to check if the observed properties can be explained by a simple model or if any additional consideration is required. Here, in order to recover the intrinsic position (therefore, magnification) of the core and jet, we model JVAS B0218+357 adopting a parametric lens modeling approach using gravlens (Keeton 2001). We backward-ray trace the emission from this blazar by using average positions and flux ratios from the KaVA observations at 43 GHz (those at the highest possible angular resolution; see Table 3), because these data clearly resolve the subcomponents (see Figure 3, second row). We keep the lens mass model parameters fixed to the values obtained by Wucknitz et al. (2004), who inferred the gravitational potential at high precision using the entire extended emission from the lensed images and Einstein ring (see also Wucknitz 2002; Biggs et al. 2003). In particular, we adopt their singular elliptical power-law model, without external shear,⁸ where the mass strength is 01616; the ellipticity and its position angle are 0071 and −47° (east of north), respectively; the slope of the power-law is −1.04; and the lensing galaxy is at (0.2582, 0.1210) arcsec with respect to image A1, which is at (0, 0) arcsec. In order to estimate the uncertainty on the recovered parameters (source position, time delays, and image magnifications), we use Monte Carlo simulations using the same methodology adopted by Spingola et al. (2019) and Spingola & Barnacka (2020). This method takes into account both the observational uncertainties on positions and fluxes (reported in Table 3), and those of the lens mass model parameters (Wucknitz et al. 2004). Lens and source plane obtained using this method are shown in Figure 7 and the fitted parameters are summarized in Table 4.

Zoom In Zoom Out Reset image size

We find that (in the source plane) the core and jet features of JVAS B0218+357 are separated by 1.4 ± 0.2 mas (10.4 ± 1.6 pc) in projection, and they are located at about 110 mas from the center of mass of the lensing galaxy, as shown in Figure 7. The expected time delays between images A1–B1 and A2–B2 are found to be consistent with those measured by Biggs & Browne (2018) and Cheung et al. (2014) within the uncertainties. As for the flux magnification, the lens model finds ${\mu }_{{\rm{A}}1}=2.40\pm 0.05$ and ${\mu }_{{\rm{B}}1}=0.48\pm 0.01$ for A1 and B1, respectively. Therefore, the model-predicted magnification ratio for the lensed images of the core results in 5. This is broadly in agreement with the observed values (A1/B1 ∼ 3.3–4.5; Section 3.3), although a mean value in Table 3 is somewhat smaller than the predicated one. As for the jet part (lensed images A2 and B2), the model-predicted magnification ratio is 4.8, consistent with that of the core. However, the observed ratios for the jet component (A2/B2 ∼ 3–3.3) appear to be quite smaller than the predicted one. Therefore, the actual spatial gradient in the magnification ratio seems to be larger than predicted in the simple lens model applied here.

Given the absorption effects being negligible at our observing frequencies, an intriguing possibility for explaining the difference between predicted and observed magnification ratios is "substructure lensing" by a compact clump in the lens galaxy (especially toward A). Such substructure lensing can locally affect individual parts of the background differently, leading to an additional magnification (or de-magnification) of the lensed source. For the lens geometry of JVAS B0218+357, a compact structure with projected size of the order of ∼10 pc in the lens plane is enough to affect the core and jet components of a radio lensed image. Sizes of tens of parsecs are typical for a giant molecular cloud (GMC; see, e.g., Chevance et al. 2020). Possible clumps in GMC can result in additional moderate amplifications by a factor of 1.5 at perpendicular distance scales of parsecs (see Sitarek & Bednarek 2016). This would be able to explain the apparent discrepancy between the lens model predictions and observations.

Alternatively, we also note that since image A is highly distorted by the lensing effect, the surface brightness in some stretched part of the jet (A2) could become either too faint to be detected by our VLBI sensitivity, or too extended to be sampled by the shortest baselines of the array. In this case we would recover only a fraction of the extended flux density that should be associated with A2, while the core (A1) should be unaffected because of its intrinsic compactness.

Modeling the whole details of the observed lens properties is beyond the scope of this paper. Nevertheless, the discussion presented here indicates a good capability of high-frequency VLBI data for studying GL at the smallest scales. A sophisticated lens modeling, including the absorption and substructure effects at the lens, together with deeper millimeter-VLBI images will allow us to better constrain the lensing potential and accordingly the physical properties of the innermost jet regions in distant AGN.