Free Access
Issue
A&A
Volume 649, May 2021
Article Number A153
Number of page(s) 9
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202039616
Published online 01 June 2021

© ESO 2021

1. Introduction

The most accepted scenario to explain the formation of relativistic jets in active galactic nuclei (AGNs) involves the extraction of energy from the innermost regions of the black hole (BH) accretion disk (Blandford & Payne 1982) and/or the ergosphere of the BH itself (Blandford & Znajek 1977) by large-scale magnetic fields anchored into the magnetosphere of the rotating BH. This energy, in the form of Poynting flux, is then converted into the kinetic energy of a collimated jet, which gradually accelerates until it reaches its asymptotic Lorentz factor where the electromagnetic and kinetic energy fluxes reach equipartition (Vlahakis & Königl 2004). In the jet acceleration region, usually denoted as the acceleration and collimation zone (ACZ; e.g., Marscher et al. 2008), the jet assumes a parabolic shape (or sometimes even a more collimated cylindrical shape).

Such a seemingly simple model, however, has been subject to continuous debate, in part due to insufficient observations but also due to the complexity introduced when the mean particle energy is relativistic. How and where the electromagnetic energy is transferred to the jet is an open question. Investigating the transverse jet size and velocity structure at the jet base through Very Long Baseline Interferometry (VLBI) studies is an important tool for studying the acceleration mechanisms internal to jets since the acceleration and the geometry of the jet are strictly connected (e.g., Vlahakis & Königl 2004; Beskin & Nokhrina 2006; Lyubarsky 2009).

In particular, millimeter-VLBI studies (see Boccardi et al. 2017 for a review) provide high-resolution images of regions deep in the jet that are otherwise optically thick at longer centimeter wavelengths. Millimeter-wave VLBI images, therefore, give us a direct probe of the inner geometry and velocity structure of jets in different AGN classes: radio galaxies (Cygnus A – Boccardi et al. 2016a,b; M 87 – Hada et al. 2013, Kim et al. 2018; NGC 4261 – Nakahara et al. 2018; NGC 1052 – Baczko et al. 2016, 2019; Nakahara et al. 2020), a narrow-line Seyfert 1 (NLS1) galaxy, 1H 0323+342 (Hada et al. 2018), and a few blazars (Mrk 501 – Giroletti et al. 2004; 3C 273 – Akiyama et al. 2018; TXS 2013+370 – Traianou et al. 2020). The picture is further complicated by synchrotron opacity effects (Blandford & Königl 1979) that can obscure the inner jet and shift the apparent location of the VLBI radio core at large distances from the BH, ∼10(4 − 6) Schwarzschild radii (Rs; e.g., Marscher et al. 2008), as well as by relativistic effects when the jet is oriented at small viewing angles (as is the case in blazars, i.e., BL Lacs and flat spectrum radio quasars) and the VLBI core emission typically dominates over the fainter extended jet emission.

Relativistic jets should also be considered in the context of their surrounding environments, the accretion flow, and the surrounding disk winds and interstellar medium. The wind (which may be magnetized) coming from the outer parts of the accretion disk may play a crucial role in the confinement of jets (McKinney 2006; Tchekhovskoy et al. 2008). Boccardi et al. (2020) show that in high-excitation radio galaxies (HERGs), the ACZ, expressed in Schwarzschild radii, extends over larger distances than it does in low-excitation radio galaxies (LERGs) and that this may be associated with the different accretion mechanisms in these two classes of objects. A HERG, hosting radiatively efficient standard disks, would have a stronger disk wind component that keeps the jet collimated for larger distances. This result is in agreement with the behavior observed in the beamed subset of HERGs and LERGs by Potter & Cotter (2015). They find that flat spectrum radio quasars (FSRQs), the beamed subgroup of HERGs, accelerate over longer distances than BL Lacs, which are instead considered the beamed subgroup of LERGs.

Recently, Kovalev et al. (2020), and previously Pushkarev et al. (2017), investigated the jet geometry of a large sample of AGNs using mainly 15 GHz Very Long Baseline Array (VLBA) data from the MOJAVE program. In Kovalev et al. (2020), they found that ten nearby AGNs (z < 0.07) show a transition from a parabolic to a conical shape in the scales probed at 15 GHz. Among these sources is the jet of BL Lacertae (or BL Lac, the eponymous blazar of the BL Lac objects class), whose transition occurs at a linear distance of ∼2.5 mas from the BH (∼24.6 pc de-projected, considering a viewing angle of 7.6°). In contrast, the earlier study by Pushkarev et al. (2017) presented a different result: The jet expands conically from 0.6 mas to 10 mas, thus up ∼100 pc de-projected, and no break is observed. In the Pushkarev et al. (2017) case, the jet is not accelerating and has probably ceased to be Poynting flux dominated. Intending to clarifying our understanding of the collimation profile in BL Lac, we performed a detailed study of the jet geometry, now also adding the information from higher-frequency VLBI observations (43 and 86 GHz). Thanks to the reduced opacity and the higher resolution gained by our data set with respect to the 15 GHz data (approximately three and six times higher for the 43 and 86 GHz observations, respectively), we can investigate the innermost portion of the jet in more detail. We demonstrate how higher-resolution millimeter VLBI data (with a high dynamic range) are essential for determining the jet’s collimation profile in the innermost jet regions.

