© The Authors 2023

1 Introduction

Collimated jets and outflows from low-mass young stars and protostars are ubiquitous and are now believed to be intimately linked with accretion discs. The most widely accepted paradigm is that they originate from magnetically collimated disc winds and help remove angular momentum from the disc, allowing the gas to accrete onto the star surface (see, e. g. Ray & Ferreira 2021 and references therein). Another well-established observational fact in low-mass star formation is the so-called Luminosity Problem, whereby protostars are usually much fainter than predicted for expected accretion rates of ~10⁻⁴–10⁻⁵ M_⊙ yr⁻¹ (see, e.g. McKee & Offner 2011), pointing to a much lower accretion activity. A proposed scenario to circumvent the problem is episodic accretion, namely protostars grow through short episodes of intense matter inflows. If accretion and outflows are linked together, this would in fact agree with the morphology of many optical and near-infrared (NIR) jets that appear composed of aligned knots (emitting in several optical and NIR spectral lines), which, in this view, would be relics of past accretion bursts. In this scenario, both protostars and very young stars should also exhibit photometric variability with augmented luminosity occurring during phases of intense accretion, which has actually been observed in both low-mass (see, e. g. Yoo et al. 2017; Caratti o Garatti & Eislöffel 2019), intermediate-mass (Benisty et al. 2010) and high-mass young stars (see, e. g. Caratti o Garatti et al. 2017; Hunter et al. 2017).

In high-mass star formation, the picture is observationally more problematic as high-mass protostars are rare and usually far away (at kilo-parsec distances), and embedded in crowded regions where low-mass protostars are also present. This highlights the need for high-spatial-resolution observations. In addition, studying photometric variability of a few selected high-mass-star precursors and their associated outflows can be a valuable tool to understand how accretion occurs in these objects, circumventing the need to resolve the driving protostars. This is of course better achieved when multi-epoch high-spatial-resolution observations are available.

In this respect, IRAS 20126+4104 is a source of great interest. This massive protostar has a luminosity of ~10⁴ L_⊙ (Hofner et al. 2007; Johnston et al. 2011) and is surrounded by a disc; this was first detected by Cesaroni et al. (1997) and confirmed by subsequent higher-resolution studies (Zhang et al. 1998; Cesaroni et al. 1999, 2005, 2014). It appears to be powering a parsec-scale jet-outflow undergoing precession about the disc axis (Shepherd et al. 2000; Cesaroni et al. 2005; Caratti o Garatti et al. 2008). This jet-outflow has been imaged on scales from a few 100 au (see Fig. 1; Moscadelli et al. 2005; Cesaroni et al. 2013) to 0.5 pc (Shepherd et al. 2000; Lebrón et al. 2006) and its 3D expansion velocity, close to the source, has been measured through maser proper motions (Moscadelli et al. 2005). Two radio sources (see Fig. 1) have been detected at high-spatial resolution and are associated with the protostar (Hofner et al. 1999, 2007). The brighter one, N1, has been detected at 3.6, 1.3, and 0.7 cm and exhibits a spectral index consistent with thermal emission. It is likely to be located ~0″.3 north-west of the protostar, whose position has been inferred by OH maser emission at 1665 MHz, displaying an elongated structure with a north-eastern-south-western velocity gradient reminiscent of Keplerian rotation (Edris et al. 2005). The fainter one, N2, has only been detected at 3.6 cm, is located ~1″ north-west of N1 and its morphology is reminiscent of a bow shock. Radio sources N1 and N2 are roughly aligned in the direction of the jet, are associated with the water maser system, and have been interpreted as a manifestation of gas that has been shocked by the jet from the protostar (Hofner et al. 2007). Another feature found is a faint radio source, S (see Fig. 1, top panel), which has only been detected at 3.6 cm, and lies ~1″ south of N1 and is roughly elongated in a direction parallel to that of the jet. This has been interpreted as a distinct radio jet and its contribution to the outflow should be negligible (Shepherd et al. 2000; Hofner et al. 2007); Sridharan et al. (2005, 2007) found a point-like source in L- and M-band high-resolution images coinciding with S, proving that the latter has a stellar counterpart. A clear X-ray detection was obtained with Chandra and the X-ray spectrum is consistent with an embedded early B-type star (Anderson et al. 2011).

The distance to IRAS 20126+4104 (1.64±0.05 kpc) has been accurately determined from parallax measurements of water masers (Moscadelli et al. 2011), which prove that this is one of the closest discs around a B-type protostar and this allows one to observe its environment with excellent linear resolution. More recently, Nagayama et al. (2015) repeated the parallax measurement with VERA, obtaining a distance of ${1.33}_{-0.092}^{+0.19}$ kpc. As the difference with the value of Moscadelli et al. (2011) is not statistically significant, in the following, we continue to use a distance of 1.64±0.05 kpc.

Sub-arcsecond imaging (Cesaroni et al. 2014) and model fitting (Chen et al. 2016) have proven that the disc is undergoing Keplerian rotation around a ~12 M_⊙ protostar, as suggested by the butterfly-shaped position-velocity plot obtained in almost every hot molecular core tracer observed. Evidence of accretion onto the protostar has also been reported (Cesaroni et al. 1997; Keto & Zhang 2010; Johnston et al. 2011). The jet traced by the water masers over a few 100 au from the protostar appears to be co-rotating with the disc (Trinidad et al. 2005; Cesaroni et al. 2014).

Imaging at infrared wavelengths (Qiu et al. 2008) suggested that IRAS 20126+4104 was associated with an anomalously poor stellar cluster. However, subsequent X-ray observations by Montes et al. (2015) revealed an embedded stellar population that hints at the existence of a richer (and possibly very young) cluster.

Measurements of the magnetic field orientation both at the clump scale (Shinnaga et al. 2012) and at the disc scale (Surcis et al. 2014) show that the direction of the field lines is basically the same from ~0.2 pc to ~1000 au and roughly perpendicular to the outflow axis. This suggests that the accreting material could be flowing onto the star through the disc along the magnetic field lines.

With all the above in mind, we think that IRAS 20126+4104 is a suitable target for high-resolution imaging at NIR wavelengths, in order to measure the physical properties of its protostellar jet, as previously undertaken, in part, by Caratti o Garatti et al. (2008) and Cesaroni et al. (2013). These works detected a well-collimated and structured precessing jet, with two lobes, a south-eastern one (red-shifted) and a north-western one (blue-shifted), composed of several knots and bow-shock structures (see top panel of Fig. 1), with each lobe spanning an arc of ~10″ (~0.08 pc). Larger field images show further distant knots which are not aligned with the lobes, but they are reminiscent of a wiggled configuration. Wiggling is also evident in the two lobes, suggesting that the outflowing direction is in fact precessing with a period whose estimates range from ~1100 yr (for the closest knots; Caratti o Garatti et al. 2008) to ~64 000 yr (for the outer flow; Shepherd et al. 2000; Cesaroni et al. 2005). A puzzling feature that stands out in the images is a small (~1″) ring-like patch (hereby, the ring-like feature or X, see Fig. 1) that mainly emits in the H₂ lines. This is not the circumstellar disc, which appears to be heavily obscured and invisible in the NIR, but a structure located a few tens of arcsecs south-east of the disc, hence facing the red-shifted lobe of the jet (see Fig. 1).

We therefore targeted IRAS 20126+4104 for observations with the second generation adaptive optics (AO) system SOUL that is nearly commissioned at the Large Binocular Telescope (LBT). The complex jet structure and the presence of a suitable guide star are optimal features to test the performances of an AO system. In addition, we aimed to complement the new observations with other high-resolution images available from 2000 on to derive, for the first time, proper motions, 3D kinematics and photometric variability along the NIR jet of a massive protostar.

The paper layout is the following. Section 2 describes the observations and the data reduction, including all the ancillary data. The results of our data analysis are dealt with in Sect. 3. The implications of our results are discussed in Sect. 4.

2 Observations and data reduction

2.1 SOUL AO at the LBT

SOUL (Pinna et al. 2016, 2021) is the upgrade of the FLAO systems (Quirós-Pacheco et al. 2010) at the LBT with an electro-multiplied camera for the wavefront sensor. As FLAO, SOUL is a single conjugated AO system, using the adaptive secondary mirror as a corrector and one natural guide star for a pyramid wavefront sensor. The upgrade allows us to fully exploit the light flux from the reference star, improving the correction. In this observation, we used the R = 13.49 star (according to UCAC5; Zacharias et al. 2017) as the AO reference, which is located at the western edge of the field (Fig. 1, bottom), correcting 91 modes at 850 Hz.

2.2 LBT/SOUL NIR observations

The new AO-assisted observations were carried out with the NIR imager LUCI1 at the beginning of November 2020, during the commissioning time of SOUL. On-source and off-source dithered images were obtained through the narrow-band filters H2 and Brγ (centred at 2.124 µm and 2.170 µm, respectively), and the broad-band filter K_s (centred at 2.163 µm) with the N30 camera (plate scale 0″.015 pix⁻¹, field of view 30″ × 30″). LUCI1 is equipped with a HAWAII-2RG HgCdTe detector with 2048 × 2048 pixels. The K_s images were obtained using the standard double-correlated sample readout scheme (indicated as LIR in the instrument handbook), whereas the H2 and Brγ images used a readout scheme called MER, consisting of a set of five readings at the beginning and at the end of each integration. In addition, each set of NDIT¹ exposures through the H2 and Brγ filters was saved as a 3D FITS file with NDIT planes. A log of the observations is given in Table 1.

The target field was imaged through a sequence of dithered exposures lasting DITxNDIT each, alternating an on-source integration and an off-source one. The centre of field of each on-source frame was randomly shifted with offsets of <3″ in the x and y directions of the detector reference frame. A dark area ~1′ away was selected as an off-source position and a dithering scheme with random offsets <1′ was adopted. A standard star (FS149, Hawarden et al. 2001) was imaged in K_s with an open loop (i.e. with the AO offline) when observations were carried out in this band (see Table 1) using a five-position dithering scheme with the star alternatively in the frame centre and in each of the four quadrants (x and y offsets of 5″). On November 6, only two shifted images of FS149 were taken due to poor weather conditions.

Data reduction was performed using standard IRAF routines. First, the frames were flat-field corrected using dome flat images. All images in each off-source 3D FITS file were averaged together in a single frame. Then, for each on-source 3D FITS file, a sky image was obtained as the median of the four nearest (in time) NDIT-averaged off-source images. For each on-source 3D FITS file, the single exposures were extracted, sky-subtracted, corrected for bad pixels, and checked for any change in point spread function (PSF) width and in flux. Deviant frames were discarded and the remaining ones were registered and averaged together. This was to compensate for any internal flexure from the LUCI instrument during the observations. The final averaged, sky-subtracted on-source images obtained from each 3D FITS file were registered and summed together. The same scheme was used for the K_s frames, except that in this case the NDIT exposures were not saved as 3D FITS files. In the end, the final H2 image was the sum of 65 exposure (total integration time DIT × 65 = 1950 s) with a PSF width ⪝0.1″ (mostly less than 0.083″; measured with IRAF, for an accurate determination of the PSF widths see Appendix A), and the final Brγ image is the sum of 40 exposures (total integration time DIT × 40 = 1200 s) with a PSF width ⪝0.09″. The final K_s images are the sum of five DIT×NDIT exposures (total integration time 300 s, final PSF width <0.11″; for an accurate assessment of the PSF widths see Appendix A) for November 2, of five DIT×NDIT exposures (total integration time 300 s, final PSF width ⪝0.08″) for November 4, of three DIT×NDIT exposures (total integration time 180 s, final PSF width ⪞0.12″) for November 6 (used only for calibration purposes). A three-colour image (red H2 filter, green K_s, blue Brγ filter) of IRAS 20126+4104 is shown in Fig. 1.

To perform photometry of the jet features, we first calibrated the stars in the small field of view (FoV). This was done using the K_s images and the daophot package in IRAF. For each night with K_s observations, we measured the brightness of the standard star FS149 in each position of the frame with different aperture radii. As the instrumental magnitude, we adopted the value at which the measured brightness does not change significantly any more (aperture radius ~1.5″). The zero point was obtained by averaging the instrumental magnitudes obtained at each position in the frame (and the standard deviation was taken as the error on the zero point), assuming - as K_s standard magnitude - the value given in the 2Mass PSC (K_s = 10.052 ± 0.017).

