The SOUL view of IRAS 20126+4104. Kinematics and variability of the H 2 jet from a massive protostar

Context. We exploit the increased sensitivity of the recently installed adaptive optics SOUL at the LBT to obtain new high-spatial-resolution near-infrared images of the massive young stellar object IRAS20126 + 4104 and its outﬂow. Aims. We aim to derive the jet proper motions and kinematics, as well as to study its photometric variability by combining the novel performances of SOUL together with previous near-infrared images. Methods. We used both broad-band ( K s , K (cid:48) ) and narrow-band (Br γ , H2) observations from a number of near-infrared cameras (UKIRT / UFTI, SUBARU / CIAO, TNG / NICS, LBT / PISCES, and LBT / LUCI1) to derive maps of the continuum and the H 2 emission in the 2 . 12 µ m line. Three sets of images, obtained with adaptive optics (AO) systems (CIAO, in 2003; FLAO, in 2012; SOUL, in 2020), allowed us to derive the proper motions of a large number of H 2 knots along the jet. Photometry from all images was used to study the jet variability. Results. We derived knot proper motions in the range of 1 . 7–20 . 3 mas yr − 1 (i.e. 13–158 km s − 1 at 1 . 64 kpc), implying an average outﬂow tangential velocity of ∼ 80 km s − 1 . The derived knot dynamical age spans a ∼ 200–4000 yr interval. A ring-like H 2 feature near the protostar location exhibits peculiar kinematics and may represent the outcome of a wide-angle wind impinging on the outﬂow cavity. Both H 2 geometry and velocities agree with those inferred from proper motions of the H 2 O masers, located at a smaller distance from the protostar. Although the total H 2 line emission from the knots does not exhibit time variations at a (cid:101) > 0 . 3 mag level, we have found a clear continuum ﬂux variation (radiation scattered by the dust in the cavity opened by the jet) which is anti-correlated between the blue-shifted and red-shifted lobes and may be periodic (with a period of ∼ 12 − 18 yr). We suggest that the continuum variability might be related to inner-disc oscillations which have also caused the jet precession.


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 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 4 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. 1999Cesaroni et al. , 2005Cesaroni et al. , 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(Hofner et al. , 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. (2005Sridharan et al. ( , 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.05kpc) 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.19  −0.092 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.05kpc.
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 2 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.

Observations and data reduction
2.1.SOUL AO at the LBT SOUL (Pinna et al. 2016(Pinna et al. , 2021) )  Notes. (a) H2 and Brγ were obtained using the detector reading mode MER and all single DITs in a set of NDIT exposures are saved as a 3D FITS file, so that each frame consists of a data-cube of NDIT planes. (b) Taken with an open loop.
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.49star (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.

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 −1 , 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 1 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 A113, page 3 of 26 A&A 672, A113 (2023) 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.
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.

Continuum subtraction
To obtain H 2 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 2 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).

Jet photometry
We performed photometry of the jet emission, both in the continuum and in the H 2 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 2 line emission and where Brγ emission has not been detected.Only in two cases did we use K images (UKIRT 2000 andTNG 2006) after subtracting the derived H 2 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 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 G, 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 2 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 smallerscale 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.

Proper motion analysis
We used the H2 continuum-subtracted high-angular resolution images from CIAO, PISCES and SOUL after flux calibration 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 H 2 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.
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 crosscorrelation 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 A113, page 5 of 26 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.

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 2 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.24 pc), 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).

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.