In Sect. 2 we present the data set and the methods used for the analysis, and in Sect. 3 we discuss our results, contextualizing them within the jet kinematics and the physical conditions of the external medium. Our conclusions are summarized in Sect. 4. We adopt a flat cosmology with Ωm = 0.27, ΩΛ = 0.73, and H0 = 71 km s−1 Mpc−1 (Komatsu et al. 2009). Assuming these values, at the redshift of BL Lac (z = 0.0686), 1 mas corresponds to a linear distance of 1.296 pc.

2. Data and methods

We analyzed 86 GHz Global Millimeter VLBI Array (GMVA) data of BL Lac obtained from May 2009 to April 2017 (11 observing epochs) and 43 GHz VLBA data from the VLBA-BU-BLAZAR program1 covering 10 years of observations (94 epochs), obtained between May 2009 and March 2019. The GMVA data are part of an ongoing monitoring program 2 (PI: A. Marscher) of a sample of γ-ray bright blazars and radio galaxies. The data reduction of both data sets is available in Casadio et al. (2019) and references therein.

We used Difmap to fit the complex visibilities with a model source consisting of components described by circular Gaussian brightness distributions. We obtained, for each epoch, a model fit that provides information about the full width at half maximum (FWHM) size, the flux density, the distance, and the position angle relative to the core (considered stationary) of each component. For the relative uncertainties, we followed Casadio et al. (2019). Model-fit parameters for all components and epochs at 86 GHz, as used in the analysis, are reported in Table 1.

Table 1.

Model-fit parameters for the 86 GHz GMVA observations.

Blazars are highly variable sources, and this is mostly reflected in their morphology and behavior at millimeter wavelengths. Therefore, to obtain a higher dynamic range image and to better recover the entire jet cross section, we stacked all the total intensity images together. Before stacking, we convolved the 86 and 43 GHz images with 0.10 mas and 0.19 mas FWHM restoring circular beams, respectively. The alignment of the images was done using the position of the intensity peak. Our radio map stacking, therefore, implicitly assumes that the radio core in BL Lac is consistently the brightest feature in the jet at all epochs. While this would not necessarily be true for an FSRQ (e.g., 3C 273), this is indeed generally the case in BL Lac type objects. The resulting stacked images at both frequencies are shown in Figs. 1 and 2, respectively.

thumbnail Fig. 1.

Total intensity stacked image of BL Lac at 43 GHz. A total of 94 images obtained between May 2009 and March 2019 are used. The restoring circular beam (lower left circle) is 0.19 mas. Contours are drawn at 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, and 64% of the peak (1.71 Jy beam−1). The total intensity ridge line is shown in red.

thumbnail Fig. 2.

Total intensity stacked image of BL Lac at 86 GHz. Eleven images obtained between May 2009 and April 2017 are used. The restoring circular beam (lower left circle) is 0.10 mas. Contours are drawn at 0.35, 0.5, 1, 2, 4, 8, 16, 32, and 64% of the peak (0.847 Jy beam−1). The total intensity ridge line is shown in red.

We measured the radial jet width in the stacked images using two different methods: a Python (Van Rossum & Drake 2009) script where the jet ridge line is obtained first and perpendicular slices are then taken to obtain the transverse jet brightness profile at every cut; and the tasks SLICE and SLFIT in the Astronomical Image Processing System (AIPS; Greisen 1990), where slices instead have all the same orientation, perpendicular to the main jet direction. The Python script (the following packages are included: matplotlib 3, numpy 4, and scipy 5) should be more sensitive to jet bending, which is only marginally observed in the considered images. The two methods provide consistent results, and the comparison between them was used to test the robustness of the obtained jet profile. Here we will only focus on the first method; the results obtained using the AIPS procedure are outlined in Appendix A.

We fit every slice brightness profile with a Gaussian distribution, and we obtained the deconvolved jet width (w) by subtracting from the Gaussian FWHM the FWHM of the restoring beam (b) in quadrature (). Following Pushkarev et al. (2017), we obtained the uncertainties on the jet width measurements by rotating the position angle of each slice by ±15° in steps of 2° and using the new Gaussian profiles to compute the root-mean-square noise.

In Marscher et al. (2008), the authors located the 43 GHz VLBI core at ∼0.91 pc (de-projected, assuming a 7° viewing angle) from the BH, just after the end of the ACZ that terminates at ∼0.65 pc from the BH. Given the inverse dependence of the position of the radio core from the jet base with the frequency (rcore(ν)∝ν−1) found in BL Lac, and a core shift of ∼0.04 mas between 22 and 43 GHz (O’Sullivan & Gabuzda 2009), the core position offset between 43 and 86 GHz is expected to be small. Hence, assuming a negligible core shift between these two frequencies, we displaced the VLBI cores at both frequencies by 0.91 pc, and we estimated the radial de-projected distance of each slice from the BH using the same viewing angle as in Marscher et al. (2008).

3. Results and discussion

3.1. The jet collimation profile