In addition, we must take into account the fact that in the target images the Strehl ratio decreases with the distance to the AO guide star (the bright star west of the target in Fig. 1, i.e. star 8 in Fig. C.1). To minimise this, we analysed the radial profile of the stellar PSF in the frames. This is composed of a narrow peak on top of a broader fainter feature, as expected (see Fig. A.1). The broader feature exhibits a width of ~0″.3 (November 2 and 4) and ~0″.45 (November 6), which was taken as the aperture radius. We constructed the growth curve for the bright star lying to the north of the target (star 4 in Fig. C.1), which is isolated enough, and we derived a correction term for the aperture photometry. Using the zero point estimated from the standard star, we derived the K_s magnitudes of the stars in the target image. As the airmass difference between the target field and standard was always <0.1, no correction for airmass was needed.

Fig. 1 Three colour image (red H2 filter, green KS, blue Brγ filter) of IRAS 20126+4104 obtained with PISCES and FLAO at the LBT in 2012 (top) and LUCI1 and the new AO system SOUL at the LBT in November 2020 (bottom). The KS image used for the bottom image is that taken on November 4, that is the one exhibiting the best spatial resolution. Approximate positions of radio sources N1 (blue square), N2 (yellow square), and S (orange square, see text) are marked in the top panel with small coloured squares.

2.3 Ancillary data

In order to derive the jet proper motion and its photometric variability, we retrieved all previous H2, Brγ and K images, both with low- and high-spatial resolution, from a number of archives. A K image, presented in Sridharan et al. (2005), was obtained with UFTI at UKIRT from observations carried out from August 1518, 2000 (seeing ~0″.3). Narrow-band H2, Brγ, and K_cont images, presented in Cesaroni et al. (2005), were obtained with NICS at the TNG in June 2001 (seeing ~1″.2). Caratti o Garatti et al. (2008) studied a further series of H2, Brγ, and K′ images; namely, a dataset (first presented in Sridharan et al. 2007) acquired at the Subaru telescope with the Coronagraphic Imager with Adaptive Optics (CIAO; Murakawa et al. 2003), used as a high-resolution NIR imager, on July 10,2003 (pixel scale of ~0″.022; seeing ~0″.9; PSF width ~0″.15), and a second dataset obtained with NICS at the TNG on August 6, 2006 (seeing ~1.2″). Finally, high-spatial-resolution H2, Brγ, and K_s images, analysed by Cesaroni et al. (2013), were obtained with PISCES and the AO system FLAO at the LBT on June 21, 2012 (PSF width ~0″.09).

Data reduction is described in the cited papers; in a few cases we redid it to improve the data quality and test photometric stability. Notably, the CIAO, PISCES and SOUL maps roughly cover the same jet region around the source extending ~30″ (see Fig. 1).

Photometry on the large FoV images (UFTI, TNG) was performed with daophot using aperture radii ~1 full width at half maximum (FWHM) of the PSF. Calibration was obtained matching the imaged stars and the 2Mass PSC and performing a linear fit. The fit formal error is typically <0.2 mag (~0.08 for UFTI), but one has to consider that this includes intrinsic variability in a number of stars and the zero point is then expected to be more accurate than ~0.1 mag. As neither the CIAO images nor those obtained with PISCES are associated with observations of standard stars, in these cases we only performed relative photometry on the stars contained in the small FoV, with aperture radii of the order of 1 FWHM of the broad component of the PSF to minimise the effects of a varying Strehl ratio.

2.4 Continuum subtraction

To obtain H₂ and Brγ pure line emission images, we first used the broad-band K_s and K′, and the narrow-band K_cont images, depending on the dataset. These were first rotated and shifted to the corresponding H2 and Brγ images using the IRAF routines geomap and geotran. Then we performed aperture photometry on the common stars (>600 for the large FoV images, 4–15 for the small FoV AO images) in the frames to scale the images by assuming that the intrinsic stellar fluxes do not vary with wavelength, a zeroth-order approximation. We note that this implies that images through different filters do not need to be taken during the same night, as telluric extinction differences are automatically taken into account. Short time stellar variability might affect the subtraction only for the small FoV images, where the number of stars is low. However, the r.m.s. of the average flux ratio is ⪝ 10% in these cases, which indicates no large stellar photometric variations between pairs of images. Finally, the scaled images were subtracted to yield continuum-corrected images. We also note that when the continuum contribution is estimated from K, both images have to be scaled before subtraction to take into account that both include line emission as well as continuum emission.

We did not notice any significant Brγ emission, so in the end we used the available Brγ images to estimate the continuum contribution in the H2 filter. In fact, the broad-band K filters encompass other H₂ emission lines and they are therefore bound to overestimate the continuum contribution. In principle, this overestimate should be ~ 10%; in fact we found larger line fluxes when the continuum contribution was subtracted using a narrow-band filter rather than the K band. However, we also found that the ratio of line flux (measurement with the continuum estimated using a narrow-band filter over the measurement with the continuum estimated using the K band) increases when the ratio of the line to continuum (K) flux in a knot decreases. This indicates that the line flux may be significantly underestimated for faint line emission superimposed on a strong continuum (K) flux when the continuum contribution is evaluated using the K band (see Appendix G).

2.5 Jet photometry

We performed photometry of the jet emission, both in the continuum and in the H₂ 2.12 µm line, using the task polyphot in IRAF. We divided the jet into smaller components and defined a polygon for each one, taking care of encompassing most of its emission (down to a ~5σ level of the background as determined far from the jet). Before doing the photometry, we registered and projected all of the images onto the SOUL H2 frame grid as explained in Appendix B, so we did not have to adapt the polygons to each image. The photometric zero points have been computed for each frame based on two bright stars which are common to all images. The variability of the two stars, the stability of the zero points, and the consistency of the calibration are discussed in Appendices C–E.

As for continuum emission, the selected polygons are described in Appendix F and shown in Fig. 2. Where possible, it was measured on the Brγ-filter images, which do not include the H₂ line emission and where Brγ emission has not been detected. Only in two cases did we use K images (UKIRT 2000 and TNG 2006) after subtracting the derived H₂ line emission (we used TNG 2001 to estimate the correction for UKIRT 2000). For the lower-resolution images (UKIRT and TNG), we also subtracted the stars in the frames by using daophot before performing the photometry. The sources of error involved in the various steps of the photometry are assessed in Appendix F.

Photometry of the jet emission in the H₂ 2.12 µm line was performed in the same way as was done for continuum emission. However, in this case we used two different sets of polygons. First, we defined a set of larger polygons (see Appendix G and Fig. 2.) to exploit both the higher and the lower spatial resolution images. Where possible we used the pure line emission images obtained by subtracting the continuum contribution as estimated from adjacent narrow-band filters. In only one case was the correction estimated from the broad-band K′ filter. As shown in Appendix A, while narrow-band filters (namely Brγ and K_cont) yield a consistent continuum correction, using broadband K images results in more discrepant continuum corrections compared to Brγ, leading to increasingly fainter magnitudes in the line emission with increasing continuum contamination.

Next, we analysed the high-spatial-resolution pure H₂ line emission images from the AO-assisted observations to identify single emission knots and defined a set of smaller polygons encompassing each knot (see Fig. 3 for designation). In this case, we only used the high spatial resolution images (CIAO 2003, PISCES 2012, and SOUL 2020) with the continuum correction derived from the Brγ images. To focus only on the smaller-scale structures, we adapted the multi-scale image decomposition described in Belloche et al. (2011), with seven levels of decomposition, to filter out large-scale emission structures. This should have resulted in the efficient filtering out of ⪞2″ structures and a better background estimate. We note that the shape and location of a few polygons of this second set had to be slightly changed and adapted depending on the image because of the knots’ proper motions.

2.6 Proper motion analysis

We used the H2 continuum-subtracted high-angular resolution images from CIAO, PISCES and SOUL after flux calibration based on the photometric zero points (see Sect. 2.5), and registration to and projection onto the SOUL H2 frame grid (Appendix B) to infer proper motions of knots and bow shocks along the jet. This dataset provides a time baseline of more than 17 yr. We note that the proper motions are essentially derived in the reference frame of the protostar (see Appendix B). To compute the proper motions (PMs), we used the same strategy and our own developed software in python as described in Caratti o Garatti et al. (2009).

Briefly, single shifts were computed between image pairs (namely the SOUL image versus previous PISCES and CIAO images) keeping the latest epoch as a reference and using a cross-correlation method. After identifying each knot or structure in each image, a sub-image, enclosing the structure 5σ contour, was selected and cross-correlated with the corresponding sub-images obtained from the earlier epochs. In a given pair of sub-images, the earlier one was then shifted in steps of 0.1 pixel in x (RA) and y (Dec), and for each shift (∆x, ∆y) a product image was created. The highest-likelihood shift for each epoch is that yielding the maximum correlation signal f (x, y) integrated over the considered structure on the corresponding product image. To quantify the systematic errors for the shift measurements, the size and shape of each sub-image that encloses the considered structure was varied. The resulting range of shifts gives the systematic error, which depends on the structure signal-to-noise ratio (S/N), shape, and time variability. This uncertainty is typically much larger than the register error. Finally, the PM value of each structure was derived from a weighted least square fit of the shifts derived for each epoch, fitting the motion in R.A. and Dec. simultaneously. This fit also provides the associated errors. The results of the fits are detailed in Appendix I.

Fig. 2 Jet decomposition in polygons. Top panel: Brγ-filter image of IRAS 20126+4104 obtained with LUCI1 and the AO system SOUL, overlaid with the polygons used for the photometry in the continuum. The adopted designation is labelled. We note that polygon S is composed of S1 and S2, and S1 is further subdivided into S1a and S1b. Bottom panel: pure H2 line emission image of IRAS 20126+4104 obtained with LUCI1 and the AO system SOUL, overlaid with the polygons used for the jet photometry. The adopted designation is labelled. We note that polygon LLS is composed of LLS1 and LLS2, and LLN2 is subdivided into LLN2a and LLN2b.

Fig. 3 Pure H2 line emission image of IRAS 20126+4104 obtained with LUCI1 and the AO system SOUL, overlaid with the polygons used for the high-spatial-resolution jet photometry.

3 Results

3.1 Morphology of the jet and outflow cavities

The three-colour image (red H2 filter, green K_S, blue Brγ filter) obtained from the SOUL data, shown in the bottom panel of Fig. 1, is quite similar to the PISCES one presented in Fig. 1 of Cesaroni et al. (2013), which has been adapted in the top panel of Fig. 1. The continuum emission (Brγ and K_s filters) delineates the outflow cavities (in blue and cyan), which scatter the radiation from the central source, which is undetected at 2 µm. The H₂ shocked emission traces a precessing jet (in red; see also Caratti o Garatti et al. 2008), which is highly fragmented, especially towards the north-western blue-shifted side, displaying a large number of knots and bow shocks. It is worth pointing out that our high-resolution images show only the inner portion of the flow (~30 ″ in size or ~0.24pc), which extends further north and south (see, e. g. Shepherd et al. 2000, Lebrón et al. 2006).

To study the jet variability and kinematics, each sub-structure (knots and bow shocks) showing local emission peaks above 5σ was identified and labelled in the SOUL H2 continuum-subtracted image. In our kinematical analysis we followed the general nomenclature used by Caratti o Garatti et al. (2008), who divided the flow regions into six different groups: A1, A2, X, B, C1, and C2 (see also Fig. 8d).

3.2 Continuum variability

We found the most remarkable flux variations in the continuum emission. To analyse such variability, we selected seven continuum emitting regions: three in the north-western and southeastern cavity, respectively, and one enclosing the ring-like feature close to the source position (see Appendix F and Fig. 2). The light curves of these regions are displayed in Fig. 4. The figure shows a clear decrease of the continuum emission in the north-western lobe between 2000 and 2006 (see also Fig. 5, upper and bottom left panels, and the short movie in Fig. H.1), with a subsequent increase. Conversely, the south-eastern lobe seems to have brightened up during the same time interval and then to have decreased in intensity. The ring-like feature (labelled PD in Fig. 4) appears to exhibit a slow steady flux decrease.

3.3 H₂ line variability

The line photometry on the set of larger polygons (including the low spatial resolution images; see Appendix F and Fig. 2) is shown in Fig. 6. We note that, as discussed in Appendix G, the fluxes of the data points obtained by correcting the continuum contamination using the K′ image are expected to be underestimated, so the dip displayed in the plot on JD 2453953 (black solid squares) is likely not real. A very rough correction of 0.3 mag (open triangles in the figure), which should represent an upper limit, has been obtained following Appendix G. Furthermore, the obtained magnitudes need to be corrected for the different band-widths of the H2 filters (see Appendix D). We adopted the LUCI1 filter bandwidth as a standard and corrected all other measurements using Eqs. (D.10) and (D.11). The line flux can be derived from the magnitude m given in the figure through the following: $$F=1.06\times {10}^{-15}\times {10}^{-0.4m}\text{W\hspace{0.17em}}{\text{cm}}^{-2}.$$(1)