H 2 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 A113, page 6 of 26 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.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 bandwidths 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: (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.8 mag 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 −1 (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.I.1.In particular, the table provides knot IDs, coordinates (as derived from the SOUL map), proper motions (in mas yr −1 ), 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 −1 and groups X and A2 on the red-shifted side with average tg ∼ 40 and 70 km s −1 ) 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 −1 , 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 −1 (see e. g. knots B3c and B3bb) and shocks against the C1hi, C1hj, and C1hk regions ( tg ∼ 10−70 km s −1 , corresponding to polygons HLN13 and HLN14 in Fig. 3), producing a bow shock (knots C1hi, hk, A113, page 7 of 26 The average tangential velocity of the H 2 flow is 80±30 km s −1 .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 −1 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 2 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 τ (τ) 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 (τ X = 280±50 yr), whereas the τ value of the farthest group towards the north-west (C2) is roughly consistent with that of A1, namely τC2 = 1070±300 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 τ 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.

H 2 3D kinematics
By combining the tangential velocities ( tg ) derived in this work with the H 2 radial velocities ( r ) obtained at high-spectral A113, page 8 of 26  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: (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 2 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).

Bow shocks
We tested our data to check that the H 2 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 2 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 −6 W m −2 sr −1 ) roughly agree with a modelled bow shock impacting a medium of density ∼10 5 cm −3 with a speed of 50−60 km s −1 .Since the actual speed that we have derived is ∼45 km s −1 , A1 is likely to represent a terminal shock.The width (∼800 au) and average brightness (∼9 × 10 −7 W m −2 sr −1 ) of B3c again points to a bow shock propagating in a medium of density ∼10 5 cm −3 with a speed of 50−60 km s −1 .Finally, C1h-C1k has a width of ∼1600 au and an average brightness of ∼0.4−1 × 10 −5 W m −2 sr −1 , which is consistent with a bow shock propagating in a medium of density ∼10 5 −10 6 cm −3 with a speed of 50−60 km s −1 .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 −1 ).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 A1, 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 2 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 −1 , a sound speed C s = 8 km s −1 (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 A1, B2, and C1.The results of the fit are also listed in Table 2. Figure 14 shows that the morphology of knots A1, B2, and C1 is well fit by the bow-shock outer-shell model.The proper motions of the knots making up group A1 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.

Jet precession
A comparison between the continuum emission and the H 2 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 2 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 A113, page 13 of 26  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 −1 with respect to the systemic velocity (Palau et al. 2017), hence it is consistent with the radial velocities derived from the H 2 emission by Caratti o Garatti et al. (2008).A nearly symmetrical patch of SO emission between 11 and 32 km s −1 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 nonnegligible 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 −1 towards A, but still above ∼100 km s −1 towards C (see Table I.2).
The global morphology of the H 2 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 −1 , 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 2 luminosity monotonically increases with time, but apparently only less than ∼0.01 mag yr −1 , 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.

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 2 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 wideangle 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 2 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 subarcsecond resolution using the PdB interferometer.Notably, one of the detections spatially matches the H 2 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 wideangle wind that is produced by episodic accretion and ejection.This would also be consistent with knot fragmentation seen in A113, page 15 of 26 the IRAS 20126+4104 jet.Finally, by fitting an ellipse to the ringshaped 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.

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 −1 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 −1 , 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 −1 .
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.

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 H 2 and CO emission display a more extended north-south outflow centred on the protostar position (see Shepherd et al. 2000).The southern lobe is redshifted, has a jaw-like morphology, and the south-eastern lobe of the HCO + and H 2 outflow lies at its south-eastern border, indicating, on a larger scale, a different precessing geometry with respect to the closest H 2 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 2 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 2 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 A113, page 16 of 26 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 where b is the smallest encounter distance.These authors computed Γ 0 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 H 2 inner jet (the average dynamical ages of groups A1 and C2 are 2200±900 yr 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 −1 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 −1 (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.

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 2 jet and outflow.Two other observing runs with AO (2003 and2012) 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 −1 , which translate into tg = 13-158 km s −1 at a distance of 1.64 kpc.The average tg of the flow is 80±30 km s −1 .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 2 knots, we have reconstructed the 3D kinematics of the jet, which roughly expands along an axis with an average inclination angle of 8 -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 2.12 µm luminosity does not exhibit time variations at a > ∼ 0.4 mag level, which is consistent with hydrodynamical 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-spatialresolution 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.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).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 The zero magnitude fluxes from 2Mass can be used to define a source which has mag i = 0 in every NIR band linked to this standard system.By simply fitting to the 2Mass J, H, K s zero magnitude fluxes, one obtains log(A 0 ) = −12.16and γ 0 = −3.59,with F λ in W cm −2 µm −1 , 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 × log[F λ,0 (λ j )∆λ j ] (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 mag ′ j = 2.5 × log[F λ,0 (λ j )/F λ (λ j )] + 2.5 log[∆λ j /∆λ] (D.5) mag ′ j − mag j = 2.5 log[∆λ j /∆λ], (D.6) where mag j is the magnitude that would be measured in the exact photometric system of the standard star:

Fig. 1 .
Fig. 1.Three colour image (red H2 filter, green K S , 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 K S 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.

Fig. 2 .
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 H 2 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 .
Fig. 3. Pure H 2 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.
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 −1 .The results are listed in TableI.1.In particular, the table provides knot IDs, coordinates (as derived from the SOUL map), proper motions (in mas yr −1 ), 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 −1 and groups X and A2 on the red-shifted side with average tg ∼ 40 and 70 km s −1 ) 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 −1 , 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 −1 (see e. g. knots B3c and B3bb) and shocks against the C1hi, C1hj, and C1hk regions ( tg ∼ 10−70 km s −1 , corresponding to polygons HLN13 and HLN14 in Fig.3), producing a bow shock (knots C1hi, hk, 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.

Fig. 6 .
Fig.6.Photometric variability of jet areas in the H 2 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.

Fig. 7 .Fig. 8 .Fig. 9 .
Fig. 7. Photometric variability of jet knots in the H 2 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. 10 .
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. 12 .
Fig. 12.Average dynamical age vs. average distance to the source for each group of knots.

Fig. 13 .
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 protostar, the more distant (to the north-east) intersection of its mean direction is with the track.The protostar is moving towards the south-west.Table2.Shape of three possible bow shocks (from fitting Eq. (22) ofOstriker et al. 2001).
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 .
Fig. 15.Zoom-in on Fig. 8 around the source showing H 2 O maser proper motions (magenta arrows) as measured by Moscadelli et al. (2011) and H 2 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).

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 .

Fig
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.
mag j = 2.5 × log[F λ,0 (λ j )/F λ (λ j )].(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 A113, page 20 of 26
, 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; Table I.1.Derived proper motions, kinematics, and dynamical ages of the IRAS 20126+4104 knots.