In Figs. 3 and 4 we analyze the dependence between the jet width, w, and the de-projected distance, d, from the BH (in pc) for the 43 and 86 GHz profiles, respectively. It is evident from both profiles that in neither is it possible to fit the w − d dependence with a single power law. The 86 GHz jet expansion profile exhibits a decrease in the jet width around (0.3–0.4) mas followed by a hint of a more rapid expansion. Indeed, after almost 0.4 mas, as observed at 43 GHz, the jet expands more rapidly and recollimates later, around 1.5 mas, where it restarts its expansion with a rate similar to that at the beginning. Therefore, to obtain the information on the overall geometry of the jet, we fit the 86 GHz profile only up to 0.25 mas, and we fit the parts before 0.5 mas and after 1.5 mas at 43 GHz separately. Using the least-squares method, we fit the w − d dependence at 86 GHz and the first 0.5 mas at 43 GHz with a power law of the type w = a + b ⋅ dk, where a gives the jet width at the location of the core. For the 43 GHz jet beyond 1.5 mas, we used a power law fit of the type w = b ⋅ dk. The best-fit parameters and relative uncertainties obtained are listed in Table 2.

thumbnail Fig. 3.

Jet width versus radial distance at 43 GHz, calculated with respect to the VLBI core (in mas) and to the BH (in de-projected parsecs). The blue and orange lines are the best fit of assumed power laws of the type w = a + b ⋅ dk and w = b ⋅ dk, respectively. The best-fit parameters with their uncertainties are also reported in Table 2.

thumbnail Fig. 4.

Jet width versus radial distance at 86 GHz, calculated with respect to the VLBI core (in mas) and to the BH (in de-projected parsecs). The blue line is the same as in Fig. 3.

Table 2.

Best-fit parameters obtained for the power law equations w = a + b ⋅ dk and w = b ⋅ dk fitting the jet width–distance dependence at 43 and 86 GHz.

The exponent k is expected to be 0.0 for a cylindrical jet, close to 0.5 in the case of parabolic expansion, and 1.0 for a conically expanding jet. We obtained k ≃ 1 for the three fits, which means that the overall jet geometry in BL Lac, from ∼0.9 to ∼30 pc (de-projected), is conical.

The derived values for a, for the 86 and 43 GHz cases, are a = (31 ± 2) μ as and (37 ± 7) μas. If we compare these values with the core FWHM values obtained from model fits, we obtain a good agreement within the errors. Indeed, for the 86 GHz core, we obtained an average FWHM of 0.023 ± 0.016 mas, while 0.03 ± 0.02 mas is the average size of the 43 GHz core as reported in Jorstad et al. (2017). Given the lower values of the core size at 86 GHz (as inferred from both the fits and the model fits) relative to the corresponding values at 43 GHz, we deduced that in BL Lac the VLBI core at 43 GHz is not completely resolved. To avoid resolution problems, which could artificially flatten the jet collimation profile and mimic a parabolic expansion, Kovalev et al. (2020) considered data starting from 0.9 mas for BL Lac. However, the different expanding rate region observed between ∼0.5 mas and 1.5 mas (Fig. 3) could introduce a bias in the fit if it is not excluded from the analysis. Considering data from 0.9 mas without excluding the aforementioned region would flatten the jet expansion profile beyond 1.5 mas, producing a parabolic geometry instead of a conical expansion. We think this could explain the parabolic shape obtained by Kovalev et al. (2020) in the innermost 2.5 mas of the jet of BL Lac as well as the disagreement between the 15 GHz core size they obtained by fitting the radial width and the visibility data (model-fit components), the one they obtained from the fit being more than two times that obtained from MOJAVE model fits.

We used the BH mass as estimated in Woo & Urry (2002), MBH ∼ 1.6 × 108M, to convert the radial de-projected distances into physical distances (Schwarzschild radii, Rs). The VLBI core, which we assume is co-spatial at the two frequencies, is located at ∼0.59 × 105RS from the BH. In Fig. 5, we superimpose the jet collimation profiles in de-projected distances, as obtained at the two frequencies. The 3 mm data, stacked together to obtain a higher dynamic range image, can recover the same jet width obtained from stacking 7 mm images. They also provide the resolution needed to shed light on the innermost structure of the jet in BL Lac. The higher resolution at 86 GHz helps to better constrain the onset of the higher expanding rate region around 5 pc (∼3.3 × 105RS). Later, at ∼17 pc (∼2.2 × 106RS), the jet recollimates and begins to expand once again, following the initial expansion profile.

thumbnail Fig. 5.

Overlap of the jet collimation profiles in BL Lac as obtained from 43 GHz VLBA and 86 GHz GMVA observations. The 86 GHz data define the starting point of the higher expanding rate region better.

3.2. The inner jet structure and kinematics

In the previous section, we reported the existence of a region in the jet collimation profiles where the jet expands more rapidly, which we excluded from our analysis of the overall jet geometry. We focus our attention here on that region, and we try to contextualize it within the jet kinematics and the physical conditions of the external medium. The kinematics in the first 0.5 mas is rather complex, even when analyzed with the higher resolution provided by the 86 GHz data. Aside from the core, which we consider stationary, we modeled this first part of the jet with four main stationary components (Figs. 6 and 7). In Rani et al. (2016), a similar source modeling for BL Lac has been presented for a GMVA data set obtained in May 2013, hence within our observing period.