As we cannot rule out systematic errors in the zero points of ~0.1–0.2 mag in an unknown direction, the maximum variation intervals of ⪝0.4 mag per polygon indicate that the total energy emitted in the 2.12 µm line is probably almost constant in time within a few tenths of magnitude in the sampled regions.

The knot line variability on the smaller scale (excluding the lower spatial resolution images; see Fig. 3 for polygon designations) is shown in Fig. 7. Most knots exhibit a constant flux within ~0.2–0.4 mag from 2003 to 2020 in a similar trend as found from the low spatial resolution photometry. Nevertheless, even assuming zero point systematic errors of ~0.2 mag in opposite directions, variations >0.4 mag should be considered as real. In this respect, an emission increment ≥0.8mag is evident for knot HLPD2 in the ring-like feature, the outermost north-western knot (HLN25), knot HLN11 adjacent to star 6, knot HLN13 connected with bow-shock C1h and C1k (see Sect. 4.2), and knot HLN18 in the middle section of C1 (see Fig. 8). Another five knots exhibit an increase of ~ 0. 5 mag, namely HLN2, HLN14, and HLN19 and, in the south, HLS1 and HLS3 (see Fig. 7). On the other hand, two knots exhibit a decrease of at last ~0.5 mag, that is HLN16c and HLN23c, and they are located in outer areas of C1.

An increase in line emission implies an increase in the column density of molecules in the upper level of the transition integrated over the knot area. In bow shocks, this column density increases with increasing shock speed in the range ~5–15 km s⁻¹ (Tram et al. 2018). A generalised trend of brightening therefore suggests that more and more gas is being entrained into the shocked regions.

Due to the shape of the PSF in the AO-assisted images, the filtering is bound to cause some of the smaller-scale flux to be missed since a residual small fraction of the PSF flux is spread on a larger size compared to the PSF diffraction-limited peak. Thus, it is essential to ensure that this does not affect the time variations displayed by the single knots as well. By redoing the same photometry on the unfiltered images, we found the same trends as from the filtered images, with the main difference being that the knots are a few tenths of magnitude fainter in the filtered images (down to ~ 1 mag fainter in only a few cases). So we can confirm that the flux from most knots is stable within a few tenths of magnitude, with some exceptions.

Fig. 4 Photometric variability of the continuum features defined in Fig. 2. We note that N and S indicate the northern and southern lobe, respectively, and PD is the ring-like feature. The numbers increase from south-east to north-west. We note that S is further subdivided into S1 and S2, and S1 is in turn subdivided into S1a and S1b. The symbol size is ~0.2 mag in the vertical direction, which is comparable with the photometric errors. A sinusoid of period ~12 yr is overlaid on the N1 light curve. The same, but with a phase shift of π, is also drawn in the bottom box.

3.4 Proper motion, kinematics and dynamical age of the H₂ jet

Figure 8 (a-d panels) displays the knots identification and their measured proper motions (PMs; panel e) as derived from the analysed H2 continuum-subtracted images at high-spatial resolution (see Sect. 2.6). The shifts of some of the knots are evident in the short movie of Fig. H.2. Thanks to the extremely high-spatial resolution of our images and their 17-yr time baseline, we were able to measure knot PMs down to a few milliarcseconds (mas) per year. These values range from ~2 to ~20 mas yr⁻¹. The results are listed in Table I.1. In particular, the table provides knot IDs, coordinates (as derived from the SOUL map), proper motions (in mas yr⁻¹), tangential velocities (at a distance of 1.64 kpc) and position angles (PAs, from north to east) of the corresponding vectors, and their uncertainties. The blue arrows in Fig. 8 (panel e) show the average displacement of each knot in 100 yr along with its uncertainty (red ellipses). Overall, both blue-shifted (B and C groups) and red-shifted (X and A groups) structures along the jet move away from the source position, roughly following the precession pattern described in Caratti o Garatti et al. (2008, see also their Fig. 11).

On average, velocities increase moving farther away from the source position (groups B1, B2 and B3 on the blue-shifted side with average υ_tg ~ 69, 107, and 115 km s⁻¹ and groups X and A2 on the red-shifted side with average υ_tg ~ 40 and 70 km s⁻¹) until they drop as the flow encounters a slower group moving ahead (i. e. group C1 and A1 with average υ_tg ~ 64 and 55 km s⁻¹, respectively) and collide against it (see Fig. 9). Our analysis of the 3D velocities (see Sect. 4.1) confirms that this is not just a projection effect.

The clearest example is shown in Fig. 10. The fast flow in the B3 region is moving at υ_tg ~ 100–120 km s⁻¹ (see e. g. knots B3c and B3bb) and shocks against the C1hi, C1hj, and C1hk regions (υ_tg ~ 10–70 km s⁻¹, corresponding to polygons HLN13 and HLN14 in Fig. 3), producing a bow shock (knots C1hi, hk, hh, k, and k1), which increases in brightness from the 2003 to the 2020 epoch (middle and bottom panels of Fig. 10, see also the short movie of Fig. H.2). Indeed, one can see from Fig. 7 that these two H₂ emitting regions are among those displaying the largest increase in luminosity (about 1.6 mag).

The average tangential velocity of the H₂ flow is 80±30 km s⁻¹. Here the reported uncertainty is the standard deviation over the whole sample. In fact, by analysing the velocity vectors of all sub-structures in a group, it is clear that tangential velocities and position angles exhibit large spreads (namely tens of km s⁻¹ and degrees, i. e. much larger than the velocity uncertainty of each single knot) within the same group (see bottom and middle panels of Fig. 11). In addition, the scatter in velocity is larger in those regions where knots are more crowded and shocks are therefore expected to interact with each other. This reflects and possibly causes the observed fragmentation along the flow.

There are other interesting features arising from the proper motion analysis of the different regions along the jet. Group X, corresponding to the ring-like feature (PD) and as the H₂ region closest to the source position, looks like a structured expanding ring with the exception of knot X2 (polygon HLPD2 in Fig. 3), which first appeared in the 2012 image and increased in luminosity by 2020. The expansion is clearly visible from the vectors in Fig. 8 and from the large scatter in their position angles (see middle panel of Fig. 11). On the other hand, knot X2, which first appeared in the 2012 image as well, is not co-moving with the expanding ring-like feature but it is rather moving straight away from the source position.

By combining proper motions and projected distance (d) to the source, we also estimated the dynamical age (τ = d/PM; in yr) of each knot. Values are reported in Col. 8 of Table I.1 and plotted against their projected distance to IRAS 20126+4104 in the top panel of Fig. 11. The dynamical ages provide a raw estimate for the ejection time. Actually they provide a lower limit, if matter actually accelerates soon after its ejection, as Fig. 9 would suggest. Figure 12 shows the average dynamical age of each group of knots versus their average projected distance to the protostar. Notably, all the knots have been ejected recently. Their average τ $\left(\overline{\tau }\right)$ ranges from 220±50 yr for group B1 along the jet to 2200± 900 yr for the farthest group A1 towards the south-east. Group X (the ring-shaped region close to the source, labelled as PD in our photometry analysis) has a similar dynamical age as B1 $\left({\overline{\tau }}_X=280\pm 50\text{\hspace{0.17em}yr}\right)$, whereas the $\overline{\tau }$ value of the farthest group towards the north-west (C2) is roughly consistent with that of A1, namely ${\overline{\tau }}_{C2}=1070\pm 300\text{\hspace{0.17em}yr}$. This might indicate that the A1 group decelerated, as is also hinted at by the presence of a large bow shock at the front of the group. The large uncertainty on $\overline{\tau }$ of some groups reflects the large scatter in tangential velocities of the groups as seen in Fig. 11, in particular for groups A1, A2, B3, C1, and C2, namely where the dynamical interaction between knots may be more important.

One puzzling feature displayed by the PM vectors is that the average direction of the farthest knots (A1, C1, and C2) does not intersect the current location of the protostar, but passes northeast of it. If a line coinciding with the proper motion direction of the protostar is drawn, the average directions of the knot proper motions cross it at the earliest locations (C1 and C2 ~ 1″ – 2″ north-east of the protostar and A1 and A2 ~ 3″ north-east of the protostar), in accordance with them being ejected at earlier epochs (see Fig. 13). However, the knot PMs have been derived in a reference frame that appears to be similar to that of the protostar (see Appendix B), so such an effect should not be detected in our analysis. One possible explanation is that the protostar proper motion has decreased (in absolute value) with time.

Fig. 5 Comparison of the continuum emission throughout the various observed epochs. This has been approximated with the Brγ-filter images, except for the runs of July 2000 and August 2006 where broad-band K images corrected for line emission have been used. The flux levels have been adjusted so that the mean of the counts of stars 2 and 4 is the same in all frames. The dimming of the northern lobe (especially region N1) and the simultaneous brightening of the southern lobe in July 2003 are evident.

4 Discussion

4.1 H₂ 3D kinematics

By combining the tangential velocities (υ_tg) derived in this work with the H₂ radial velocities (υ_r) obtained at high-spectral resolution (R = 18 500) with UKIRT/CGS4 by Caratti o Garatti et al. (2008), we are able to infer both the total velocity (υ_tot) and the inclination of its vector with respect to the line of sight (i) for the knots encompassed by the CGS4 slit. It is worth noting, however, that the spectral images of CGS4 have a spatial resolution one order of magnitude worse (slit width 0″.5, seeing ~1″) than the high-resolution images used here. As a result, the spectra extracted by Caratti o Garatti et al. (2008) typically encompass more than one substructure and more than one velocity component is detected as well (see their Fig. 8 and Table 4). To associate radial and tangential velocities of each region, the spatial resolution of the CIAO image (closest in time to the spectral data) was first degraded to 1″ to match that of the spectral images of UKIRT/CGS4. Then, slit-like areas were cut out of the degraded image and the corresponding 1D profiles were extracted and matched to those extracted from the CGS4 spectral images. We then assumed that the peak intensity of the extracted spectra originated from the brightest knot-substructure encompassed by the slit and combined its radial velocity with the tangential velocity of that knot or substructure (given in Table I.1). The total velocity and inclination of each knot are then as follows: $${\upsilon }_{\text{tot}}=\sqrt{{\upsilon }_{\text{tg}}^2+{\upsilon }_r^2}$$(2) $$i=\arctan \left(\frac{{\upsilon }_{\text{tg}}}{{\upsilon }_{\text{r}}}\right).$$(3)

The results are shown in Table I.2, where the identified substructures with total velocities and inclination i with respect to the line of sight are listed. The inclination i of vectors in both lobes ranges from 76° to 98° with respect to the observer, yielding a weighted-mean inclination with respect to the plane of the sky of i_sky = 8°±1°. This result is consistent with the value i ~ 80° derived from the H₂O maser measurements obtained (towards the blue-shifted lobe) close to the source by Moscadelli et al. (2011). Our results therefore confirm that the outflow axis lies close to the plane of the sky at spatial scales ranging from a few hundred mas to ~10″. The ~10° spread in the i values of the knots might arise from the jet precession or it might indicate that our original assumption in associating each peak υ_r to the υ_tg of the corresponding brightest substructure is not completely correct. In any case, such a spread does not really affect the fact that the reported tangential velocities of the knots can be assumed equal to the total velocities (υ_tg/ sin i ≃ υ_tg ≃ υ_tot) within a 3% uncertainty in the worst-case scenario (i. e. i = 76^° towards group C2).

Fig. 6 Photometric variability of jet areas in the H2 2.12 µm line emission. We note that LLN and LLS indicate the northern and southern lobe, respectively, and PD is the ring-like feature (see Fig. 2). In addition, LLS is further subdivided into LLS1 and LLS2, and LLN2 is in turn subdivided into LLN2a and LLN2b. The small open triangles (connected with dotted lines) indicate the fluxes with continuum subtraction from the K′ band when corrected following Appendix G. The photometric errors are smaller than the symbol sizes. The magnitudes have been corrected for bandwidth differences.

4.2 Bow shocks

We tested our data to check that the H₂ emission (morphology, fluxes, and projected velocities) of at least some knots is consistent with being originated in bow shocks. We performed two tests based on a comparison with model predictions. In the first test, the morphology, size, and brightness of the H₂ knots were compared with magneto-hydrodynamic (MHD) models of bow shocks. As an example, we considered knots A1, B3c, and C1h-C1k, as well as the 3D bow-shock models of Gustafsson et al. (2010). We converted the knot photometry into brightness by using the 2MASS zero magnitude flux at K_s, a filter width of 0.023 µm (see Table D.1), and the solid angle of the relevant polygons (correcting for extinction by adopting the A_V values derived by Caratti o Garatti et al. 2008). The morphology of the selected knots is clearly reminiscent of that of bow shocks from a flow parallel to the plane of the sky. As for A1, both its width (~4000 au) and its average brightness (~1–2 × 10⁻⁶ W m⁻² sr⁻¹) roughly agree with a modelled bow shock impacting a medium of density ~10⁵ cm⁻³ with a speed of 50–60 km s⁻¹. Since the actual speed that we have derived is ~45 km s⁻¹, A1 is likely to represent a terminal shock. The width (~ 800 au) and average brightness (~9 × 10⁻⁷ W m⁻² sr⁻¹) of B3c again points to a bow shock propagating in a medium of density ~10⁵ cm⁻³ with a speed of 50–60 km s⁻¹. Finally, C1h-C1k has a width of ~1600 au and an average brightness of ~0.4–1 × 10⁻⁵ Wm⁻² sr⁻¹, which is consistent with a bow shock propagating in a medium of density ~10⁵–10⁶ cm⁻³ with a speed of 50–60 km s⁻¹. We note that the models of Gustafsson et al. (2010) are able to reproduce some of the asymmetry exhibited by the knots on the basis of the magnetic field direction both on the plane of the sky and along the line of sight. The magnetic field strength in the models of Gustafsson et al. (2010) selected here range from between 500 and 5000 µG. It is worth noting that the most remarkable discrepancy between observations and models in the specific cases of B3c and C1h-C1k is the measured knot velocity, which exceeds the hydrogen dissociation limit (~60 km s⁻¹). Indeed, most of the observed knot velocities exceed such a value (see Table I.1). The most likely explanation, as already mentioned, is that these features represent internal shocks, in other words that the flow here is impacting a parcel of gas already in motion, resulting in a lower relative velocity. As a consequence, the measured tangential velocities (and thus total velocities) can be greater than the actual shock velocities.

In the second test, we compared knots Al, B2, and C1 with the ballistic bow-shock model of Ostriker et al. (2001). In order to derive the proper coordinate system used by those authors (i.e., distance z from the head of the shock along the jet axis and distance R from the jet axis), we have developed a software routine that fits their Eq. (22) to the outermost shape of H₂ emission structures resembling bow shocks. For the sake of simplicity, the fit assumes the jet to lie in the plane of the sky (a good assumption, as discussed in Sect. 4.1), a bow-shock speed υ_s = 150 km s⁻¹, a sound speed C_s = 8 km s⁻¹ (which is appropriate for a gas temperature of 10 000 K), and β= 4.1 (Sanna et al. 2012). The fit yields the coordinates of the head of the shock (the centre of the working surface), the projected orientation of the shock axis, and R_j, the radius of the inner driving jet (and of the working surface). Using Table I.1, one can now compute the coordinates (R, z) of the associated knots in the bow-shock reference frame and the components of their projected speed parallel (longitudinal) and perpendicular (transverse) to the jet axis. In turn, these can be compared with the ones predicted by the model (Eqs. (18)–(21) of Ostriker et al. 2001). In particular, given a knot associated with a bow shock, its projected velocity component perpendicular to the jet axis should fall between zero and the predicted radial speed of the outer surface layer at the same distance z from the shock head, whereas the knot-projected velocity component parallel to the jet axis should fall between the predicted mean shell and outer surface layer longitudinal speeds at the same distance z from the shock head. The comparison and the spatial distribution of the shock envelope from the fit is displayed in Fig. 14 for knots Al, B2, and C1. The results of the fit are also listed in Table 2. Figure 14 shows that the morphology of knots Al, B2, and C1 is well fit by the bow-shock outer-shell model. The proper motions of the knots making up group Al appear to be consistent with the model prediction. As for C1, the velocity plots indicate that the location of the working surface at the head of the shock has probably not been estimated correctly and that it should be shifted a little backwards along the jet. After this correction, the proper motions of the knots making up C1 are consistent with the model predictions as well (possibly except for knot C1qd). Some of the proper motions of the knots making up B2 are consistent with the model predictions, while others are not. In this case, one has to take into account that the region around B2 displays a complex pattern of shocks, thus the matter along the jet axis is probably soon entrained in other internal shocks due to subsequently ejected matter. In this respect, knot B2b, the most diverging one, is clearly a bright isolated patch of emission (see Fig. 8) located at the border of the fitted shock outer shell, marking a different shock episode with high probability.

Interestingly, the method also provides an estimate of the aperture angle of the jet from the radius of the inner jet driving the bow shock. The values obtained range from between 2–5.5° (half aperture, see Table 2), which is slightly less than the value obtained from the water masers (~9º; Moscadelli et al. 2011). In addition, the bow-shock axes derived from A1 and C1 intersect the proper motion track of the protostar north-east of it, as found for the average proper motion directions.

Fig. 7 Photometric variability of jet knots in the H2 2.12 µm line emission. We note that LHN and LHS indicate the northern and southern lobe, respectively, and LHPD is the ring-like feature (see Fig. 3). The magnitudes have been corrected to the LUCI1 H2 filter bandwidth.

Fig. 8 Knot identification in the SOUL H2 continuum-subtracted image of the IRAS 20126-4104 flow and knot proper motions. Panel a: zoom-in on the A and X knot regions (red-shifted lobe). Contour levels at 5, 7, 10, 15, 20, and 25σ are overlaid. We note that 1σ corresponds to 8×10−23 W cm−2 arcsec−2. Panel b: zoom-in on the B and C1 knot regions (blue-shifted lobe). Contour levels at 5, 7, 10, 15, 30, and 40σ are overlaid. Panel c: zoom-in on the C1 and C2 knot regions (blue-shifted lobe). Contour levels at 5, 10, 20, 30, 40, 50, 60, 80, and 120σ are overlaid. Knot peaks are indicated by white dots. Panel d: overall view of the IRAS 20126+4104 flow close to the source. The main structures as reported in Caratti o Garatti et al. (2008) are labelled. Panel e: proper motions (PMs) with their uncertainties (blue arrows and red ellipses) in 100 yr of structures and sub-structures along the H2 jet in IRAS20126+4104. The actual observed shifts are approximately one fourth the length of the corresponding arrow. The red circle marks the position (along with its uncertainty) of the protostellar continuum emission at 1.4 mm (Cesaroni et al. 2014). The main structures as labelled in Caratti o Garatti et al. (2008) are also indicated.

Fig. 9 Average tangential velocity vs average projected distance to the source for each group of knots.

Fig. 10 Evolution of the bright bow shock in the C1h and C1k regions. The top, middle and bottom panels show the same area of the sky at the three different epochs (2003, 2012, and 2020). The fast moving flow (υtg ~ 100–120 km s−1) corresponding to B3 shocks against the C1h and C1k regions (υtg ~ 10–70 km s−1) producing a bright bow shock. The PMs of knots are shown as in Fig. 8.

Fig. 11 Tangential velocity (bottom panel), position angle (middle panel), and dynamical age (top panel) vs. distance to the source for each knot. Red labels show knots in the C1 group that might belong to a different flow (see text).

Fig. 12 Average dynamical age vs. average distance to the source for each group of knots.

4.3 Jet precession

A comparison between the continuum emission and the H₂ line emission in Fig. 1 clearly indicates that continuum emission (outlined by the blue-coloured emission) is detected on both the blue-shifted and the red-shifted sides of the outflow, meaning that scattered dust emission from the cavity dug by the protostellar jet is visible on both sides. This is consistent with a jet lying almost on the plane of the sky (as confirmed by our measured 3D velocities, see Table I.2), with a large enough opening angle (half opening ~9º according to the water maser motion modelling by Moscadelli et al. 2011). The relatively bright double-sided cavity, compared to the nearby fainter source S suggests that the protostar may be located near to the front surface of its parental molecular clump, which could explain the low extinction towards both cavities.

The proper motion analysis confirms that all the H₂ features in the north-western lobe are moving away from areas close to the protostar location. In particular, the proper motions of group B knots agree well both in speed and direction with those of the cluster of water masers, which are closer to the protostar location (Moscadelli et al. 2011), confirming that the water masers and group B knots are tracing different episodes of mass ejection from the same source (see Fig. 15).

The knots of group A are roughly symmetrical to those of group C with respect to the protostar location and the velocities are red-shifted as expected. In contrast, their morphology is different from that of knots in group C, possibly due to one or more of the following reasons: (i) the mass ejection is asymmetric; (ii) the south-eastern lobe is impinging on denser ambient gas; and (iii) the south-eastern lobe is more extincted than the north-western one. The lower velocity of knots A (see Table I.2) would suggest deceleration, which is consistent with case ii.

The knots of group B have no symmetrical counterparts in the south-eastern lobe, except for two tiny features south of the ringlike emission X (see Fig. 8). However, a symmetrical counterpart has been detected at millimetre wavelengths. In fact, B1 coincides with a patch of SO emission between −26 and −12 km s⁻¹ with respect to the systemic velocity (Palau et al. 2017), hence it is consistent with the radial velocities derived from the H₂ emission by Caratti o Garatti et al. (2008). A nearly symmetrical patch of SO emission between 11 and 32 km s⁻¹ south-east of the protostar was detected by Palau et al. (2017). This of course discards the possibility of an asymmetric mass ejection for group B.

The HCO⁺(1−0) outflow emission mapped with the IRAM Plateau de Bure (PdB) interferometer by Cesaroni et al. (1997) displays the same south-eastern-north-western morphology as seen in the K band. In particular, the south-eastern lobe, when mapped at red-shifted velocities, exhibits the ‘jaw-like’ morphology found by Smith & Rosen (2005) in their hydrodynamic simulations of slow-precessing protostellar jets. The knots of group A are located at the northern edge of the south-eastern molecular lobe, which corresponds to the most recent ejection episode when compared to the morphology of the hydrodynamical models. Interestingly, the knots of group A coincide with a blue-shifted patch of molecular emission (see Fig. 7 of Cesaroni et al. 1997), which confirms that the outflow is currently almost perpendicular to the line of sight. The reversal of blue-shifted and red-shifted lobes, when the HCO⁺(1−0) radial velocities are closer to the systemic radial velocity, clearly indicates both the jet high inclination and its non-negligible opening angle. This is widely discussed in Cesaroni et al. (2005) and our proper motion analysis definitely confirms this scenario. The north-western outflow lobe does not exhibit a ‘jaw-like’ morphology as the south-eastern one does. This may be due to the fact that the outflow is still eroding the innermost parts of the parental cloud in the south-east, whereas it has partly emerged from the cloud in the north-west. This is confirmed by the 3D velocities of groups A and C, which are 100 km s⁻¹ towards A, but still above ~100 km s⁻¹ towards C (see Table I.2).

The global morphology of the H₂ emission resembles the models by Smith & Rosen (2005) if one assumes that in the case of IRAS 20126+4104, the mass ejection is episodic rather than periodic. The interactions between different bow shocks and the sudden changes of direction predicted by the model are clearly visible. The models of Smith & Rosen (2005) assume a precession angle of 10°, a precession period of 400 yr, and a jet velocity of 100 km s⁻¹, which are close to the values of ~7.6° and 1000 yr derived by Caratti o Garatti et al. (2008), and the velocities obtained in this work. Smith & Rosen (2005) also found that the total H₂ luminosity monotonically increases with time, but apparently only less than ~0.01 mag yr⁻¹, which is consistent with the results of Sect. 3.3. Thus, the similarity between models and observations further supports the idea of a pulsed precessing jet fed by episodic ejection in IRAS 20126+4104.

Fig. 13 Mean directions of the proper motion vectors in the blue-shifted lobe (top panel) and in the red-shifted lobe (bottom panel). The proper motion track of the protostar is marked by the red line and its current location by the red cross. Clearly, the more distant the knot group is from the protostar, the more distant (to the north-east) the intersection of its mean direction is with the track. The protostar is moving towards the south-west.

Fig. 14 Ballistic bow-shock model of Ostriker et al. (2001) compared with knot chains A1 (top row), B2 (middle row), and C1 (bottom row). In the column on the left 6″ × 6″ zoom-ins of the SOUL H2 image are shown with the bow-shock model shape overlaid. In the middle column, the transverse velocities of the involved knots are compared with the mean shell (dotted line) and outer surface layer (solid line) values predicted by the model. The same comparison is made in the column on the right for the longitudinal velocities. The knots whose speed departs from the model predictions the most are marked in red and labelled. For the sake of simplicity, the outflow is assumed to lie on the plane of the sky.

Fig. 15 Zoom-in on Fig. 8 around the source showing H2O maser proper motions (magenta arrows) as measured by Moscadelli et al. (2011) and H2 proper motions (blue arrows) from this paper. Red ellipses mark the uncertainties on PMs. The black circle shows the position of the protostellar continuum at 1.4 mm from Cesaroni et al. (2014).

4.4 A possible wide-angle wind shocking the outflow cavity