thumbnail Fig. 6.

Total intensity images (black contours) of BL Lac obtained with 86 GHz GMVA observations. Contours are drawn at 0.77, 1.53, 3.02, 5.96, 11.74, 23.15, 45.65, and 90% of the peak (1.528 Jy beam−1) registered in May 2011. The restoring beams, shown at the bottom-left corner, range from 0.04 × 0.1 mas to 0.08 × 0.22 mas in FWHM. Model-fit components (blue circles) are overlaid on the maps. In May 2016 and March 2017, an additional component appears upstream of the core.

thumbnail Fig. 7.

Temporal evolution of the distances from the core of the 86 GHz model-fit components. Symbol sizes are proportional to flux density.

A quasi-stationary feature, located at 0.26 mas from the 15 GHz VLBI core, has been reported within the MOJAVE survey and interpreted as a recollimation shock (Cohen et al. 2014). If we assume a core shift of 0.1 mas between the 15 GHz core and the 43 GHz core (O’Sullivan & Gabuzda 2009), the stationary feature should be located around 0.36 mas in our 86 GHz images. Indeed, the 0.3–0.4 mas region in the 86 GHz jet is modeled with two components: d3, located at ∼0.3 mas, and d4, which shifts in position between 0.35 and 0.4 mas. The two components are never observed at the same time, with the exception of the May 2009 epoch, though the identification of d3 at that time is in doubt. Similar behavior is observed for the emission around 0.1–0.2 mas. It is hard to say if components d1 − d2 as well as components d3 − d4 should be considered separately or if they instead describe the same emitting feature, which shifts in position over time due to core-shift variability. The core-shift variability is strongly related to the flux density variability of the core (Plavin et al. 2019), and the 86 GHz core is observed in a rather variable state (Fig. 8). The same radio core displacement effect was adopted by Arshakian et al. (2020) to explain the motions of the 15 GHz quasi-stationary feature.

thumbnail Fig. 8.

Light curves of the 86 GHz GMVA model-fit components.

The location of the stationary emitting region around 0.35 mas coincides with a decrease in the jet width (recollimation), as shown in Fig. 4. This would support the nature of the stationary feature as a recollimation shock. Moreover, the fact that the jet, starting a small distance downstream, expands at a different rate is a hint of change in the physical conditions of either the jet or the external medium (or both).

While the innermost 0.5 mas of the jet is better explored by the higher resolution provided by the 86 GHz data, the jet geometry and brightness distribution until ∼2.9 mas (31 pc, de-projected) are provided by the 43 GHz data set. As visible from Fig. 9, the rapid expansion that the jet undergoes from ∼0.4 mas to 1 mas is associated with a sharp decrease in total intensity. The flux density goes from ∼100 mJy beam−1 at 0.5 mas to 4.5 mJy beam−1 at 1 mas, and then it increases again until 1.5 mas, at which point the flux density is back at ∼20 mJy beam−1.

thumbnail Fig. 9.

Jet width (top) and flux density (bottom) versus distance along the ridge line at 43 GHz. The jet width and flux density are computed for each transverse slice, as explained in Sect. 2.

Bach et al. (2006) suggested that the low brightness region – visible between ∼0.7 mas and 1 mas in the 15, 22, and 43 GHz VLBI images – is the result of reduced Doppler boosting associated with either a bending of the jet, which becomes misaligned with our line of sight, or a change in the jet speed. We instead find that the region of lower brightness usually observed in high-frequency VLBI images of BL Lac is naturally explained by the rapid expansion of the jet and the associated adiabatic expansion losses and decrease in magnetic field strength. This is also supported by the fact that the low brightness region has a dependence on frequency, being less pronounced at low frequencies (e.g., Pearson & Readhead 1988; Bach et al. 2006). When adiabatic losses are the dominant energy loss mechanism, a shift of the turnover frequency toward lower frequencies is indeed expected (Marscher & Gear 1985).

Between 1 mas and 1.5 mas, the jet undergoes a compression: The jet width decreases and the flux density increases. The flux density increase could be explained by a local increase in particle density at the site where the jet recollimates or by reacceleration of relativistic electrons (Mizuno et al. 2015; Nishikawa et al. 2020).

3.3. Jet expansion and the Bondi radius

In a conical jet, as is the case for BL Lac, the jet freely expands in the external medium, and the magnetic field and the electrons’ Lorentz factor decrease with the expansion. To have a change of the status quo, a pressure mismatch between the jet and the ambient medium should occur (see, for example, Fig. 26 in Fromm et al. 2013 for a comparison between a conical and an over-pressured jet). A pressure imbalance naturally arises if the ambient pressure decreases, and, as a result, a shockwave forms to reestablish the equilibrium between the jet and the external medium (e.g., Gómez et al. 1995; Komissarov & Falle 1997; Agudo et al. 2001).