A clue as to the nature of the ring-like feature X can be found in the models of Smith & Rosen (2005). These authors find a ‘stagnation point’ at the apex of the cone of matter which is not affected by the jet (inside the precession cone). Gas accumulates at the stagnation point originating a disc-like feature which then releases matter with large lateral speed. This might resemble the expansion of the ring-like feature shown in Fig. 8. Alternatively, an explanation for this expanding ring feature, here traced by H₂ shocked emission, is that we are actually seeing a wide-angle wind, which shocks the outflow-cavity walls and blows away the outer envelope, further opening the cavity. These wide-angle winds are predicted by MHD disc-wind models (see e. g. Kwan & Tademaru 1995; Shang et al. 2006) and have recently been detected at millimetre wavelengths with the ALMA interferometer in a few low-mass protostars (e. g. HH47 and DG Tau B; see Zhang et al. 2019; de Valon et al. 2020). Indeed, such models predict that MHD disc winds have both a fast-moving collimated component (the jet itself) and a slow-moving less dense wide-angle component that blows into the ambient matter, forming an expanding outflow shell. In our case, both the ring-shaped geometry and the lower velocities observed in the ring-shaped structure of group X are consistent with this picture. Additionally, the fact that group X has lower H₂ column density and mass flux (up to an order of magnitude) with respect to the main jet (see Table 5 of Caratti o Garatti et al. 2008) further supports this idea. Palau et al. (2017) detected complex molecules at sub-arcsecond resolution using the PdB interferometer. Notably, one of the detections spatially matches the H2 emitting region of group X. By modelling the gas emission and from the X-shaped morphology of the SO low-velocity component and HNCO, Palau et al. (2017) conclude that such molecules are tracing shocks along the outflow cavity walls of IRAS 20126+4104, which is consistent with our analysis. Furthermore, Zhang et al. (2019) suggest that these shells are the results of an intermittent wide-angle wind that is produced by episodic accretion and ejection. This would also be consistent with knot fragmentation seen in the IRAS 20126+4104 jet. Finally, by fitting an ellipse to the ring-shaped feature, we find that its major axis has a position angle of 29°± 5°, namely orthogonal to the jet position angle (−60.9°) derived by Caratti o Garatti et al. (2008). If one assumes circular symmetry for the ring-shaped feature, it is also reasonable to expect its symmetry axis to be roughly parallel to the jet axis. In this case, the inclination of the ring-shaped feature to the line of sight can be derived from the ratio of the ellipse axes, which is also equal to the symmetry axis inclination to the plane of the sky. The feature inclination with respect to the observer ranges from between 68° and 54°, meaning that under the hypothesis of circular symmetry the inclination of the feature symmetry axis with respect to the plane of the sky would be quite different from ~8°, as estimated for the jet in Sect. 4.1. A possible explanation to this conundrum might be that we are not seeing the whole shocked region along the red-shifted outflow cavity, namely that part or most of the lobe surface towards the observer is obscured by high visual extinction.

4.5 Disc oscillations and the anti-correlated continuum variability

The continuum variability observed in the outflow cavity could in principle have been produced by a short accretion increase that occurred between 2000–2003. However, such an event would have either illuminated both outflow cavities or left one of them unaffected in the unlikely event that the burst had just occurred on one side of the disc. This scenario is therefore difficult to reconcile with the clear anti-correlation of the K-band continuum light curves between the north-western (dimming lobe, see Fig. 3) and the south-eastern lobe (brightening lobe, see Fig. 3). Thus we have explored the alternative view of periodic flux oscillations related to the disc and jet precession. As the jet precession originates from a simultaneous disc precession, the inner circumstellar disc might periodically eclipse one hemisphere of the central protostar at different times, as seen from the north-west and south-east, producing such an out-of-phase variability. As shown in Fig. 4, an anti-correlated sinusoidal variation with a period of ~12–18 yr roughly overlays the K-band continuum light curves. Unfortunately, no clear-cut conclusions can be drawn concerning periodicity due to the sparse photometric data coverage after 2003. Nevertheless, in the following we assume that a periodic scenario holds to develop a picture that can be easily checked by future observations.

By modelling the dynamics of a protostellar disc in a binary system, Bate et al. (2000) found that a circumstellar disc misaligned with the orbital plane of the binary system exhibits a wobbling with a period about half the orbital period and precesses with a period nearly equal to 20 orbital periods. If the continuum variability is caused by the disc wobbling, this would indicate the presence of a companion with an orbital period of ~24–36 yr, which would also produce a precession with a period of ~480–720 yr, not too far from the inner-jet precession period estimated by Caratti o Garatti et al. (2008), that is ~1000 yr. It is worth noting, however, that these authors assumed a jet velocity of 80 km s⁻¹ and a distance to the source of 1.7 kpc. By using the revised, more accurate distance of 1.64 kpc and a jet velocity of 110 km s⁻¹, their precession period would be about 700 yr, which nicely matches the period estimated by us using the recipe of Bate et al. (2000).

As the central protostellar mass derived by Chen et al. (2016) from the Keplerian velocity pattern of the disc roughly agrees with that derived from the system bolometric luminosity, the companion, if any, must be a low-mass star. Assuming a central protostellar mass of ~12 M_⊙, it should be located at a distance of ~ 19–25 au (i.e. subtending an angle <0″.02). It is also worth noting that, according to Bate et al. (2000), the plane of the orbit of the supposed companion must be misaligned with respect to the disc plane. We can use Eq. (1) of Terquem et al. (1999) to obtain some more constraints. Of course, this equation cannot be applied to the circumstellar disc as a whole in our case, with a precession period ~1000 yr. In fact, it assumes that the disc precesses as a rigid body and this requires that the time needed to cross the disc at the speed of sound is much less than the precession period. However, a disc of ~1000 au would be crossed in ~5000 yr at a sound speed of 1 km s⁻¹.

By assuming a primary star mass of 12 M_⊙, an orbit inclination of 30° with respect to the disc, an orbital radius of 19 au, and Eq. (1) of Terquem et al. (1999) yields a disc radius of ~15 au for a secondary star of 1 M_⊙ and a precession period of 1000 yr or a secondary star of 2 M_⊙ and a precession period of 500 yr. Vaidya et al. (2009) estimated that the dust sublimation radius is between 3–17 au for a 10 M_⊙ protostar. However, Gravity Collaboration (2020) and Koumpia et al. (2021) measured a smaller range of values (from 0.6 to 10 au) for a sample of massive protostars (~ 10–20 M_⊙) using NIR interferometry. In the end, a jet precession with a period of ~700 yr is consistent with the scenario of a nearby low-mass companion inducing oscillations in the innermost region of the circumstellar disc, which might also explain the anti-correlated K-band continuum brightness variations. An analysis of the NIR fluctuations could then become a powerful tool to understand the circumstellar structure of the protostar and would call for a more frequent monitoring.

4.6 Origin of the parsec-scale outflow morphology

The scenario pictured above does not entirely match the data on larger scales. On a scale of ~1 pc, both H2 and CO emission display a more extended north-south outflow centred on the protostar position (see Shepherd et al. 2000). The southern lobe is red-shifted, has a jaw-like morphology, and the south-eastern lobe of the HCO⁺ and H₂ outflow lies at its south-eastern border, indicating, on a larger scale, a different precessing geometry with respect to the closest H₂ flow, which represents the most recent mass ejection episodes. The northern flow (blue-shifted) is smaller than the southern one, confirming that the outflow is still moving inside the cloud in the south, but has probably emerged outside in the north. The jet axis position angles of the closest (−60.9^°) and farthest (−37°) H₂ features differ by 22.1°. This is difficult to reconcile with the precession period (~ 1000 yr) and opening angle (~14°) of the inner jet derived by Caratti o Garatti et al. (2008), in that the CO morphology, along with the presence of optical and NIR line emitting knots further from the protostar, would imply a much larger precession period (≥64 000 yr) and opening angle (37^°; Cesaroni et al. 2005). Indeed, Caratti o Garatti et al. (2008) were unable to find a superposition of two different precessions by simultaneously fitting the positions of the H₂ knots along the whole flow. Nevertheless, the flow agrees with a jet whose projected axis was initially in a roughly north-south direction with an inclination angle of ~45^° with respect to the observer, which has then moved to the current north-western-south-eastern direction (roughly on the plane of the sky).

A possible explanation is that there are two (or more) different outflows, a younger one and an older one. In fact, two other signatures of outflowing matter have been evidenced near IRAS20126+4104. However, one of them, source S, is unlikely to be the driving source of an older outflow, as jet signatures could not have remained obscured in the NIR after some 10 000 yr of mass ejection activity out of the plane of the sky.

An alternative explanation is that a close encounter with a passing-by massive body has caused the disc axis to twist by ~22^° on the plane of the sky and by ~35^° away from the observer line of sight. According to Bate et al. (2000), a close encounter, particularly if occurring during a phase of intense accretion, can cause the outflow axis to tilt, mimicking a precession. The likelihood of such an encounter is of course higher if the protostar is embedded in a young star cluster. In this case, Proszkow & Adams (2009) found that the number of close encounters per cluster member is given by $$\Gamma ={\Gamma }_0{\left[b/1000\text{au}\right]}^{\gamma }$$(4)

where b is the smallest encounter distance. These authors computed Γ₀ and γ for a wide parameter space populated by young stellar clusters, including the number of members. As for IRAS20126+4104, Montes et al. (2015) estimate that the cluster it is embedded in is made up of at least 43 young stars within a radius of 1.2 kpc. Given the low X-ray sensitivity this is bound to be a lower limit, so that we can assume that the cluster has at least 100 members. As they provide Γ per million years, this needs to be multiplied by 0.1 since a time span of 100 000 yr seems more appropriate for this case. Under these assumptions, Γ ranges from between ~1–2 % for b = 1000 au, which, although small, indicates that this scenario is not implausible.

It is unlikely that this body is the hypothesised nearby low-mass companion (see previous section), because a similar precession pattern is also detected in the southern and farthest part of jet (see left panel of Fig. 7 in Caratti o Garatti et al. 2008). This suggests that if the jet precession is caused by a companion, the latter has to be orbiting around the central source ever since. A likely candidate for such a close encounter could then be the intermediate-mass medium infrared- and radio-source S, located 1″ away. The observed tilt of the disc axis is compatible with the system geometry and with source S moving from the northeast to south-west and passing by IRAS 20126+4104. The close encounter might also have caused the similar orientation of the jet-like radio structure associated with S. Indeed the magnitude of the disc axis tilt depends on the angle between the disc plane and the orbital plane, the mass of the star and the distance between the disc and the passing-by object (see e. g. Papaloizou & Terquem 1995).

If a close encounter with source S has caused the jet axis to tilt, then source S must have passed by the protostar near enough to affect the disc plane but at least at the escape velocity. Another obvious constraint is that the close encounter must have occurred earlier than the ejection of the outer knots of the H2 inner jet (the average dynamical ages of groups A1 and C2 are 2200±900yr and 1070±300, respectively, see Sect. 3.4 and Fig. 12), but later than a CO outflow dynamical time ago (~64 000 yr, Shepherd et al. 2000). Assuming a closest distance of 1000 au, the escape velocity would be ~6 km s⁻¹ even for a total central mass of 20 M_⊙ (star plus disc). Typical stellar 1D velocity dispersions in young clusters are ~1–3 km s⁻¹ (e. g. Kuhn et al. 2019), thus such a value would not be implausible for two young stellar objects that originated in the same region, whereas much higher values are unlikely. The protostar and S are now at a projected distance of ~1640 au, which is a lower limit for the actual distance, and would then be crossed at the escape speed in >1300 yr. This timing is consistent with the two age constraints.

A less likely alternative is that S has been trapped in an elliptic orbit, namely IRAS 20126+4104 would then be a triple system, which, typically, is not stable. In this respect, we note that Shepherd et al. (2000) found that a star roughly the same mass as the protostar in a circular orbit with R ~ 1400 au (i.e. similar to the projected distance between the protostar and S) would be able to cause the slower precession. In fact, Sridharan et al. (2005) proposed that source S is the stellar companion responsible for the slower precession. In this case, the orbital motion would imprint the outflow as well. The effects of orbital motions on outflows were discussed by Masciadri & Raga (2002). By assuming that source S has a mass equal or less than the protostar (12 M_⊙) and a circular orbit with a radius greater than the projected distance of source S to the protostar, we can use Eq. (1) of Masciadri & Raga (2002) to estimate a period >38 000 yr, which is of course roughly the same as that of the induced precession for any outflow oscillations due to the source orbital motion. However, from Eq. (13) of Masciadri & Raga (2002), one can see that the oscillations would be confined to a cone of aperture ~1^°, which is negligible compared to the observed wiggling cone. Thus, the existence of an orbital motion cannot be inferred from the outflow morphology. Another possibility is that accretion of matter with different angular momenta has caused the disc inclination to change in time (Matsumoto et al. 2017).