In order to explain the wider jet expansion and further recollimation observed in BL Lac, we have to assume a change in the external pressure, most probably occurring at the location where the jet changes its expansion rate (∼5 pc, 3.3 × 105RS). The Bondi radius, which delimits the sphere of gravitational influence of the BH and usually extends to a distance of 105–106RS (Asada & Nakamura 2012, Boccardi et al. 2020), could be responsible for this apparent change in the surrounding pressure of the environment into which the jet propagates. Inside the Bondi sphere, the external pressure is expected to decrease with distance from the BH. Beyond the Bondi radius, the pressure profile instead becomes flatter, corresponding to the host galaxy’s general pressure profile. The jet’s inertia causes it to overexpand so that it becomes under-pressured relative to the external medium, producing the condition for a recollimation shock.

The Bondi radius, rB, is obtained by setting the escape velocity equal to the sound speed, cs (rB = ). Given the dependence of sound speed on temperature, cs ∼ 104T1/2 cm s−1, the Bondi radius can be expressed as follows (Russell et al. 2015):

(1)

Since we do not have information about the temperature of the accretion flow, but we do have clues regarding the location of the Bondi radius, we can invert the above formula to obtain a range of temperatures, associated with possible locations of the Bondi radius, to compare with values reported in the literature for other sources (e.g., Allen et al. 2006; Balmaverde et al. 2008). Nevertheless, we are aware that the resulting temperatures are sensitive to the assumed BH mass estimate (MBH ∼ 1.6 × 108M; Woo & Urry 2002).

In the radio galaxy M 87, it has been proposed that a change in the external pressure profile associated with the Bondi sphere is responsible for the jet recollimation occurring at the location of the stationary jet feature HST-1 and for the jet geometry transition (from parabolic to conical) observed at a similar distance (Asada & Nakamura 2012). A possible location for the Bondi radius in BL Lac is the millimeter-VLBI core, if it is also associated with a recollimation shock (Gómez et al. 2016; Marscher et al. 2008), or even upstream, at the end of the ACZ (where the jet should go from a parabolic to a conical shape), which Marscher et al. (2008) located 0.65 pc upstream from the millimeter-VLBI core. Therefore, if we place rB at the site of the millimeter-VLBI core (∼0.9 pc), the electron temperature would be 5.45 keV, and even higher if we instead consider the end of the ACZ. The obtained temperature, in this case, is too high compared to values found in the literature. If rB is instead located at around 5 pc (where the jet starts expanding faster) and 10 pc (where the jet experiences the maximum width and the minimum flux density), then the electron temperature in the accretion material would be between 0.99 and 0.33 keV. Balmaverde et al. (2008) inferred the temperature at the accretion radius for 44 low power radio galaxies with values ranging between 0.39 ± 0.10 keV and 1.5 ± 0.1 keV. Since BL Lacs are also low power sources, considered the beamed subset of low power radio galaxies, it is natural to suppose that the temperature of the accretion flow in BL Lac objects should not deviate much from the values reported in Balmaverde et al. (2008). We therefore conclude that the Bondi sphere should extend much farther than the millimeter-VLBI core and that the change we observe in the jet expansion, around 5 pc de-projected (3.3 × 105RS), could be caused by a change in the external pressure profile in the proximity of the Bondi radius.

3.4. The recollimation shock at ∼1.5 mas

The change in the ambient pressure profile, which the jet experiences crossing the Bondi sphere, causes a pressure mismatch between the jet and the external medium. In these circumstances, numerical hydrodynamical and magnetohydrodynamic simulations predict the formation of a recollimation shock (e.g., Gómez et al. 1995; Nishikawa et al. 2020) and the subsequent acceleration of the bulk flow (Mizuno et al. 2015). The recollimation shock should reduce the Lorentz factor while compressing the plasma. As the plasma flows downstream from the shock, the Lorentz factor increases (Gómez et al. 1995; Mizuno et al. 2015).

We investigated the kinematics at 43 GHz, but the identification of components in the region around 1.5 mas is complicated, mainly due to the low brightness region before it (discussed previously), which makes it difficult to track the same component in time, as already reported in Jorstad et al. (2017). Hence, we scrutinized the kinematics at 15 GHz provided by the MOJAVE program because there the low brightness region is less pronounced. We found that six components (components 4, 5, 6, 11, 16, and 20)6 out of the 15 components identified starting from at least 1 mas and moving until at least 2 mas exhibit a clear acceleration parallel to the velocity vector after ∼1 mas. This may support the presence of a recollimation shock around 1.5 mas, as already suggested by the jet collimation profile (Figs. 3 and 5).

3.5. Global jet shape

In Fig. 10 we compare the overall jet shape obtained from our 43 and 86 GHz data sets with the one obtained using 15 GHz MOJAVE data. We downloaded the publicly available stacked image of BL Lac, which is also the one presented in Pushkarev et al. (2017), and we obtained the jet collimation profile in the same way as we did at the other two frequencies.

thumbnail Fig. 10.

Global jet profile in BL Lac as obtained from 86 GHz (GMVA), 43 GHz (VLBA-BU-BLAZAR), and 15 GHz (VLBA MOJAVE) data. The overall jet shape is conical and in agreement at the three frequencies.

We find agreement in the global jet profiles at the three frequencies. The conical shape is also well reproduced by the 15 GHz data, as already found in Pushkarev et al. (2017). The higher expansion rate region seen at 43 GHz is also present at 15 GHz (as is a hint of recollimation), although the region is less pronounced than that observed at 43 GHz.

4. Summary

We studied the jet collimation profile in BL Lac using 3mm GMVA observations from May 2009 to April 2017 as well as 7mm VLBA data from the VLBA-BU-BLAZAR program obtained between May 2009 and March 2019. We built the stacked images at the two frequencies, and we analyzed the jet width at different separations from the VLBI core. Knowing the BH mass in BL Lac (Woo & Urry 2002), the BH–VLBI core distance, and the viewing angle (Marscher et al. 2008), we also translated the radial distance in de-projected Schwarzschild radii and parsecs from the BH. Our results can be summarized as follows.

  • The jet expands with an overall conical geometry, which is also confirmed by the analysis of 15 GHz MOJAVE data.

  • The overall jet expansion is interrupted by a different expanding region between ∼5 pc and 17 pc, where the jet expands faster and then recollimates. The higher expansion rate is accompanied by a severe decrease in brightness, which is explained by the adiabatic cooling being the main energy loss mechanism there. The jet recollimation, which occurs around 17 pc downstream, is instead associated with an increase in flux density, which is explained either by particle re-acceleration or by a local increase in particle density.

  • The change in the jet expansion profile could be associated with a change in the external pressure profile occurring where the BH ceases its gravitational influence on the surrounding material, at the Bondi radius (∼3.3 × 105RS).

  • Both the kinematics and the collimation profile at 3 mm support the existence of a stationary recollimation shock just before the jet starts expanding faster. The location coincides with the recollimation shock observed at 15 GHz within the MOJAVE survey (Cohen et al. 2014, 2015; Arshakian et al. 2020).

  • The parallel acceleration displayed by some model-fit components at 15 GHz (MOJAVE survey) after ∼1 mas (11 pc) is in agreement with the jet physical conditions downstream from a recollimation shock, as predicted by numerical simulations (e.g., Mizuno et al. 2015). This could prove the existence of a recollimation shock around 17 pc, as already observed in the jet expansion profile and supported by the radial flux density profile.

In this scenario, the Bondi sphere, and hence the external medium, play the main role in shaping the jet of BL Lac, which is already freely expanding in the region explored by millimeter-VLBI data. We think that the jet acceleration region predicted by theoretical models (Vlahakis & Königl 2004) ends before the millimeter core, as proposed by Marscher et al. (2008). There, the large-scale helical magnetic field observed by Gómez et al. (2016) could be responsible for the jet collimation.

The case of BL Lac resembles the situation in the radio galaxy M 87, where the stationary feature HST-1 is observed in the proximity of the Bondi radius (Asada & Nakamura 2012). However, in BL Lac we do not observe a change between a parabolic and a conical expansion, as observed in M 87.


Acknowledgments

We thank Eduardo Ros for the valuable suggestions that helped improve this manuscript. This research has been consulting data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2009). C. Casadio acknowledges support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation programme under the grant agreement No. 771282. The research at Boston University was supported by NASA through Fermi Guest Investigator program grants NNX17K0649 an 80NSSC20K1567. This research has made use of data obtained with the Global Millimeter VLBI Array (GMVA), which consists of telescopes operated by the MPIfR, IRAM, Onsala, Metsahovi, Yebes, the Korean VLBI Network, the Greenland Telescope, the Green Bank Observatory and the Very Long Baseline Array (VLBA). The VLBA is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The data were correlated at the correlator of the MPIfR in Bonn, Germany. We thank W. Alef, A. Bertarini, H. Rottmann and I. Wagner for their support at the MPIfR VLBI correlator. We also thank Pablo Torne for his help at the IRAM 30m telescope.