We can conclude that the jet from IRAS 20126+4104 displays one faster precession motion and exhibits the signature of either a close stellar encounter that occurred less than ~64 000 yr ago, or of a slower precession motion. These scenarios need to be investigated through more detailed theoretical models of the circumstellar and stellar environments, in particular, by simulating the effects of a multiple system composed of a nearby low-mass star affecting the innermost region of the disc (this work) and a more massive, further companion or cluster member affecting the whole disc (Shepherd et al. 2000), and of asymmetric accretion (Shepherd et al. 2000) as well.

5 Conclusions

By means of new high-spatial-resolution NIR images (K_s, H2, and Brγ filters) obtained at the LBT with the LUCI1 instrument and the new AO system SOUL, combined with published imaging and spectroscopic data, we have presented a multi-epoch kinematic and photometric study of the IRAS 20126+4104 H₂ jet and outflow. Two other observing runs with AO (2003 and 2012) providing high-spatial-resolution data have been used for our analysis.

The main results of this paper are the following:

• Thanks to the extremely high-spatial resolution of our images and their 17-yr time baseline, we have been able to measure knot proper motions down to a few mas per year, with values ranging from 1.7 to 20.3 mas yr⁻¹, which translate into υ_tg = 13–158 km s⁻¹ at a distance of 1.64 kpc. The average υ_tg of the flow is 80±30 km s⁻¹. The large scatter in velocity is likely due to the strong interaction between knots along the flow.

• The dynamical age of the knots increases moving away from the source, ranging from ~200 to ~4000 yr, indicating that the inner jet is relatively young.

• By combining tangential and radial velocities (the latter inferred from Caratti o Garatti et al. 2008) of the H₂ knots, we have reconstructed the 3D kinematics of the jet, which roughly expands along an axis with an average inclination angle of 8°±1° with respect to the plane of the sky. Overall the inferred 3D kinematics is consistent with the jet precession model presented by Caratti o Garatti et al. (2008) with a ~700 yr period and agrees with the motion of the water masers studied by Moscadelli et al. (2011).

• Some of the larger knot groups are consistent with the prediction of bow-shock models. A fit to a ballistic bow-shock model provides estimates for the half-aperture of the inner driving jet in the 2–5.5° range.

• Both the average proper motion directions and the bow-shock axes of the outer knot groups do not intersect the protostar position, but cross its proper motion projected track ~1″–3″ to the north-east. As the knot proper motions are roughly computed in the protostar framework and the protostar moves towards the south-west, this may suggest that the protostar proper motion has decreased in time.

• Our multi-epoch photometry of the north-eastern and southwestern outflow cavities shows an anti-correlated variability of the NIR continuum with a possible periodicity of 12–18 yr. We argue that such a variability could be caused by the wobbling of the inner disc, which would alternately eclipse the two hemispheres of the central protostar. Such wobbling could be caused by a nearby low-mass star companion with an orbital period in the range of 24–36 yr, which would also produce the precession of the jet.

• The total H₂ 2.12 µm luminosity does not exhibit time variations at a ~0.4 mag level, which is consistent with hydro-dynamical models of outflows. However, some knots exhibit clear line flux increases (of ~1 mag), notably one associated with the bow-shock C1, and two knots exhibit a flux decrease >0.5 mag.

• The longer precession time of 64 000 yr found on the basis of CO emission imaged on the parsec scale by other authors can be reconciled with our estimate by assuming either a past close encounter with a relatively massive object or the presence of a second precession caused by a companion orbiting outside the disc. We have shown that in both cases the nearby intermediate-mass source S is a suitable candidate.

Our analysis demonstrates that multi-epoch high-spatial-resolution observations are a powerful tool to unveil the physical properties of highly embedded massive protostars. In this respect, a photometric monitoring of IRAS 20126+4104 could provide valuable insights into its circumstellar environment.

Movies

Movie 1 associated with Fig. H.1 (aa45235_fig21_movie) Access here

Movie 2 associated with Fig. H.2 (aa45235_fig22_movie) Access here

Acknowledgements

A.C.G. has been supported by PRIN-INAF MAIN-STREAM 2017 “Protoplanetary disks seen through the eyes of new- generation instruments”. F.M. and A.C.G. also acknowledge support from PRIN-INAF 2019 ”Spectroscopically tracing the disk dispersal evolution (STRADE)”. The LBTO AO group would like to acknowledge the assistance of R.T. Gatto with night observations. This work is based on observations made with the Large Binocular Telescope. The LBT is an international collaboration among institutions in the United States, Italy, and Germany. LBT Corporation partners are: The University of Arizona on behalf of the Arizona Board of Regents; Istituto Nazionale di Astrofisica, Italy; LBT Beteiligungsgesellschaft, Germany, representing the Max-Planck Society, The Leibniz Institute for Astrophysics Potsdam, and Heidelberg University; The Ohio State University, representing OSU, University of Notre Dame, University of Minnesota, and University of Virginia. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work has also made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Appendix A SOUL performances

The seeing correction achieved with SOUL was assessed by fitting two Moffat distributions (along two perpendicular axes) to stars 2, 4, 5, 6, and 8 (see Fig. C1) in the final images. This allowed us to determine the variation of the PSF width as a function of distance from the AO guide star (star 8). The same was done on the CIAO and FLAO(PISCES) images to evaluate the performance improvement. An example of Moffat fit is shown in Fig. A.1.

Fig. A.1 Profile of the PSF of the AO guide star along the right ascension axis (black line) in the SOUL H2 image overlaid with the best-fitting Moffat distribution (blue line). The plate scale is 0″.015 pix−1.

The maximum and minimum FWHMs obtained from the five stars in the H2 images from SOUL and FLAO are shown in Fig. A.2 as a function of distance from the AO guide star. The improvement in the correction degree achieved with SOUL compared to FLAO(PISCES) is evident. As expected, the FWHM increases with distance.

The overall results for the three sets of images (CIAO, FLAO, and SOUL) are listed in Table A.1. As for the H2 and Brγ filters, the correction has clearly improved from 2003 to 2013, with SOUL achieving the best results (also considering that FLAO and SOUL operated roughly in the same seeing conditions). As for the K_s filter, our SOUL images were obtained in nights with seeing conditions not as good as those in the other available nights. The values listed were obtained on September 4, with the best seeing of the three nights occurring when K_s was acquired.

Appendix B Astrometry and image mapping

The astrometric calibration of the SOUL H2 image was achieved by a two-step process, using the iraf tasks ccxymatch and ccmap. First, we calibrated the large field H2 image obtained with the TNG in 2001 by matching the stars in the frame and the Gaia EDR3 catalogue (Gaia Collaboration 2016, 2021). From this calibrated image we derived the astrometric information for all the stars in the small FoV of the SOUL images. To minimise any degradation of the astrometric accuracy due to the proper motion of the stars in the two fields, we first noted that Gaia EDR3 provides the motion with good accuracy for two stars in the SOUL field, namely number 4 and 8 of Fig. C1. These are −4.307 ± 0.048 mas/yr (RA) and −4.503 ± 0.051 mas/yr (DEC) for star 4, and −4.150 ± 0.015 mas/yr (RA) and −4.544 ± 0.017 mas/yr (DEC) for star 8. The proper motion of IRAS 20126+4104 can be obtained from water masers (L. Moscadelli, private communication) and is ~ − 4.1 mas/yr (RA) and ~ −4.1 mas/yr (DEC), which is similar to that of the two stars. Therefore, we derived the astrometric solution for the TNG image using only stars with proper motions between 0 and −8 mas/yr both in RA and DEC from the Gaia catalogue, so that the reciprocal positions of the star selected should not have changed by more than 0.1″ from 2001 to 2020.

Fig. A.2 Maximum and minimum FWHMs obtained from the Moffat fits to the PSFs as a function of distance from the AO guide star in the SOUL (red triangles) and FLAO (blue triangles) H2 images.

We checked that none of the eight stars in the SOUL image (Fig. C1) had shifted too much from each other by constructing triangles from each triplet of stars and measuring the relative changes of their sides between the CIAO and the SOUL image. We found that the relative changes imply side length differences ⪝0.05″, confirming that the stellar reciprocal positions have changed by 0.05″ at most, making them an ideal reference grid for the astrometric calibration of all the small field images. The absolute positions in the calibrated H2 SOUL image are then in the ICRS standard (epoch 2015.5) and should be more accurate than ~ 0.1″.

Finally, we registered and re-gridded all the other images to the H2 SOUL one by computing the transformations with the iraf tasks geomap (using the eight stars as the reference grid), and applying them with geotran so that all registered images matched each other. In the transformations we only considered a translation, a scale factor, and a rotation to avoid any possible deformation due to the small differences in the proper motions of the reference stars (although, as discussed above, they are probably not significant).

Appendix C Photometric stability

As not all the images collected were obtained along with nearby standard stars’ observations, a photometric monitoring of both line and continuum jet features must rely on stars common to all the frames, used as local standards. In turn, one has to check the photometric stability of these local stars to select the most stable ones. The candidates suitable to be used as secondary standards are labelled in Fig. C.1. It is important to note that star 8 is the AO guide star in all AO-assisted observations examined in this work and is out of the PISCES frames. Our analysis, summed up in Appendix E, shows that stars 2 and 4 are probably stable within ~ 0.1 mag, although they may have entered the non-linear regime in the CIAO frames. The jet photometry is therefore referenced to stars 2 and 4.

An effect to take into account in using only two local standards is the error in the obtained zero points caused by small differences in the PSF, the aperture, and the sky estimate from image to image. We have adopted apertures whose radii encompass the position where the signal is ⪝5 – 10 % of the peak, which means ⪞1 FWHM of the PSF for the lower-resolution images (UKIRT, TNG) and a few FWHMs of the PSF for the AO assisted images (to compensate for the varying Strehl ratios typical of these). To estimate the error due to any difference in the fractions of flux recovered in different images, we repeated the photometry doubling the aperture size, which should largely overestimate the effect we are concerned with. In fact this causes differences ⪝0.2 mag in all cases and bands, except for star 2 in the low-resolution images ( ⪝0.5 mag), as in this case a larger aperture encompasses more diffuse emission from the south-eastern lobe of the jet resulting in a flux increase. We can therefore expect in this case as well that the fraction of stellar flux recovered by doubling the aperture would be ⪝0.2 mag if there was no diffuse emission.

Appendix D Band effects on the photometry

The photometric variability of the stars in the small FoV of the AO-assisted images (see Fig. C.1) was measured in all available bands (K_s, K′, H2, Brγ, and K_cont; see Table D.1) by simply adopting the K_s values provided by (or computed from) the 2Mass PSC, irrespectively of the band, to derive the zero point. A simple procedure can be used to estimate the effects of this simplified zero point estimate on a source magnitude in each band, at least at the first order.

Given a band centred at λ_c with an FWHM of ∆λ, it can be shown that a second-order approximation (in ∆λ/λ_c) is given by $${\text{mag}}_{\text{i}}=C_i-2.5\times \log \left[F_{\lambda }\left({\lambda }_c\right)\text{Δ}\lambda \right],$$(D.1)

Fig. C.1 Brγ-filter image of IRAS 20126+4104 obtained with LUCI1 and the AO system SOUL, with the identification number of the stars available as photometric references.

where C_i depends on the band i in the adopted standard system. This requires that the source has a continuum spectrum with no lines that is well approximated by F_λ (λ) = A × λ^γ inside the band being considered.

In the NIR, ∆λ/λ_c < 1 so the second-order approximation is valid. The terms C_i can be obtained if the zero magnitude flux F_λ0(λ_c) is known, as in this case $$C_i=2.5\times \log \left[F_{\lambda ,0}\left({\lambda }_c\right)\text{Δ}\lambda \right].$$(D.2)

The zero magnitude fluxes from 2Mass can be used to define a source which has magi = 0 in every NIR band linked to this standard system. By simply fitting $$F_{\lambda ,0}\left({\lambda }_c\right)=A_0\times {\lambda }_c^{{\gamma }_0}$$(D.3)

to the 2Mass J, H, K_s zero magnitude fluxes, one obtains log(A₀) = − 12. 16 and γ₀ = −3. 59, with F_λ in W cm⁻² µm⁻¹, which approximates the zero magnitude fluxes by better than 5 %.

This allows one to address the first issue, that is to say how the magnitude of a standard star should be changed to account for a different passband overlapping the original one. Given a filter j with central wavelength λ_j inside the K band and FWHM = ∆λ_j, its zero point in this generalised 2Mass standard system is $$C_j=2.5\times \log \left[F_{\lambda .0}\left({\lambda }_j\right)\text{Δ}{\lambda }_j\right]$$(D.4)