References

  1. Agudo, I., Gómez, J.-L., Martí, J.-M., et al. 2001, ApJ, 549, L183 [NASA ADS] [CrossRef] [Google Scholar]
  2. Akiyama, K., Asada, K., Fish, V. L., et al. 2018, Galaxies, 6, 15 [NASA ADS] [CrossRef] [Google Scholar]
  3. Allen, S. W., Dunn, R. J. H., Fabian, A. C., Taylor, G. B., & Reynolds, C. S. 2006, MNRAS, 372, 21 [NASA ADS] [CrossRef] [Google Scholar]
  4. Arshakian, T. G., Pushkarev, A. B., Lister, M. L., & Savolainen, T. 2020, A&A, 640, A62 [EDP Sciences] [Google Scholar]
  5. Asada, K., & Nakamura, M. 2012, ApJ, 745, L28 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bach, U., Villata, M., Raiteri, C. M., et al. 2006, A&A, 456, 105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Baczko, A. K., Schulz, R., Kadler, M., et al. 2016, A&A, 593, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Baczko, A. K., Schulz, R., Kadler, M., et al. 2019, A&A, 623, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Balmaverde, B., Baldi, R. D., & Capetti, A. 2008, A&A, 486, 119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Beskin, V. S., & Nokhrina, E. E. 2006, MNRAS, 367, 375 [NASA ADS] [CrossRef] [Google Scholar]
  11. Blandford, R. D., & Königl, A. 1979, ApJ, 232, 34 [NASA ADS] [CrossRef] [Google Scholar]
  12. Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] [CrossRef] [Google Scholar]
  13. Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [Google Scholar]
  14. Boccardi, B., Krichbaum, T. P., Bach, U., Bremer, M., & Zensus, J. A. 2016a, A&A, 588, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Boccardi, B., Krichbaum, T. P., Bach, U., et al. 2016b, A&A, 585, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Boccardi, B., Krichbaum, T. P., Ros, E., & Zensus, J. A. 2017, A&ARv, 25, 4 [NASA ADS] [CrossRef] [Google Scholar]
  17. Boccardi, B., Perucho, M., Casadio, C., et al. 2020, A&A, 647, A67 [Google Scholar]
  18. Casadio, C., Marscher, A. P., Jorstad, S. G., et al. 2019, A&A, 622, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Cohen, M. H., Meier, D. L., Arshakian, T. G., et al. 2014, ApJ, 787, 151 [NASA ADS] [CrossRef] [Google Scholar]
  20. Cohen, M. H., Meier, D. L., Arshakian, T. G., et al. 2015, ApJ, 803, 3 [NASA ADS] [CrossRef] [Google Scholar]
  21. Fromm, C. M., Ros, E., Perucho, M., et al. 2013, A&A, 557, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Giroletti, M., Giovannini, G., Feretti, L., et al. 2004, ApJ, 600, 127 [NASA ADS] [CrossRef] [Google Scholar]
  23. Gómez, J. L., Marti, J. M. A., Marscher, A. P., Ibanez, J. M. A., & Marcaide, J. M. 1995, ApJ, 449, L19 [Google Scholar]
  24. Gómez, J. L., Lobanov, A. P., Bruni, G., et al. 2016, ApJ, 817, 96 [Google Scholar]
  25. Greisen, E. W. 1990, Acquisition, Processing and Archiving of Astronomical Images, 125 [Google Scholar]
  26. Hada, K., Kino, M., Doi, A., et al. 2013, ApJ, 775, 70 [NASA ADS] [CrossRef] [Google Scholar]
  27. Hada, K., Doi, A., Wajima, K., et al. 2018, ApJ, 860, 141 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  29. Jorstad, S. G., Marscher, A. P., Morozova, D. A., et al. 2017, ApJ, 846, 98 [NASA ADS] [CrossRef] [Google Scholar]
  30. Kim, J. Y., Krichbaum, T. P., Lu, R. S., et al. 2018, A&A, 616, A188 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Komatsu, E., Dunkley, J., Nolta, M. R., et al. 2009, ApJS, 180, 330 [NASA ADS] [CrossRef] [Google Scholar]
  32. Komissarov, S. S., & Falle, S. A. E. G. 1997, MNRAS, 288, 833 [NASA ADS] [CrossRef] [Google Scholar]
  33. Kovalev, Y. Y., Pushkarev, A. B., Nokhrina, E. E., et al. 2020, MNRAS, 495, 3576 [CrossRef] [Google Scholar]
  34. Lister, M. L., Cohen, M. H., Homan, D. C., et al. 2009, AJ, 138, 1874 [Google Scholar]
  35. Lyubarsky, Y. 2009, ApJ, 698, 1570 [NASA ADS] [CrossRef] [Google Scholar]
  36. Marscher, A. P., & Gear, W. K. 1985, ApJ, 298, 114 [NASA ADS] [CrossRef] [Google Scholar]
  37. Marscher, A. P., Jorstad, S. G., D’Arcangelo, F. D., et al. 2008, Nature, 452, 966 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  38. McKinney, J. C. 2006, MNRAS, 368, 1561 [NASA ADS] [CrossRef] [Google Scholar]
  39. Mizuno, Y., Gómez, J. L., Nishikawa, K.-I., et al. 2015, ApJ, 809, 38 [NASA ADS] [CrossRef] [Google Scholar]
  40. Nakahara, S., Doi, A., Murata, Y., et al. 2018, ApJ, 854, 148 [NASA ADS] [CrossRef] [Google Scholar]
  41. Nakahara, S., Doi, A., Murata, Y., et al. 2020, AJ, 159, 14 [Google Scholar]
  42. Nishikawa, K., Mizuno, Y., Gómez, J. L., et al. 2020, MNRAS, 493, 2652 [CrossRef] [Google Scholar]
  43. O’Sullivan, S. P., & Gabuzda, D. C. 2009, MNRAS, 400, 26 [NASA ADS] [CrossRef] [Google Scholar]
  44. Pearson, T. J., & Readhead, A. C. S. 1988, ApJ, 328, 114 [NASA ADS] [CrossRef] [Google Scholar]
  45. Plavin, A. V., Kovalev, Y. Y., Pushkarev, A. B., & Lobanov, A. P. 2019, MNRAS, 485, 1822 [NASA ADS] [CrossRef] [Google Scholar]
  46. Potter, W. J., & Cotter, G. 2015, MNRAS, 453, 4070 [Google Scholar]
  47. Pushkarev, A. B., Kovalev, Y. Y., Lister, M. L., & Savolainen, T. 2017, MNRAS, 468, 4992 [Google Scholar]
  48. Rani, B., Krichbaum, T., Hodgson, J., et al. 2016, Galaxies, 4, 32 [Google Scholar]
  49. Russell, H. R., Fabian, A. C., McNamara, B. R., & Broderick, A. E. 2015, MNRAS, 451, 588 [NASA ADS] [CrossRef] [Google Scholar]
  50. Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2008, MNRAS, 388, 551 [NASA ADS] [CrossRef] [Google Scholar]
  51. Traianou, E., Krichbaum, T. P., Boccardi, B., et al. 2020, A&A, 634, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [CrossRef] [Google Scholar]
  53. Van Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: CreateSpace) [Google Scholar]
  54. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  55. Vlahakis, N., & Königl, A. 2004, ApJ, 605, 656 [Google Scholar]
  56. Woo, J.-H., & Urry, C. M. 2002, ApJ, 579, 530 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: The jet collimation profile in BL Lac using AIPS

Figures A.1 and A.2 show the jet collimation profiles at 43 and 86 GHz, respectively, obtained using the tasks SLICE and SLFIT in AIPS. Both profiles as well as the fit parameters, obtained using the least-squares method, are in agreement with that obtained using the ridge line method (see Sect. 3 and Figs. 3 and 4). The difference here is that we could not fit the w − d dependence after 1.5 mas in the 43 GHz profile because of some missing values around 2 mas, where AIPS was not able to properly fit the brightness profile associated with some slices.

thumbnail Fig. 11.

Jet width versus distance along the jet in the 43 GHz VLBA stacked image of BL Lac, obtained using AIPS.

thumbnail Fig. 12.

Jet width versus distance along the jet in the 86 GHz VLBA stacked image of BL Lac, obtained using AIPS.

All Tables

Table 1.

Model-fit parameters for the 86 GHz GMVA observations.

Table 2.

Best-fit parameters obtained for the power law equations w = a + b ⋅ dk and w = b ⋅ dk fitting the jet width–distance dependence at 43 and 86 GHz.

All Figures

thumbnail Fig. 1.

Total intensity stacked image of BL Lac at 43 GHz. A total of 94 images obtained between May 2009 and March 2019 are used. The restoring circular beam (lower left circle) is 0.19 mas. Contours are drawn at 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, and 64% of the peak (1.71 Jy beam−1). The total intensity ridge line is shown in red.

In the text
thumbnail Fig. 2.

Total intensity stacked image of BL Lac at 86 GHz. Eleven images obtained between May 2009 and April 2017 are used. The restoring circular beam (lower left circle) is 0.10 mas. Contours are drawn at 0.35, 0.5, 1, 2, 4, 8, 16, 32, and 64% of the peak (0.847 Jy beam−1). The total intensity ridge line is shown in red.

In the text
thumbnail Fig. 3.

Jet width versus radial distance at 43 GHz, calculated with respect to the VLBI core (in mas) and to the BH (in de-projected parsecs). The blue and orange lines are the best fit of assumed power laws of the type w = a + b ⋅ dk and w = b ⋅ dk, respectively. The best-fit parameters with their uncertainties are also reported in Table 2.

In the text
thumbnail Fig. 4.

Jet width versus radial distance at 86 GHz, calculated with respect to the VLBI core (in mas) and to the BH (in de-projected parsecs). The blue line is the same as in Fig. 3.

In the text
thumbnail Fig. 5.

Overlap of the jet collimation profiles in BL Lac as obtained from 43 GHz VLBA and 86 GHz GMVA observations. The 86 GHz data define the starting point of the higher expanding rate region better.

In the text
thumbnail Fig. 6.

Total intensity images (black contours) of BL Lac obtained with 86 GHz GMVA observations. Contours are drawn at 0.77, 1.53, 3.02, 5.96, 11.74, 23.15, 45.65, and 90% of the peak (1.528 Jy beam−1) registered in May 2011. The restoring beams, shown at the bottom-left corner, range from 0.04 × 0.1 mas to 0.08 × 0.22 mas in FWHM. Model-fit components (blue circles) are overlaid on the maps. In May 2016 and March 2017, an additional component appears upstream of the core.

In the text
thumbnail Fig. 7.

Temporal evolution of the distances from the core of the 86 GHz model-fit components. Symbol sizes are proportional to flux density.

In the text
thumbnail Fig. 8.

Light curves of the 86 GHz GMVA model-fit components.

In the text
thumbnail Fig. 9.

Jet width (top) and flux density (bottom) versus distance along the ridge line at 43 GHz. The jet width and flux density are computed for each transverse slice, as explained in Sect. 2.

In the text
thumbnail Fig. 10.

Global jet profile in BL Lac as obtained from 86 GHz (GMVA), 43 GHz (VLBA-BU-BLAZAR), and 15 GHz (VLBA MOJAVE) data. The overall jet shape is conical and in agreement at the three frequencies.

In the text
thumbnail Fig. 11.

Jet width versus distance along the jet in the 43 GHz VLBA stacked image of BL Lac, obtained using AIPS.

In the text
thumbnail Fig. 12.

Jet width versus distance along the jet in the 86 GHz VLBA stacked image of BL Lac, obtained using AIPS.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.