and the magnitude of a continuum source with flux density F_λ(λ_j) measured through a slightly different bandwidth ∆λ (standard star and target observed through slightly different passbands with the same central wavelength but a different width) is given by $${\text{mag}}_{\text{j}}^{\prime }=2.5\times \log \left[F_{\lambda .,0}\left({\lambda }_j\right)/F_{\lambda }\left({\lambda }_j\right)\right]+2.5\log \left[\text{Δ}{\lambda }_j/\text{Δ}\lambda \right]$$(D.5) $${\text{mag}}_{\text{j}}^{\prime }-{\text{mag}}_{\text{j}}=2.5\log \left[\text{Δ}{\lambda }_j/\text{Δ}\lambda \right],$$(D.6)

where mag_j is the magnitude that would be measured in the exact photometric system of the standard star: $${\text{mag}}_{\text{j}}=2.5\times \log \left[F_{\lambda ,0}\left({\lambda }_j\right)/F_{\lambda }\left({\lambda }_j\right)\right].$$(D.7)

If the standard star and target are measured through the same filter, then Eq (D.7) can be used to recompute the standard star magnitude once a new F_λ,0 is determined at the modified central wavelength from Eq. D.3. The effects of absorption or emission lines in the spectra can be neglected as long as the flux removed or added by the lines is negligible compared to the total flux in the passband.

The standard star imaged in the SOUL run, FS149, is an A2 star and it can be seen that the slightly different K_s band of the LUCI1 filter and even the H2 LUCI1 band do not introduce a significant difference compared to its 2Mass K_s value (less than 0.01 mag). However, since the jet photometry has been done using stars 2 and 4 as local standards by assuming their K_s values in all bands, the band effects must also be estimated for these stars. By fitting Eq (D.3) to the 2Mass J, H, K_s fluxes of star 4 we obtained log(A₀) = −17.62 and γ = −1.47 (with F_λ in W cm⁻² µm⁻¹). Unfortunately, this procedure does not yield a good approximation for star 2 (meaning that its spectral energy distribution is more complex), so, from H and K_s only, we obtained log(A₀) = −20.19 and γ₀ = 5.31.

By using Eq. (D.7), it can easily be shown that the corrected magnitudes are within ~ 0.05 mag of the K_s value for star 4 (except for K_cont which is ~ 0.1 mag brighter) and within ~ 0.2 mag for star 2 (except for K_cont which is ~ 0.5 mag brighter). The corrected magnitudes and the measured instrumental magnitudes of the two stars were then used to compute the zero points for each image. We note that the zero points obtained from each of the two stars in the TNG and SOUL images, assuming the 2Mass K_s values in all bands, always differ by ~0.05 mag (except for K_cont for which the difference is ~ 0.1 mag), indicating that the spectral energy distribution (SED) of star 2 is not as steep as in the above approximation when limited inside the K band. Conversely, the same zero points estimated from the CIAO and PISCES images (i.e. derived separately from star 2 and 4) differ by ~ 0.1 – 0.2 mag in the H2, Brγ, and K_s bands. As for CIAO, the PSF of the two stars suggests they are in a slightly non-linear regime. In addition, star 4 lies very close to the frame edge in the PISCES images.

Line emission is measured with photometry on narrow-band images encompassing the line after subtraction of the continuum contribution estimated from a nearby, ideally line-free, band. As continuum emission (mostly due to dust scattered emission) changes with wavelength, one expects a systematic error in the continuum estimate depending on the central wavelength of the band used. From Eq. (D.5), one can derive a correction to apply to the estimated continuum magnitudes to obtain the right values in the 2Mass K_s band. By assuming a continuum flux increasing as λ⁴, this correction would be 0.04 mag for the Brγ filters, ~ −0.14 mag for the H2 filters, ~ −0.15 mag for the K′ and K98 filters, and ~ 0.43 mag for the K_cont filter.

In conclusion, at least for filters K_s, K98, K′, H2, and Brγ the continuum photometry can be carried out simply by assuming the 2Mass K_s values for the two local standard stars with errors < 0.2 mag. We also note that, provided the two stars do not vary significantly, the choice of the standard magnitudes assigned to the reference star for the zero point computation does not affect the relative photometry in the field for a given band.

Nevertheless, it can be shown that to derive the line emission inside a passband one needs to account for the band width. If L = ∫ F_λdλ is the line flux, Eq. D.7 can be generalised as $${\text{mag}}_{\text{p}}=-2.5\times \log \left[L/\left(F_{\lambda ,0}\text{Δ}\lambda \right)\right].$$(D.8)

Since the width of the H2 passband of the LUCI filter is significantly different from that of the TNG and CIAO filters (see Table D.1), this means that $${\text{mag}}_{\text{LBT}}={\text{mag}}_{\text{TNG,CIAO}}+2.5\log \left(\text{Δ}{\lambda }_{\text{LBT}}/\text{Δ}{\lambda }_{\text{TNG,CIAO}}\right)$$(D.9) $${\text{mag}}_{\text{LBT}}={\text{mag}}_{\text{TNG,CIAO}}-0.36$$(D.10) $${\text{mag}}_{\text{LBT}}={\text{mag}}_{\text{PISCES}}+\mathrm{0.18.}$$(D.11)

It is therefore important that this correction is applied to the H2 line magnitudes before they can be compared.

Appendix E Photometric variability of local standard stars

As not all the observations were performed providing a suitable nearby standard star, one needs to resort to local stars that are common to all frames to do photometry. Thus, only seven stars are available (see Fig. C.1, as star 8 was outside the PISCES frame). To be used as secondary standards, their photometric variability must be checked. In fact, stars 2, 4, 6, and 8 were detected in X-rays with Chandra (sources 50, 49, 41, and 37, respectively, in the catalogue of Montes et al. 2015; and sources 192, 187, 175, and 166, respectively, in the catalogue of Townsley et al. 2018). This is consistent with these stars being Weak-line T-Tauri Stars (WTTSs) associated with the protostar cluster. Typically, WTTSs exhibit very stable long-time variability at optical wavelengths rarely exceeding a few tenths of magnitude (Grankin et al. 2008). On the other hand, a small fraction of classical T-Tauri stars (CTTSs) can exhibit larger long-time photometric variations (Grankin et al. 2007). All four stars did not display X-ray variations, but star 8 is listed as ’possibly variable' (Townsley et al. 2018). The lack of X-ray flares during the Chandra exposure may lead one to discard the identification as T-Tauri stars or magnetically active dwarfs; nevertheless, one has to keep in mind that this holds only for the ~ 39 Ks of integration. In any case, as photometric variations of a few tenths of magnitude in the secondary standard stars can significantly bias the light curves of the jet features, the brightness stability of the seven available stars must be checked.

We show our photometry of the eight stars in Fig. E.1, as summarised in Table E.1. The zero points to calibrate them were derived in different ways. The SOUL Ks image was calibrated using standard star observations. The TNG and UKIRT images were calibrated by cross-correlating the relatively large field photometry and the 2Mass PSC. We note that for all bands we adopted the 2Mass K_s magnitudes (for an estimate of the error implied, see Appendix D). Finally, as the small field PISCES and CIAO images do not have standard star observations associated with them, we used star 4 as a reference assuming its brightness is K_s = 11.977, the 2Mass value.

In the last run (SOUL), the target and the standard star FS 149 were imaged in K_s during three different nights (see Tables 1 and E.1), so the different values shown in Fig. E.1 are essentially due to errors on the zero point caused by non-perfectly photometric nights. On JD 2459155.6 and 2459157.6 the standard star was imaged in five different areas of the detector and its instrumental magnitudes remained within ~ 0.03 and ~ 0.2 mag, respectively. On JD 2459159.6, the standard star was imaged in only two positions on the detector and its photometry differs by ~ 0.02 mag.

Figure E.1 indicates that stars 2 and 4 should be stable within ~ 0. 1 mag and these were therefore used as the local standards. Nevertheless, we caution that stars 2, 4, and 8 may be in a nonlinear regime in the CIAO images, leading to an overestimate of the fluxes measured for the other stars on that date.

Appendix F Photometry of the continuum-emitting jet regions with polygons

The photometric variability of the jet in the NIR continuum was analysed by performing aperture photometry. The used images were first registered to and projected onto the SOUL H2 frame grid. The jet was decomposed into different parts and, due to the image transformation, this only required the definition of a polygon for each part in the reference grid of the SOUL H2 image. The adopted polygons and their names are shown in Fig. 2.

The selected polygons are large enough to be insensitive to changes in the morphology of the emitting regions due to their proper motions. The continuum flux inside each polygon was measured using the task polyphot in IRAF. This was done on the Brγ images; when unavailable, we used the K images (UKIRT 2000 and TNG 2006) after subtracting the derived H2 line emission (we used TNG 2001 to estimate the correction for UKIRT 2000). For the epochs when only lower spatial resolution images are available, we also subtracted the stars from the frames by using daophot. We note that using the transformed images also allows one to select the same sky regions to estimate the background to subtract from the raw flux inside the polygons.

We carried out a few tests to assess the errors affecting our photometry. First, we checked on the UKIRT K98 image that performing the photometry on the transformed image would be equivalent to transforming the polygons to fit the frame grid of the original image and perform the photometry on this image (i.e. the task geomap conserves the flux), with only the photometric errors slightly underestimated as a result. Then, we measured the differences between fluxes in the low spatial resolution images with and without star subtraction, which are ⪝0.05 mag, except for N2 in the TNG 2001 image where the difference is ~ 0.1 mag. We also estimated the effects of the different spatial resolutions by convolving the SOUL Brγ image with a Gaussian distribution to obtain a spatial resolution similar to that of the TNG images and repeating the photometry. The flux differences between the two images are ⪝0.05 mag, except for the ring-like feature whose flux is overestimated by 0.2 mag in the resolution-degraded image. Finally, we measured the flux differences for TNG 2003 between the Brγ and the K_cont images, which are 0.3 ± 0.1 mag, hence they are consistent with the estimate given in Sect. D (continuum flux increasing with wavelength).

Fig. E.1 Photometric variability of the stars in the SOUL FoV. The star ID (with reference to Fig. C.1) is indicated on the right side of the panel. The filters are indicated with symbols as explained at the top of the panel. The data taken on JD 2452830.5 (CIAO) and JD 2456099.8 (PISCES) have not been calibrated and the magnitudes were computed by setting star 4 to 11.977 mag in all bands.

As for the K98 (UKIRT 2000) and K′ (TNG 2006) images used to measure the continuum flux due to the unavailability of Brγ images, it is important to determine how accurate the subtraction of the H₂ line emission falling in the filter band is. Firstly, the flux differences before and after H₂ line subtraction are always ~0.1 mag, except for N3 (~ 0.3 – 0.4 mag). Then, the flux differences between the Brγ and the K_s images are 0.24 ± 0.11 (PISCES 2012) and −0.07 ± 0.13 (SOUL 2020), respectively. The differences are dominated by S1a and N3, which enclose an area of fainter continuum superimposed by strong H2 line emission (compare Fig. 2 and Fig. 1).

This analysis suggests that the derived continuum fluxes may be affected by a systematic epoch-to-epoch error (i.e. constant in each epoch) of up to 0.1–0.2 mag related to the zero point error. In addition, the fluxes from CIAO 2003 are possibly being overestimated due to the slightly non-linear regime of the secondary standard stars (see Sect. E). Finally, a photometric error (i.e. depending on the polygon inside each epoch) of up to ~ 0.1 mag, which is slightly larger for the data points obtained from K-corrected images (UKIRT 2000 and TNG 2006), has to be taken into account as well.

Appendix G Line photometry of the jet with polygons: The low- plus high-spatial-resolution case

We also investigated the flux variability in the H₂ line by adopting the same technique as used for the continuum emission. First, photometry was carried out with polyphot in IRAF by defining a set of larger polygons to allow data points from the low spatial resolution images to be included. The adopted polygons and their labelling are shown in Fig. 2.

As for the photometric errors, one has to take into account an additional source of uncertainty compared to what is discussed in Appendix F, namely how good the subtraction is of the continuum contamination due to dust scattered emission. In fact, two approaches have been followed to estimate the continuum contamination. The first consists in using a narrow-band filter near the H2 passband which does not encompass the H₂ line, namely Brγ and K_cont. The narrow-band images taken into the same runs as H2 exhibit roughly the same spatial resolution as the H2 image and were just scaled to the H2 image using the stars in the field (assuming their spectrum is essentially a continuum), then they were registered and subtracted from the H2 image. This appears as the more accurate correction, particularly when using the Brγ filter, which is very close to the H2 filter passband. One can use the data taken with the TNG in 2001 to compare the correction obtained with Brγ and that obtained with K_cont. The difference between the corrected line fluxes from the various polygons is 0.07 ± 0.09 mag, with the correction with K_cont resulting in fainter magnitudes. This is consistent with the fact that the dust scattered emission is expected to be more intense in the K_cont than in the Brγ band. Nevertheless, the two corrections appear to be equivalent.

The second approach consists in estimating the continuum contamination from the broad-band K_s images. In this case the problem is that the band includes the H₂ 2.12 µm line itself along with other H₂ lines. One has to assume that the continuum flux density is constant inside the band and that the band encompasses only the 2.12 µm line. Then, from the H2 to K_s flux ratios of the stars in the field, one can derive flux scaling factors for the H2 and K_s images and subtract the latter from the former after flux scaling. We used the PISCES 2012 and SOUL 2020 data to compare the continuum corrections obtained from the Brγ and K_s band. In Fig. G.1 we plotted a raw estimate of the continuum contamination in the polygons (given by the difference of the total flux in magnitudes measured in the H2 band and the line flux corrected using the Brγ image) versus the difference (in magnitudes) of the pure line flux from the polygons corrected with the Brγ image and the K_s image. It can be seen that the correction from the K_s image tends to overestimate the continuum contribution compared to that from narrow band images and the difference (in magnitudes) is of the same order of the continuum contamination. So the uncertainty appears to be in the 0. 2 – 0. 5 mag range, increasing with the continuum contamination level.

Fig. G.1 Difference (in magnitudes) between the pure line flux in the polygons after continuum correction from the Brγ image, (H2 - Brγ) and the Ks image, (H2 - Ks) vs the continuum contamination (in magnitudes) in the polygons estimated as the difference between the total flux measured in the H2 band and the line flux corrected using the Brγ image (H2 - Brγ). The symbols refer to the polygons in Fig. 2, as explained at the top of the box. Two datasets have been used (PISCES 2012 and SOUL 2020).

By using the spectra in Caratti o Garatti et al. (2008), we found that the inclusion of various H2 lines in the K_s band cannot contribute more than 0.1 – 0.2 mag, even assuming that the continuum flux density increases as λ⁴ would not affect the continuum correction significantly. This agrees with the results from the PISCES photometry, whereas the SOUL photometry is affected by larger mismatches. We have found that the source of the error lies in the accuracy of the ratio of the filter+atmospheric transmission (narrow band over broad band). An overestimate of this ratio by only 10–20 % accounts for a correction error increasing with an increasing level of contamination as shown in Fig. G.1. As we derived this ratio from the count ratios of the stars in the field, an error of 10–20 % can be easily caused by their spectra containing absorption lines or other irregularities and/or significant differences in the PSFs of the two images (considering only a few stars are available in the AO fields). This highlights the need to check at least visually the level of continuum subtraction when using broad-band images.

Appendix H Short movies

To highlight the time variations in the field of IRAS20126+4104 through the various NIR images available from 2000 to 2020, we have made two short movies by flux-scaling and aligning those images. Figure H.1 shows the flux variations in the continuum K band, obtained from the Brγ filter frames (TNG 2001, CIAO 2003, PISCES 2012, and SOUL 2020) and the line-corrected K frames (UKIRT 2000 and TNG 2006). The resolution was degraded to the worst case and the photometric reference stars in the field were scaled to the same integrated fluxes. To slow down the time variations, we added synthetic images in between each pair of consecutive epoch frames obtained through an epoch-to-epoch linear interpolation of pixel values. No assumptions were made as to possible periodical oscillations of the light curves, so the flux changes were just linearly interpolated between 2006–2012 and 2012–2020. The year is always indicated on the top of the frame as a rough clock marker.

Fig. H.1 (Online movie) Continuum flux variability in the K band.

Fig. H.2 (Online movie) H2 line structure variations between 2003 and 2020.

The movie in Fig. H.2 only uses the high-resolution images (CIAO 2003, PISCES 2012, and SOUL 2020) to show the variations of the H₂ line-emitting structures from 2003 to 2020. The photometric reference stars were scaled to the same integrated fluxes, but the original resolution of the frames has not been modified. The proper motions of knots B are evident, along with the development of bow-shock C1 (HLN13, HLN14) and knot X2 (HLPD2).

Appendix I Kinematic measurements

In Table I.1 we list the proper motions of the jet knots derived in this work. By combining these with the radial velocities measured by Caratti o Garatti et al. (2008), we obtained the 3D velocities listed in Table I.2.

References

All Tables

All Figures

Fig. 1 Three colour image (red H2 filter, green KS, blue Brγ filter) of IRAS 20126+4104 obtained with PISCES and FLAO at the LBT in 2012 (top) and LUCI1 and the new AO system SOUL at the LBT in November 2020 (bottom). The KS image used for the bottom image is that taken on November 4, that is the one exhibiting the best spatial resolution. Approximate positions of radio sources N1 (blue square), N2 (yellow square), and S (orange square, see text) are marked in the top panel with small coloured squares.

In the text

Fig. 2 Jet decomposition in polygons. Top panel: Brγ-filter image of IRAS 20126+4104 obtained with LUCI1 and the AO system SOUL, overlaid with the polygons used for the photometry in the continuum. The adopted designation is labelled. We note that polygon S is composed of S1 and S2, and S1 is further subdivided into S1a and S1b. Bottom panel: pure H2 line emission image of IRAS 20126+4104 obtained with LUCI1 and the AO system SOUL, overlaid with the polygons used for the jet photometry. The adopted designation is labelled. We note that polygon LLS is composed of LLS1 and LLS2, and LLN2 is subdivided into LLN2a and LLN2b.

In the text

	Fig. 3 Pure H2 line emission image of IRAS 20126+4104 obtained with LUCI1 and the AO system SOUL, overlaid with the polygons used for the high-spatial-resolution jet photometry.
In the text

Fig. 4 Photometric variability of the continuum features defined in Fig. 2. We note that N and S indicate the northern and southern lobe, respectively, and PD is the ring-like feature. The numbers increase from south-east to north-west. We note that S is further subdivided into S1 and S2, and S1 is in turn subdivided into S1a and S1b. The symbol size is ~0.2 mag in the vertical direction, which is comparable with the photometric errors. A sinusoid of period ~12 yr is overlaid on the N1 light curve. The same, but with a phase shift of π, is also drawn in the bottom box.

In the text

Fig. 5 Comparison of the continuum emission throughout the various observed epochs. This has been approximated with the Brγ-filter images, except for the runs of July 2000 and August 2006 where broad-band K images corrected for line emission have been used. The flux levels have been adjusted so that the mean of the counts of stars 2 and 4 is the same in all frames. The dimming of the northern lobe (especially region N1) and the simultaneous brightening of the southern lobe in July 2003 are evident.

In the text

Fig. 6 Photometric variability of jet areas in the H2 2.12 µm line emission. We note that LLN and LLS indicate the northern and southern lobe, respectively, and PD is the ring-like feature (see Fig. 2). In addition, LLS is further subdivided into LLS1 and LLS2, and LLN2 is in turn subdivided into LLN2a and LLN2b. The small open triangles (connected with dotted lines) indicate the fluxes with continuum subtraction from the K′ band when corrected following Appendix G. The photometric errors are smaller than the symbol sizes. The magnitudes have been corrected for bandwidth differences.

In the text

	Fig. 7 Photometric variability of jet knots in the H2 2.12 µm line emission. We note that LHN and LHS indicate the northern and southern lobe, respectively, and LHPD is the ring-like feature (see Fig. 3). The magnitudes have been corrected to the LUCI1 H2 filter bandwidth.
In the text

Fig. 8 Knot identification in the SOUL H2 continuum-subtracted image of the IRAS 20126-4104 flow and knot proper motions. Panel a: zoom-in on the A and X knot regions (red-shifted lobe). Contour levels at 5, 7, 10, 15, 20, and 25σ are overlaid. We note that 1σ corresponds to 8×10−23 W cm−2 arcsec−2. Panel b: zoom-in on the B and C1 knot regions (blue-shifted lobe). Contour levels at 5, 7, 10, 15, 30, and 40σ are overlaid. Panel c: zoom-in on the C1 and C2 knot regions (blue-shifted lobe). Contour levels at 5, 10, 20, 30, 40, 50, 60, 80, and 120σ are overlaid. Knot peaks are indicated by white dots. Panel d: overall view of the IRAS 20126+4104 flow close to the source. The main structures as reported in Caratti o Garatti et al. (2008) are labelled. Panel e: proper motions (PMs) with their uncertainties (blue arrows and red ellipses) in 100 yr of structures and sub-structures along the H2 jet in IRAS20126+4104. The actual observed shifts are approximately one fourth the length of the corresponding arrow. The red circle marks the position (along with its uncertainty) of the protostellar continuum emission at 1.4 mm (Cesaroni et al. 2014). The main structures as labelled in Caratti o Garatti et al. (2008) are also indicated.

In the text

	Fig. 9 Average tangential velocity vs average projected distance to the source for each group of knots.
In the text

	Fig. 10 Evolution of the bright bow shock in the C1h and C1k regions. The top, middle and bottom panels show the same area of the sky at the three different epochs (2003, 2012, and 2020). The fast moving flow (υtg ~ 100–120 km s−1) corresponding to B3 shocks against the C1h and C1k regions (υtg ~ 10–70 km s−1) producing a bright bow shock. The PMs of knots are shown as in Fig. 8.
In the text

	Fig. 11 Tangential velocity (bottom panel), position angle (middle panel), and dynamical age (top panel) vs. distance to the source for each knot. Red labels show knots in the C1 group that might belong to a different flow (see text).
In the text

	Fig. 12 Average dynamical age vs. average distance to the source for each group of knots.
In the text

Fig. 13 Mean directions of the proper motion vectors in the blue-shifted lobe (top panel) and in the red-shifted lobe (bottom panel). The proper motion track of the protostar is marked by the red line and its current location by the red cross. Clearly, the more distant the knot group is from the protostar, the more distant (to the north-east) the intersection of its mean direction is with the track. The protostar is moving towards the south-west.

In the text

Fig. 14 Ballistic bow-shock model of Ostriker et al. (2001) compared with knot chains A1 (top row), B2 (middle row), and C1 (bottom row). In the column on the left 6″ × 6″ zoom-ins of the SOUL H2 image are shown with the bow-shock model shape overlaid. In the middle column, the transverse velocities of the involved knots are compared with the mean shell (dotted line) and outer surface layer (solid line) values predicted by the model. The same comparison is made in the column on the right for the longitudinal velocities. The knots whose speed departs from the model predictions the most are marked in red and labelled. For the sake of simplicity, the outflow is assumed to lie on the plane of the sky.

In the text

	Fig. 15 Zoom-in on Fig. 8 around the source showing H2O maser proper motions (magenta arrows) as measured by Moscadelli et al. (2011) and H2 proper motions (blue arrows) from this paper. Red ellipses mark the uncertainties on PMs. The black circle shows the position of the protostellar continuum at 1.4 mm from Cesaroni et al. (2014).
In the text

	Fig. A.1 Profile of the PSF of the AO guide star along the right ascension axis (black line) in the SOUL H2 image overlaid with the best-fitting Moffat distribution (blue line). The plate scale is 0″.015 pix−1.
In the text

	Fig. A.2 Maximum and minimum FWHMs obtained from the Moffat fits to the PSFs as a function of distance from the AO guide star in the SOUL (red triangles) and FLAO (blue triangles) H2 images.
In the text

	Fig. C.1 Brγ-filter image of IRAS 20126+4104 obtained with LUCI1 and the AO system SOUL, with the identification number of the stars available as photometric references.
In the text

	Fig. E.1 Photometric variability of the stars in the SOUL FoV. The star ID (with reference to Fig. C.1) is indicated on the right side of the panel. The filters are indicated with symbols as explained at the top of the panel. The data taken on JD 2452830.5 (CIAO) and JD 2456099.8 (PISCES) have not been calibrated and the magnitudes were computed by setting star 4 to 11.977 mag in all bands.
In the text

Fig. G.1 Difference (in magnitudes) between the pure line flux in the polygons after continuum correction from the Brγ image, (H2 - Brγ) and the Ks image, (H2 - Ks) vs the continuum contamination (in magnitudes) in the polygons estimated as the difference between the total flux measured in the H2 band and the line flux corrected using the Brγ image (H2 - Brγ). The symbols refer to the polygons in Fig. 2, as explained at the top of the box. Two datasets have been used (PISCES 2012 and SOUL 2020).

In the text

	Fig. H.1 (Online movie) Continuum flux variability in the K band.
In the text

	Fig. H.2 (Online movie) H2 line structure variations between 2003 and 2020.
In the text