© The Authors 2023

1 Introduction

In recent years, observations have provided us with increasing evidence of circumstellar rotating structures around B-type (proto)stars, especially since the advent of the Atacama Large Millimeter/submillrmeter Array (ALMA; see e.g. the review by Beltrán & de Wit 2016). This strongly suggests that disk-mediated accretion could be a viable mechanism to feed even the most massive stars. However, we still do not understand the physical properties of these rotating structures, nor how such an accretion proceeds, whether through a smooth, continuous flow or episodically with parcels of material falling onto the star. The recent detection of outbursts (Stecklum 2016, 2021; Caratti o Garatti et al. 2017; Hunter et al. 2017, 2021; Burns et al. 2020, 2023; Chen et al. 2021) in a few luminous young stellar objects (YSOs) >8 M_⊙ provides us with the intriguing possibility that these phenomena could be the consequence of episodic accretion events, akin to those commonly observed in low-mass YSOs as FU Orionis and EX Orionis events (Audard et al. 2014; Fischer et al. 2023).

In this study, we focus on the outburst from the massive (proto)star S255IR NIRS 3, located at a distance of ${1.78}_{-0.11}^{+0.12}$ kpc (Burns et al. 2016). This object is unique because the outburst has been observed not only in the emission of some maser species (Fujisawa et al. 2015; Moscadelli et al. 2017; Hirota et al. 2021) and in lines and continua at IR and (sub-)millimetre wavelengths (Caratti o Garatti et al. 2017; Uchiyama et al. 2020; Liu et al. 2018), but also in the centimetre domain (Cesaroni et al. 2018; hereafter Paper I). A time delay of ~1 yr was found between the onset of the IR and radio outbursts, consistent with the different mechanisms at the origin of the two: in fact, the IR outburst is based on radiative processes propagating at velocities comparable to the speed of light, whereas the radio outburst is due to shocks expanding approximately at the speed of sound. In Paper I, we present compelling evidence of an exponential increase in the radio emission from the thermal jet associated with this source. Given the existence of a disk-jet system in S255IR NIRS 3 (Boley et al. 2013; Wang et al. 2011; Zinchenko et al. 2015; Liu et al. 2020), we believe that we are witnessing an episodic accretion event mediated by the disk, where part of the infalling material has been diverted into the associated jet Fedriani et al. (2023). In Paper I, we show that a simple model of an expanding jet can satisfactorily reproduce the increase in the radio flux observed in four bands. Although the emission was basically unresolved, we predicted that in a few years the jet expansion should make it possible to resolve its structure. With this in mind, we performed both monitoring of the radio emission and sub-arcsecond millimetre imaging at two epochs, several years after the beginning of the outburst. In this article we report on the results of these observations.

2 Observations

The radio emission was monitored with the Karl G. Jansky Very Large Array (VLA) and, at a later time, with ALMA. In the following we describe the observational setup and data reduction separately for the two datasets. In both cases, for the phase centre we chose the position α(J2000)=06^h 12^m 54^s.02, δ(J2000)=17° 59′ 23″.1.

2.1 Very large array

S255IR NIRS 3 was observed at 11 epochs from January 2017 to January 2018 (project codes: 16B-427 and 17B-045). The observing dates and array configurations are listed in Table 1, where for the sake of completeness we have also included the 2016 data already presented in Paper I. We note that the table contains as yet unpublished data in the L band obtained in the same observing run (project 16A-424) described in Paper I.

The signal was recorded with the Wideband Interferometric Digital ARchitecture (WIDAR) correlator in six bands, centred approximately at 1.5 (L band), 3 (S), 6 (C), 10 (X), 22.2 (K), and 45.5 GHz (Q), in dual polarisation mode. The total observing bandwidth (per polarisation) was 1 GHz in the L band, 2 GHz in the S band, 4 GHz in the C and X bands, and 8 GHz in the K and Q bands. The primary flux calibrator was 3C48 and the phase-calibrators were J0632+1022 in the L and S bands, J0559+2353 in the C and X bands, and J0539+1433 in the K and Q4 bands.

We made use of the calibrated dataset provided by the NRAO pipeline and subsequent inspection of the data and imaging were performed with the CASA¹ package, version 5.6.2-2. The continuum images were constructed using natural weighting to maximise flux recovery. Typical values of the 1σ RMS noise level and synthesised beam are given in Table 2 for each band and array configuration.

The flux density of the compact, variable source centred on S255IR NIRS 3 has been estimated inside a polygonal shape encompassing the compact, unresolved radio source and is given in Table 1 for each band and epoch. As explained in Paper I, in our VLA observations the continuum emission from S255IR NIRS 3 appears as an unresolved component (the variable source of interest for our study) plus two large-scale lobes separated by ~8″, or ~ 15 000 au, in projection (see Fig. 1 of Paper I) that in all likelihood are originating from a previous outburst that occurred many decades ago. At the longest wavelengths it was possible to resolve the central source from the lobes only in the most extended array configuration. In the other cases, we could only measure the total flux density from all components and the corresponding values are reported in Table 1 as upper limits. In a few cases, as indicated in the footnotes of the table, we attempted a correction to the flux, under the assumption that the flux of the lobes is constant in time and thus equal to the value measured (at a different epoch) with the A array. With this approach a caveat is in order, because a compact configuration may be sensitive to lobe structures that are resolved out in a more extended configuration. This implies that the corrected flux densities could be still an upper limit.

In order to estimate the uncertainty on the flux density measurements of S255IR NIRS 3, we took advantage of the presence of a compact, marginally resolved continuum source located approximately at α(J2000) = 06^h 12^m53^s.61, δ(J2000) = 18° 00′ 26″.4, which happens to fall in the primary beam in all bands except band Q. Under the reasonable assumption that this object is not variable, by comparing the flux density measurements obtained at different epochs at the same frequency, we estimated a relative error of 10% in all bands.

2.2 Atacama Large Millimeter/submillimeter Array

S255IR NIRS 3 was observed with ALMA in band 3 on June 6, 2019, and September 3, 2021 (project 2018.1.00864.S, P.I. R. Cesaroni). The main characteristics of the observations are summarised in Table 3. The correlator was configured with 4 units of 2 GHz, in double polarisation, centred at 85.2, 87.2, 97.2, and 99.2 GHz. The spectral resolution is 0.49 MHz (corresponding to ~ 1.5–1.7 km s⁻¹, depending on the frequency), sufficient to identify line-free channels and obtain a measurement of the continuum emission.

The data were calibrated through the ALMA data reduction pipeline. For each 2 GHz correlator unit, we created a data cube using task tclean of CASA, adopting natural weighting and a circular beam of 0″.097 for the first epoch and 0″ 087 for the second one. To create a continuum map for each 2 GHz band, we used the STATCONT software² developed by Sánchez-Monge et al. (2018). In this way we also obtained cubes of the continuum-subtracted line emission. Finally, the four continuum maps were averaged together to increase the signal-to-noise ratio. The measured fluxes at the two epochs are reported in Table 1.

We note that the derivation of the continuum images described above also provides us with continuum-subtracted channel maps. Although a study of the line emission in this region goes beyond the purposes of this article, in the following we briefly consider the maps of two molecular transitions, SO(2₂−1₁) at 86093.983 MHz and H¹³CN(1−0) at 86338.7 MHz.

3 Results

3.1 Continuum emission at λ=3 mm

The structure of the radio jet from S255IR NIRS 3 at the two epochs observed with ALMA is shown in Fig. 1. The most striking result is the clear expansion of the jet, which becomes significantly more elongated to both the NE and SW. At first look, the figure seems to outline a bipolar structure, where the star might be located close to the geometrical centre, between the two jet lobes. However, this is not the case. Careful comparison between Figs. 1a and b reveals that, while the NE peak is moving away from the centre, the SW peak stays still, consistent with this being the location of the (proto)star. This result is confirmed by the coincidence of the SW peak with the peak of the submillimetre emission (white cross) measured by Liu et al. (2020), as expected if the (proto)star lies inside the parental molecular, dusty core. Figure 2 clearly confirms this scenario by showing an overlay of the molecular emission maps of the SO(2₂−1₁) and H¹³CN(1−0) lines observed by us, with an image of the 3 mm continuum emission. It is worth noting that in both transitions a dip is seen towards the continuum peak, which suggests that part of the line emission is likely absorbed against the bright continuum. This is consistent with the high brightness temperature of both 3 mm continuum peaks, of the order of 500 K (obtained from the maps cleaned with uniform weighting).

With all the above in mind, in the following we assume that the (proto)star powering the radio jet is located at the SW peak of the 3 mm continuum emission, namely at α(J2000) = 06^h 12^m 54^s.012, δ(J2000) = 17° 59′ 23″.04. We prefer to use the peak position of our maps instead of that of the sub-millimetre maps of Liu et al. (2020), because of the higher angular resolution (0″.087 instead of 0″.14).

It is also interesting to compare the jet lobes observed in our ALMA maps with those seen on a much larger scale (see Fig. 1 of Paper I). The comparison is presented in Fig. 3, where we show the maps of the 3.6 cm and 3 mm continuum emission. Clearly, the directions of the symmetry axes of the two pairs of lobes are remarkably different, with position angle (PA) of ~70°, on the large scale, and ~48°, on the small scale. The fact that the two axes have different orientations and do not intersect at the position of the star (i.e. at the SW peak of the 3 mm continuum) can be interpreted in two ways: either we are dealing with two jets originating from two different YSOs, or the jet is precessing and the star is moving on the plane of the sky at a different speed with respect to the large-scale lobes. The latter scenario implies that either the star or the ejected material is experiencing deceleration or acceleration, such that one of the two is lagging behind the other.

Although it is impossible to rule out one of the two hypotheses, in the following we assume that all lobes belong to the same jet, because there is only one core lying along the jet axis and evidence for precession has been found in S255IR NIRS 3 (Wang et al. 2011; Fedriani et al. 2023), as well as other similar objects (see e.g. Shepherd et al. 2000, Cesaroni et al. 2005, Sánchez-Monge et al. 2014, Beltrán et al. 2016). Further support for this hypothesis is given by the progressive change of the position angle of the jet from the large to the small scale, as shown in Table 4. This trend is indeed consistent with a jet outflow undergoing precession.

Fig. 1 Maps of the 3 mm continuum emission obtained with ALMA. (a) Data acquired on June 6, 2019. Contour levels range from 0.52 to 16.12 in steps of 2.6 mJy beam−1. The cross marks the peak of the 900 μm continuum emission imaged by Liu et al. (2020). The circle in the bottom left represents the synthesised beam. (b) Same as the top panel, but for the map obtained on September 3, 2021. Contour levels range from 0.29 to 16.24 in steps of 1.45 mJy beam−1.

3.2 Continuum emission at λ ≥ 7 mm

In Fig. 4, we present the continuum spectra obtained from the flux densities in Table 1. Our VLA observations of S255IR NIRS 3 span an interval of time prior to the ALMA observations (see Table 1) and the radio emission of the variable source in S255IR NIRS 3 is basically unresolved in all of our VLA observations. However, we can study the spatial evolution of the jet by determining the position of the peak at different times. For this purpose we fitted a 2D Gaussian to the K-band maps with sub-arcsecond resolution. We prefer the 1.3 cm data to the 7 mm data, which would provide us with better resolution, because in the K band the S/N is higher and in the Q band contamination by dust thermal emission might be present. In Fig. 5, we plot the distribution of the peak positions thus obtained. To give an idea of the uncertainty on these positions, we also draw ellipses corresponding to one-fifth of the synthesised beams. For our analysis we also included the NE peak of our ALMA data and that of Obonyo et al. (2021; hereafter OLHKP). Despite the large uncertainties, it is clear that the distance of the peak from the star is increasing with time, as expected if the jet is expanding. This expansion appears to slow down with time, because the mean velocity (projected on the plane of the sky) estimated from the ratio between the separation of the peaks at the last two epochs (ALMA data) and the corresponding time interval is ~40 au/820 days ≃ 84 km s⁻¹, much less than the mean velocity from the beginning of the radio burst, namely ~472 au/1881 days ≃ 436 km s⁻¹, where ~472 au is the distance of the NE peak from the star at the last epoch. Using the same approach, Fedriani et al. (2023) estimated an expansion speed of 450 ± 50 km s⁻¹ from their IR data. We further analyse the expansion of the jet in Sect. 4.3.1.

The flux density of S255IR NIRS 3 is changing with time in all bands, as shown in Fig. 6. One can identify two phases: the first when the flux increases exponentially for ~200 days; the second when the flux remains basically constant, or slightly declines towards the end of our monitoring. This behaviour is the same at all frequencies ≥6 GHz, but is not so obvious at the two longest wavelengths, where a precise estimate of the flux density of the compact variable source is not possible for the reason explained in Sect. 2.1. Moreover, in these bands our monitoring is more limited in time. It is hence quite possible that the flux density below 6 GHz has a different behaviour than that at higher frequencies. Therefore, the radio emission at 1.5 and 3 GHz, and probably also part of the 6 GHz flux, might not be due to free-free radiation but to another mechanism, such as synchrotron emission, which has been detected towards a number of extended radio jets from YSOs (e.g. Carraco-González et al. 2010; Moscadelli et al. 2013; Brogan et al. 2018; Sanna et al. 2019 and references therein). The evident change of slope in some of the VLA spectra seems to support this possibility. In this respect, the most representative is the spectrum acquired on 2016/10/15 (see Fig. 4), where the 1.5 GHz point lies well above any plausible extrapolation of the other fluxes. For all these reasons, in the following we focus our study on the emission above 6 GHz, with the caveat that even the 6 GHz flux might be partly contaminated by a non-thermal contribution, as suggested by OLHKP. It is worth noting that the existence of synchrotron emission from S255IR NIRS 3 is supported by the recent possible detection of high-energy gamma-ray emission from this source (see de Ona Wilhelmi et al. 2023).

Fig. 2 Maps of the line emission obtained with ALMA. (a) Map of the emission averaged over the H13CN(1−0) line (white contours) overlaid on the continuum image at 3 mm obtained on June 6, 2019. Contour levels range from 1.2 to 4.56 in steps of 0.48 mJy beam−1. The circle in the bottom left represents the synthesised beam. (b) Same as the top panel, but for the SO(22−11) line. Contour levels range from 1.65 to 4.62 in steps of 0.33 mJy beam−1.

Fig. 3 Maps of the 3.6 cm continuum emission obtained with the VLA on July 10, 2016 (red) and of the 3 mm continuum emission (blue), also shown in Fig. 1b. The dashed lines denote the symmetry axes of the bipolar structures. Red contour levels range from 0.1 to 0.6 in steps of 0.1 mJy beam−1. Blue contour levels range from 0.29 to 8.99 in steps of 2.9 mJy beam−1.

4 Analysis and discussion

4.1 Nature of the 3 mm continuum emission

For the reasons presented in Sect. 3.1, we have concluded that the (proto)star should lie at the SW peak of the 3 mm continuum map. If this is the case, one expects this peak to have a significant contribution from thermal dust emission, whereas the NE peak, being part of a jet lobe, should be dominated by free-free emission. To investigate these assumptions, we computed the spectral index of the 3 mm continuum over the maps in Fig. 1. For this purpose, we used the maps obtained from the four correlator units centred at 85.2, 87.2, 97.2, and 99.2 GHz (see Sect. 2.2), created with the same clean beam. For each pixel of the maps, the spectral index was computed from a least-square fit to the four fluxes in a log S_v−log v plot. The fit was performed only in those pixels where all of the four fluxes were above the 5σ level. The result is shown in Fig. 7, where the formal error obtained from the fit ranges from 0.03 towards the emission peaks to 0.3 towards the borders.

At 3 mm it is reasonable to assume that dust emission is optically thin and the flux density is ∝ v^γ, with γ = 2−4, because the dust absorption coefficient is believed to vary as v^β with β = 0−2 (see e.g. D’Alessio et al. 2001; Sadavoy et al. 2013), where β = 0 corresponds to the case of large grains (‘pebbles’) in disks (Testi et al. 2014). The same assumption also holds for the free-free emission: the spectra in Fig. 4 flatten beyond 45 GHz, and hence the flux density is ∝ v^−0.1. Therefore, more positive spectral indices can be associated with dust emission and, conversely, more negative ones with free-free emission. At both epochs dust emission arises from the region around the SW peak, as expected for a deeply embedded star, while the NE lobe of the jet is characterised by free-free emission. Noticeably, some free-free emission is also detected towards the most south-western tip.

When estimating the flux density at wavelengths of 3 mm or shorter, it is thus necessary to distinguish between the NE lobe and the rest of the source. In Table 5 we give all these fluxes for the two epochs of the ALMA observations. It is worth pointing out that the spectral index over the SW region is mostly ≲2, whereas the typical index expected for pure dust emission, should lie approximately between 2 and 4. This suggests the presence of non-negligible free-free emission around the SW peak as well. It is possible to estimate what fraction of the total flux is due to free-free as follows. The total flux can be written as $S_v=S_v^{\text{d}}+S_v^{\text{ff}}$, where $S_v^{\text{ff}}\propto v^{-0.1}$ is the optically thin free-free flux and $S_v^{\text{d}}\propto v^{\gamma }$, with γ = 2−4, is the dust flux. As previously explained, the spectral index between v₁ = 85.2 GHz and v₂ = 99.2 GHz was estimated assuming S_v ∝ v^α, and hence we have $$\frac{S_{v_2}}{S_{v_1}}={\left(\frac{v_2}{v_1}\right)}^{\alpha }$$(1) $$\frac{S_{v_2}}{S_{v_1}}=\frac{S_{v_2}^{\text{d}}+S_{v_2}^{\text{ff}}}{S_{v_1}^{\text{d}}+S_{v_1}^{\text{ff}}}=\frac{S_{v_1}^{\text{d}}{\left(\frac{v_2}{v_1}\right)}^{\gamma }+S_{v_1}^{\text{ff}}{\left(\frac{v_2}{v_1}\right)}^{-0.1}}{S_{v_1}^{\text{d}}+S_{v_1}^{\text{ff}}}.$$(2)

After some algebra, we obtain the ratio $$R_{\text{S}}\equiv \frac{S_{v_1}^{\text{ff}}}{S_{v_1}^{\text{d}}}=\frac{{\left(\frac{v_2}{v_1}\right)}^{\gamma }-{\left(\frac{v_2}{v_1}\right)}^{\alpha }}{{\left(\frac{v_2}{v_1}\right)}^{\alpha }-{\left(\frac{v_2}{v_1}\right)}^{-0.1}}$$(3)

and, from this, the fraction of the total flux due to free-free emission, R_S/(1 + R_S). Using the values of α in Fig. 7, we find that such a fraction on average ranges from 65%, for γ=2, to 85%, for γ=4. This implies that 15–20 mJy out of 23 mJy emitted by the SW component (see Table 5), are contributed by free-free emission, also consistent with the brightness temperature of ~500 K measured at 3 mm towards the SW peak, probably too large to be due only to dust emission. We note that part of this free-free emission might be due to ionisation by the embedded star of ~20 M_⊙ (Zinchenko et al. 2015).

As already mentioned, some free-free emission is also detected towards the tip of the SW region and becomes more prominent and more extended with time, as one can see by comparing Figs. 1a to b. This hints at the existence of another jet lobe emerging from the dusty core and expanding towards the SW.

To shed light on the nature of this putative lobe and set a constraint on its age, in Fig. 8 we compare the jet structure observed by OLHKP at 1.3 cm on May 5, 2018, with our ALMA image at 3 mm. It seems that during the period of our monitoring, up to the observation of OLHKP (664 days after the onset of the radio outburst), no significant free-free emission was seen to the SW of the (proto)star at any of the wavelengths observed with the VLA. We thus conclude that in all likelihood the SW jet lobe appeared only recently, between May 2018 and June 2019. Therefore, the radio flux variations monitored by us until January 2018 are to be attributed only to the NE lobe.

Fig. 4 Continuum spectra of the variable, compact component in S255IR NIRS 3. The date of the observations is given in the top left of each panel. The vertical error bars assume a calibration error of 10% in all bands, while the horizontal bars indicate the bandwidth covered at each frequency. The upper limits denote that the flux density is also contributed by the emission from the large-scale lobes, due to the limited angular resolution. Red symbols mark data in the L (1.5 GHz) and S (3 GHz) bands. The curves are the best fits obtained with the model of a thermal jet described in Sect. 4.3.2.

Fig. 5 Distribution of the peaks of the 3 mm and 1.3 cm continuum maps obtained, respectively, with ALMA and with the A and B configurations of the VLA. The circles indicate our VLA data, the triangle the data of OLHKP, and the two squares our ALMA data. The dashed line connects the points in order of time (the first point in the bottom-right corner corresponds to the observation on July 10, 2016). The third point from the top left is obtained from the OLHKP data. The cross marks the position of the (proto)star, and the solid red line is the axis of the radio jet (PA = 48°). The dotted ellipses are the synthesised beams scaled down by a factor of 5.

Fig. 6 Flux density of the variable, compact source in S255IR NIRS 3 as a function of time from the beginning of the radio burst (July 10, 2016) in all bands observed with the VLA. The two rightmost points are the measurements at 6 and 22 GHz by OLHKP. Curves are colour-coded by frequency (in gigahertz). Arrows denote points with upper limits (see Sect. 2.1), which are connected by dashed lines. The two black arrows pointing to the right indicate the flux densities at 92 GHz of the NE lobe (see Table 1) measured 1061 and 1881 days after the radio burst. Typical error bars are shown in the bottom right of the figure. The vertical dotted lines mark the beginning and end of an array configuration, as labelled above the figure.

Fig. 7 Maps of the 3 mm continuum emission (contours) overlaid on the corresponding map of the spectral index at 3 mm. (a) Data obtained with ALMA on June 6, 2019. Black contour levels range from 10% to 90% in steps of 10% of the peak emission. White contour levels correspond to a spectral index equal to 0. The cross marks the position of the (proto)star. The circle in the bottom left represents the synthesised beam. (b) Same as the top panel, but for the map obtained on September 3, 2021.

Fig. 8 Map of the 3 mm continuum emission obtained with ALMA on September 3, 2021, overlaid on the 1.3 cm image obtained by OLHKP on May 5, 2018. Contour levels are the same as in Fig. 1b. The cross marks the position of the (proto)star. The ellipse in the bottom left represents the synthesised beam at 1.3 cm.

4.2 Origin of the SW lobe

One may wonder if the emerging SW lobe corresponds to a new radio burst or is somehow related to the accretion outburst observed by Caratti o Garatti et al. (2017). Only follow-up observations of the jet structure will allow us to establish if the new lobe is as prominent and long-lasting as the NE one. However, we point out that the IR monitoring performed by Uchiyama et al. (2020) and Fedriani et al. (2023) between November 2015 and February 2022 as well as the methanol maser observations³ of the Maser Monitoring Organization (M2O)⁴ have not revealed any other burst after that of Caratti o Garatti et al. (2017) and before the ejection of the SW lobe. We conclude that both jet lobes could arise from the same accretion event, although with a time lag between them of 22–35 months.

This hypothesis is also supported by a noticeable feature of the jet system in S255IR NIRS 3, namely that the extension of the NE lobe is about twice as much as that of the SW lobe. This is true not only for the large-scale and small-scale lobes in Fig. 3, but also for the H₂O masers distribution, as shown in Fig. 9. The existence of the same asymmetry on different scales and tracers cannot be a coincidence and is suggestive of a mechanism that causes a delay of the ejection of the SW lobe with respect to the NE lobe. In our opinion, two explanations are possible: either the accretion (and hence ejection) event is intrinsically stronger on the NE side of the disk, or the expansion towards the SW is hindered by the presence of denser material. We favour the latter hypothesis since the near-IR emission is much fainter from the SW lobe than from the NE lobe (see Caratti o Garatti et al. 2017 and Fedriani et al. 2023). This finding is surprising, because the SW lobe corresponds to the blue-shifted emission of the jet outflow (see Wang et al. 2011), which means that the jet is pointing towards the observer on that side and the extinction should be lower than on the NE side. The weakness of the IR emission to the SW is thus indicative of an asymmetry in the density distribution along the jet axis, a fact that could naturally also explain the delay of the expansion of the SW lobe with respect to the NE lobe.

Fig. 9 Water masers spots (solid circles) observed by Goddi et al. (2007; cyan), Burns et al. (2016; green), and Hirota et al. (2021; magenta) overlaid on the 3 mm continuum map (blue contours) of Fig. 1b. The dashed lines indicate the jet axes in the two tracers. The cross marks the position of the (proto)star. The circle in the bottom left represents the synthesised beam at 3 mm.

4.3 Evolution of the radio emission

In this section, we present a model fit to the observed spectra, following the approach adopted in Paper I, with some modifications. As done in Paper I, we adopted the ‘standard spherical’ model from (Reynolds 1986, see his Table 1) to describe the radio continuum emission from the jet. In practice, this means that we assume a jet where at a given time the opening angle, electron temperature, ionisation degree, and expansion speed do not depend on the distance from the star, r.

4.3.1 Expansion law of the jet

In Paper I, the jet was assumed to undergo expansion at constant velocity so that the maximum radius could be described by the simple expression r_m(t) = r_m(0) + v₀t, with v₀ = 900 km s⁻¹. While this assumption could hold for the first few months after the onset of the radio outburst, it is inconsistent with the most recent data. In fact, in Sect. 3.2 we show that the jet expansion is slowing down with time. We thus need to adopt a more realistic law for r_m(t).

For this purpose, we assume that the jet is expanding in a medium with density ∝ r⁻², with r the distance from the star. It is possible to demonstrate (see Appendix A) that applying momentum conservation one obtains $$y_{\text{m}}(t)=y_{\text{m}0}+2T{v}_0\cos \psi \left(\sqrt{1+\frac{t}{T}}-1\right),$$(4)

where we have multiplied both terms of Eq. (A.3) by cos ψ, with ψ the angle between the jet axis and the plane of the sky. For consistency with Reynolds (1986), we have indicated with y the projection of r on the plane of the sky. Here, T is a suitable timescale, and v₀ and y_m0 are, respectively, the expansion velocity and the value of y_m at the onset of the radio outburst (i.e. on July 10, 2016). As in Paper I, we chose v₀ = 900 km s⁻¹, while the two parameters T and y_m0 can be determined from the values of y_m estimated from the two ALMA maps at t = 1061 days and t =1881 days.

The problem is that y_m cannot be trivially obtained from the position of the NE peak, which corresponds to the maximum brightness temperature. This temperature is attained either at the inner radius, if the whole jet is optically thin, or at the border between the optically thin and the optically thick parts of the jet (denoted by y₁ in Reynolds’ notation). Beyond this point, the emission is optically thin and the brightness temperature scales with the opacity τ ∝ r⁻³, as from Reynolds’ Eq. (4). For this reason, the jet can extend much beyond the position of the observed peak.

In order to obtain a reliable estimate of y_m, we have fitted the NE lobe in the two ALMA maps of Fig. 1 assuming that the 3 mm emission is optically thin all over the jet surface. Under this hypothesis the brightness temperature can be expressed as $$T_{\text{B}}=T_0\tau =T_0{\left(\frac{y}{y_0}\right)}^{-3},$$(5)

where y = r cos ψ is the projection of r on the plane of the sky and y₀ = r₀ cos ψ with r₀ the inner radius of the jet. For a given set of y₀, y_m, and θ₀ (the opening angle of the jet) a map was generated, convolved with the instrumental beam, and the model brightness temperature was computed for each pixel of the observed map. The best fit was obtained by minimising the expression $${\chi }^2={\displaystyle \sum _i{\left(T_{\text{B}}^i(\text{model})-T_{\text{B}}^i(\text{data })\right)}^2},$$(6)

with i a generic pixel, after varying the three parameters over suitable ranges. The best-fit models are compared to the observed maps in Fig. 10 and correspond to θ₀ = 9°, y₀=0″.206=367 au, and y_m=0″.29=516 au, for the first map, and θ₀=2°4, y₀=0″.207=368 au, and y_m=0″.356=634 au, for the second map.

From Eq. (4) written for t=1061 days and t=1881 days, one obtains a system of two equations in the two unknowns T and y_m0, the solutions to which are T = 124 days and y_m0 = 251 au. Hence, one has $$y_{\text{m}}(\text{au})=251+127\left(\sqrt{1+\frac{t(\text{ days })}{124}}-1\right).$$(7)

The value of r_m can be computed from Eq. (4) assuming ψ =10° (see Paper I, where the complementary angle i = 90° − ψ = 80° was used). Figure 11 compares our solution, y_m, as a function of time (blue curve) with (i) the positions of the maximum brightness temperature (the same as in Fig. 5) and (ii) the value of y₁ obtained from the model fits described later in Sect. 4.3.2. The latter corresponds to the border between the optically thin and optically thick parts of the jet. Clearly, the brightness appears to peak much closer to y₁ than to y_m, as previously mentioned.

Fig. 10 Comparison between the observed maps of the NE lobe of the radio jet (left panels) with the best-fit maps obtained with the model (right panels). The top panels correspond to the first ALMA epoch, bottom panels to the second one. The horizontal axis is parallel to the jet axis. Contour levels range from 10% to 90% in steps of 10% of the peak of each image. The circle in the bottom left represents the synthesised beam.

4.3.2 Modelling the jet variability

We wish to reproduce the spectral variation shown in Fig. 4 and thus derive the values of the jet parameters as a function of time. For this purpose, we introduce some modifications with respect to the original model by Reynolds (1986), which was used in Paper I. Reynolds’ equations were derived under the approximation of small θ₀, the opening angle of the jet⁵. However, in Paper I we find that θ₀ ≃ 20°−50° is needed to fit the spectra. In order to overcome this limitation we re-wrote Reynolds’ equations under suitable assumptions, as detailed in Appendix B, so that they are now valid for any θ₀ < 90°.

All observed spectra have been fitted with Eq. (B.9) using only the measurements with v ≥ 6 GHz, for the reasons discussed in Sect. 3.2. We stress that, unlike Paper I, here we assume the jet to be mono-polar, because there is no hint of the existence of a SW lobe during the whole period of our VLA monitoring (see Sect. 4.1). The input parameters of the model are the angle between the jet axis and the plane of the sky, ψ, the ionised gas temperature, T₀, the inner radius, r₀, the projection of the outer radius on the plane of the sky, y_m, the opening angle, θ₀, and the parameter $\text{Λ}=x_0\dot{M}/v_0$, where x₀ is the fraction of ionised gas, v₀ the expansion velocity, and $\dot{M}$ the total mass loss rate (neutral plus ionised). The quantities T₀, v₀, and x₀ are assumed to be constant along the jet, while T₀ is also constant in time.

To simplify the fitting procedure as much as possible, we fixed T₀ = 10⁴ K and ψ = 10° (see Paper I), and computed y_m from Eq. (7). Unlike in Paper I, we decided to leave r₀ free, because a priori the inner radius could change while the jet is expanding. So, we are left with three free parameters: r₀, θ₀, and Λ. The best fit to each spectrum has been obtained by minimising the χ² given in Eq. (10) of Paper I, after varying the parameters over the ranges θ₀ = 3°−50°, r₀ = 50−300 au, and Λ = 10⁻¹⁰−10⁻⁷ M_⊙ yr⁻¹ (km s⁻¹)⁻¹. The best-fit spectra are represented by the solid curves in Fig. 4 and the best-fit parameters are given in Table 6. The errors on the parameters have been computed using the criterion of Lampton et al. (1976), as done in Paper I.

While most of the fits look to be in agreement with the data within the uncertainties, the 6 GHz fluxes appear to be underestimated by the model at the last five to six epochs. In our opinion, such a discrepancy could indicate contamination from non-thermal emission at this frequency, as already suggested by OLHKP. Indeed, as discussed in Sect. 3.2, non-thermal emission is very prominent at longer wavelengths in the same spectra (red points in Fig. 4) and it is hence not surprising that this type of emission can contribute significantly to the flux up to 6 GHz.

In Fig. 12, we plot the best-fit parameters as a function of time. One sees that the opening angle rapidly increases up to the maximum value allowed by us. It may seem that an opening angle of 50^° is too large for a jet, but it is similar to that predicted by theory for a jet powered by a ~20 M_⊙ YSO Zinchenko et al. (2015), namely ~52º, obtained by interpolating the values in Table 2 of Staff et al. (2019). Moreover, we can obtain a direct estimate of θ₀ from the observed peak brightness temperature in the synthesised beam, T_SB, assuming optically thick emission. Approximating the jet as a Gaussian source, one has $$T_{\text{SB}}=T_{\text{B}}\frac{{\text{Θ}}_{\text{S}}^2}{{\text{Θ}}_{\text{S}}^2+{\text{Θ}}_{\text{B}}^2},$$(8)

where T_B is the intrinsic brightness temperature of the source and Θ_S and Θ_B are, respectively, the full widths at half power of the source and synthesised beam. Consequently, $${\text{Θ}}_{\text{S}}={\text{Θ}}_{\text{B}}\sqrt{\frac{T_{\text{SB}}}{T_{\text{B}}-T_{\text{SB}}}}.$$(9)

This expression can be used to calculate a lower limit on the source diameter, ${\text{Θ}}_{\text{S}}^{\min }$, assuming that the free-free emission is optically thick, namely for T_B = T₀ = 10⁴ K. Correspondingly, a lower limit on the opening angle is obtained from $${\theta }_0^{\min }=\arcsin \left(\frac{{\text{Θ}}_{\text{S}}^{\min }}{2\text{Δ}}\right)=\arcsin \left(\frac{{\text{Θ}}_{\text{B}}}{2\text{Δ}}\sqrt{\frac{T_{\text{SB}}}{T_0-T_{\text{SB}}}}\right),$$(10)

where ∆ is the separation between the peak of T_SB and the position of the (proto)star. We estimated ${\theta }_0^{\min }$ for all of our maps, and for each epoch we computed the mean ${\theta }_0^{\min }$ obtained from the four bands. This is plotted in Fig. 13 as a function of time. For the sake of comparison we also plot θ₀ obtained from our model fits (dashed line). We conclude that, despite the crude approximations adopted to derive Eq. (10), values of θ₀ of a few times 10º seem plausible and consistent with our model fit results. It is worth noting that θ₀ computed from the map of OLHKP and from our ALMA maps (Fig. 1) using Eq. (10) turns out to be much less, namely 7º, 1º.7, and 1º.3, respectively 664, 1061, and 1881 days after the radio outburst. We speculate that the jet could be re-collimating on a timescale of a couple of years.

Another interesting result from Fig. 12 is the behaviour of r₀, which remains basically constant for ~250 days after the radio outburst, thereafter showing a systematic increase until the end of our monitoring – and even beyond that, as we estimate values of ~370 au at the time of our first ALMA observations (see Sect. 4.3.1). A straightforward interpretation is that the inner radius expands when the mechanism feeding the jet is switched off. Noticeably, Λ appears to reach a maximum just when r₀ starts increasing, supporting the idea that no more material is added to the jet. This hypothesis is confirmed by Fig. 14, which shows how the ionised jet mass, computed from Eq. (B.13), varies with time. In the same figure we also plot r₀ for the sake of comparison. It is quite clear that both quantities have a bi-modal temporal behaviour, before and after the time interval marked by the grey area. From the onset of the radio outburst up to ~250 days, r₀ remains approximately constant, whereas the jet mass increases. After that, the reverse occurs: as soon as the inner radius starts increasing, the jet stops growing in mass, with the grey area marking the transition between these two phases. This is consistent with mass conservation in an expanding jet that is not fed anymore by the outburst.

It is also worth noting that the ratio, R_e, between the mass ejected, M_jet, and that accreted during the outburst, M_acc, can be computed from the final value of x₀ M_jet in Fig. 14 and that quoted by (Caratti o Garatti et al. 2017; M_acc ≃ 3.4 × 10⁻³ M_⊙). One obtains R_e ≃ 7.5 × 10⁻⁵/(3.4 × 10⁻³ x₀) ≃ 2.2 × 10⁻²/x₀. Since by definition x₀ ≤ 1, we conclude that R_e must be greater than a few percent. Vice versa, assuming that at least 10% of the infalling material will be redirected into outflow, one has R_e > 0.1, which sets an upper limit of ~0.2 on the ionisation fraction, consistent with the estimate obtained in Paper I and the values estimated for similar sources (see Fedriani et al. 2019).

Fig. 11 Comparison of the projection on the plane of the sky of the outer radius of the jet, ym (blue), with the position of the brightness temperature peak, ypeak (black), and the border between the optically thin and optically thick parts of the jet, y1 (red). All quantities are plotted as a function of time from the onset of the radio outburst (July 10, 2016). The black circles indicate our VLA data, the triangle the OLHKP data, and the two squares our ALMA data. We remark that ypeak is measured from Fig. 5 as the separation between the star and the projection of the peak along the jet axis.

Fig. 12 Parameters of the best fits to the continuum spectra, as a function of time from the onset of the radio outburst (July 10, 2016). The plotted parameter is indicated in the bottom right of each panel.

Fig. 13 Lower limit on the opening angle estimated from Eq. (10) for all the epochs of our monitoring. The bars indicate the standard deviation of the values obtained in the different bands. The dashed line connects the values of θ0 obtained from the best fits to the spectra (see Table 6).

Fig. 14 Ionised jet mass (solid black curve) and inner radius (dashed red curve) as a function of time from the onset of the radio outburst (July 10, 2016). The grey area separates the period when the jet is still powered by the outburst from that when the feeding mechanism is quenched.

5 Summary and conclusions

As a follow-up to Paper I, we monitored the radio continuum emission of the outburst from the massive YSO S255IR NIRS 3 over ~13 months at six wavelengths with the VLA. We also imaged the radio jet at 3 mm with ALMA at two epochs separated by ~27 months, with an angular resolution ≲0″.1. Our results indicate that after an exponential increase in the radio flux in all observed bands, the intensity becomes constant or slightly decreasing. A comparison of the two ALMA maps shows that the radio jet is expanding both to the NE and SW, although only the NE lobe was present during our VLA monitoring. The SW lobe appeared between May 2018 and June 2019, namely at a much later time than the NE lobe. We believe that this ejection event is related to the same accretion outburst, which occurred in 2015. We speculate that the delay between the two lobes might be due to a greater density of the medium facing the SW lobe, which could curb the expansion of the jet on that side.

From the analysis of the continuum spectra, we infer that two mechanisms are needed to explain the observed fluxes: free-free emission at short wavelengths and non-thermal (probably synchrotron) emission at long wavelengths. We believe that the latter should become dominant at frequencies ≲6 GHz. For this reason, we fitted only the data with v ≥ 6 GHz using a slightly modified version of the jet model adopted in Paper I, which works for any opening angle of the jet <90º. We conclude that the spectra can be satisfactorily reproduced with an expanding, decelerating, mono-polar thermal jet that was actively powered by the outburst until mid-2017. After this date, no more mass is injected into the lobes and the inner jet radius expands, which, over the long term, is bound to give rise to one of the knots that characterise thermal jets from YSOs.

Acknowledgements

A.C.G. acknowledges from PRIN-MUR 2022 20228JPA3A “The path to star and planet formation in the JWST era (PATH)” and by INAF-GoG 2022 “NIR-dark Accretion Outbursts in Massive Young stellar objects (NAOMY)” and Large Grant INAF 2022 “YSOs Outflows, Disks and Accretion: towards a global framework for the evolution of planet forming systems (YODA)”. R.F. acknowledges support from the grants Juan de la Cierva JC2021-046802-I, PID2020-114461GB-I00 and CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033 and by “European Union NextGenera-tionEU PRTR”. T.P.R. acknowledges support from ERC grant 743029 EASY. This study is based on observations made under project 16A-424, 16B-427, and 17B-045 of the VLA of NRAO. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.. This paper makes also use of the following ALMA data: ADS/JAO.ALMA#2018.1.00864.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

Appendix A Jet expansion law

As discussed in Sect. 4.3.1, the jet is not expanding at constant velocity during the period of our monitoring, but it appears to slow down. It is thus necessary to adopt an expression for the maximum radius, r_m, that properly takes this effect into account. A reasonable scenario may be that of a jet confined in a solid angle Ω, with initial mass M₀, which expands with initial velocity υ₀ through a medium with density ρ ∝ r⁻², where r is the distance from the jet origin. Because of momentum conservation one can write $$\begin{array}{l}{M}_0{v}_0=M(t)\frac{\text{d}{r}_{\text{m}}}{\text{d}t}\\ =\left[M_0+{\displaystyle {\int }_{r_{\text{m}0}}^{r_{\text{m}}}{\rho }_0}{\left(\frac{r_{\text{m}0}}{r}\right)}^2\text{Ω}{r}^2 \text{d}r\right]\frac{\text{d}{r}_{\text{m}}}{\text{d}t}\\ =\left[M_0+\text{Ω}{\rho }_0{r}_{\text{m}0}^2\left(r_{\text{m}}-r_{\text{m}0}\right)\right]\frac{\text{d}{r}_{\text{m}}}{\text{d}t},\end{array}$$(A.1)

with ρ₀ density at radius r_m0. The solution of this differential equation is $$\begin{array}{l}{M}_0{v}_{0}t=M_0\left(r_{\text{m}}-r_{\text{m}0}\right)\\ +\text{Ω}{\rho }_0{r}_{\text{m}0}^2\left[\frac{r_{\text{m}}^2-r_{\text{m}0}^2}{2}-r_{\text{m}0}\left(r_{\text{m}}-r_{\text{m}0}\right)\right],\end{array}$$(A.2)

which after some algebra gives $$r_{\text{m}}(t)=r_{\text{m}0}+2T{v}_0\left(\sqrt{1+\frac{t}{T}}-1\right),$$(A.3)

where we have defined $T\equiv M_0/\left(2{\rho }_0\text{Ω}{r}_{\text{m}0}^2{v}_0\right)$.

While this expression provides us with a more realistic description of the jet expansion than the constant velocity assumption adopted in Paper I, we stress that it is not to be taken as the real equation of motion of the jet but as the simplest way to parametrise the observed deceleration of it.

Appendix B Description of the model

In Paper I we adopted the jet model by Reynolds (1986) to describe the integrated flux density of the radio jet from S255IR NIRS 3. More specifically, we assumed what is defined as the ‘standard spherical’ case in Reynolds’ Table 1, namely a jet where the opening angle (θ₀), internal velocity (υ₀), ionisation degree (x₀), and temperature (T₀) do not depend on the distance r from the star. Conservation of mass along the flow implies that the gas number density can be expressed as n = n₀(r/r₀)⁻².

Despite its simplicity, the model was successful in fitting the observed spectra in Paper I. However, Reynolds’ equations have been derived under the assumptions of small θ₀, an approximation that is not satisfied by the best fit to the radio spectra obtained in Paper I (see Table 3 there), which requires angles as large as ~50º. In order to overcome this limitation, we propose here a slightly modified version of Reynolds’ model that works for any θ₀ < 9oº.

To allow for an analytic solution of the equations, we maintain Reynolds’ assumption that the opacity depends only on r. This is equivalent to assuming that the jet is not conical but has a pyramidal shape with two faces parallel to the line of sight. We also assumed that the jet is delimited by two cylindrical surfaces, with radii r₀ and r_m, and that its axis is inclined by a small angle, ψ, with respect to the plane of the sky. Figure B.1 schematically illustrates the geometry of the jet, where the star lies at the origin of the axes and the line of sight is parallel to the ɀ-axis.

Fig. B.1 Sketch of the jet model. The x- and y-axes lie in the plane of the sky, and the ɀ-axis is parallel to the line of sight. The y-axis coincides with the projection of the jet axis on the plane of the sky. The dotted line labelled ‘l.o.s.’ denotes a generic line of sight. The star powering the jet lies at the origin of the coordinate system.

Following Reynolds, the absorption coefficient of the ionised gas and the opacity can be written as $$\kappa (r)={\kappa }_0{\left(\frac{r}{r_0}\right)}^{-4},$$(B.1)

where $${\kappa }_0=a_{\kappa }{n}_0^2{x}_0^2{T}_0^{-1.35}{v}^{-2.1}\text{, }$$(B.2)

with a_ĸ = 0.212 in CGS units, and $$\begin{array}{l}\tau (R)={\displaystyle {\int }_{-R\tan \left({\theta }_0-\psi \right)}^{R\tan \left({\theta }_0+\psi \right)}\kappa }\text{d}z={\kappa }_0{r}_0^4{\displaystyle {\int }_{-R\tan \left({\theta }_0-\psi \right)}^{R\tan \left({\theta }_0+\psi \right)}{\left(R^2+z^2\right)}^{-2}} \text{d}z\\ =\frac{{\kappa }_0{r}_0^4}{2R^3}\left[\frac{\tan \left({\theta }_0+\psi \right)}{1+{\tan }^2\left({\theta }_0+\psi \right)}+\frac{\tan \left({\theta }_0-\psi \right)}{1+{\tan }^2\left({\theta }_0-\psi \right)}+2{\theta }_0\right]\\ =\frac{{\kappa }_0{r}_0^4}{2R^3}\left[\sin \left(2{\theta }_0\right)\cos (2\psi )+2{\theta }_0\right]\\ ={\tau }_0{\left(\frac{r_0}{R}\right)}^3.\end{array}$$(B.3)

Here we have defined the two quantities $$R=\sqrt{x^2+y^2}$$(B.4)

and $${\tau }_0={\kappa }_0{r}_0\frac{\sin \left(2{\theta }_0\right)\cos (2\psi )+2{\theta }_0}{2}$$(B.5)

To ease the comparison with Reynolds’ expressions, we indicate with y₀ and y_m the projections of r₀ and r_m on the plane of the sky, namely y₀ = r₀ cos ψ and y_m = r_m cos ψ, and with y₁ the value of R at which τ(R) = 1. From Eq. (B.3) one has $$y_1=r_0{\tau }_0^{\frac{1}{3}}.$$(B.6)

For the calculation of the total flux density of the jet we follow Reynolds’ approach and assume that the emission between R = 0 and R = y₁ is optically thick, and that between R = y₁ and R = y_m is optically thin. Under this approximation the flux can be written as $$\begin{array}{l}{S}_v=\frac{1}{d^2}{\displaystyle {\int }_{y_0}^{y_{\text{m}}}{B}_v}\left(T_0\right)\left(1-{\text{e}}^{-\tau }\right)2{{\theta }^{\prime }}_{0}R \text{d}R\\ \simeq \frac{2{{\theta }^{\prime }}_0{B}_v\left(T_0\right)}{d^2}\left[{\displaystyle {\int }_{y_0}^{y_1}R} \text{d}R+{\displaystyle {\int }_{y_1}^{y_{\text{m}}}\tau }(R)R \text{d}R\right]\\ =\frac{2{{\theta }^{\prime }}_0{B}_v\left(T_0\right)}{d^2}\left[\frac{y_1^2-y_0^2}{2}+y_1^3\left(\frac{1}{y_1}-\frac{1}{y_{\text{m}}}\right)\right],\end{array}$$(B.7)

where B_v is the Planck function and θ′₀ is the projection of θ₀ on the plane of the sky. The two angles are related by the expression $$\tan {{\theta }^{\prime }}_0=\frac{\tan {\theta }_0}{\cos \psi }.$$(B.8)

We note that it is not strictly correct to integrate from y₀ to y_m because the projections of the two circles with radii r₀ and r_m on the plane of the sky are ellipses, and the integral should be made not only in the variable R but also in the azimuthal angle. However, the approximation adopted by us is acceptable for small ψ.

The expression of S_v in the case of a jet totally thin (y₁ < y₀) or thick (y₁ > y_m) can be obtained in a similar way, by considering only the relevant approximation in the argument of the integral in Eq. (B.7). In conclusion, one obtains $$\begin{array}{l}{S}_v=\frac{2{{\theta }^{\prime }}_0{B}_v\left(T_0\right)}{d^2}\\ \times \{\begin{array}{ll}{y}_1^3\left(\frac{1}{y_0}-\frac{1}{y_{\text{m}}}\right)\hfill & \iff y_1\le y_0\hfill \\ \left[\frac{y_1^2-y_0^2}{2}+y_1^3\left(\frac{1}{y_1}-\frac{1}{y_{\text{m}}}\right)\right]\hfill & \iff y_0<y_1<y_{\text{m}}.\hfill \\ \frac{y_{\text{m}}^2-y_0^2}{2}\hfill & \iff y_1\ge y_{\text{m}}\hfill \end{array}\end{array}$$(B.9)

We note that this expression gives the total flux emitted by a single jet lobe, whereas Eq. (5) of Paper I takes into account both lobes. The different approach is justified by the fact that our new findings have proved that only the NE lobe was present during our monitoring (see Sect. 4.1).

The quantity y₁ is a function of n₀, in addition to other parameters. However, following Reynolds, it may be convenient to express it in terms of the mass loss rate of the jet, $\dot{M}$, which is obtained by integrating the flux of mass through the inner surface of the jet. Obviously, $\dot{M}$ does not depend on the inclination of the jet with respect to the line of sight, so we can simplify the calculation by assuming ψ = 0. Hence, we obtain $$\begin{array}{l}\dot{M}={\displaystyle {\int }_{-r_0\tan {\theta }_0}^{r_0\tan {\theta }_0}\mu }{n}_0\frac{r_0^2}{r_0^2+z^2}{\upsilon }_0\frac{r_0}{\sqrt{r_0^2+z^2}}2{\theta }_0{r}_0 \text{d}z\\ =4\mu n_0{\upsilon }_0{\theta }_0{r}_0^4{\displaystyle {\int }_0^{r_0\tan {\theta }_0}{\left(r_0^2+z^2\right)}^{-\frac{3}{2}}} \text{d}z\\ =4\mu n_0{\upsilon }_0{\theta }_0{r}_0^2\frac{\tan {\theta }_0}{\sqrt{1+{\tan }^2{\theta }_0}},\end{array}$$(B.10)

where µ is the mean particle mass per hydrogen atom (we assume µ = 1.67 × 10⁻²⁴ g) and ${\upsilon }_0\frac{r_0}{\sqrt{r_0^2+z^2}}$ is the component of the velocity perpendicular to the inner surface of the jet. From this it is trivial to express n₀ as a function of $\dot{M}$, and from Eqs. (B.2), (B.5), and (B.6) one obtains $$y_1={\left(\frac{a_{\kappa }}{16{\mu }^2}{\Lambda }^2{T}_0^{-1.35}{v}^{-2.1}\frac{1+{\tan }^2{\theta }_0}{{\tan }^2{\theta }_0}\frac{\sin \left(2{\theta }_0\right)\cos (2\psi )+2{\theta }_0}{2{\theta }_0^2}\right)}^{\frac{1}{3}},$$(B.11)

where we have defined $\text{Λ}\equiv x_0\dot{M}/{\upsilon }_0$.

Finally, a quantity of interest for our purposes is the ionised mass of the jet, which is given by the expression $$\begin{array}{l}{M}_{\text{i}}={\displaystyle {\int }_{-{\theta }_0}^{{\theta }_0} \text{d}}\theta {\displaystyle {\int }_{r_0}^{r_{\text{m}}}R} \text{d}R{\displaystyle {\int }_{-R\tan {\theta }_0}^{R\tan {\theta }_0}\mu }{x}_0{n}_0\frac{r_0^2}{R^2+z^2} \text{d}z\\ =4\mu x_0{n}_0{\theta }_0^2{r}_0^2\left(r_{\text{m}}-r_0\right)\end{array}$$(B.12) $$=\text{Λ}{\theta }_0\frac{\sqrt{1+{\tan }^2{\theta }_0}}{\tan {\theta }_0}\left(r_{\text{m}}-r_0\right),$$(B.13)

where we have used Eq. (B.10) to replace n₀.

References

All Tables

All Figures

Fig. 1 Maps of the 3 mm continuum emission obtained with ALMA. (a) Data acquired on June 6, 2019. Contour levels range from 0.52 to 16.12 in steps of 2.6 mJy beam−1. The cross marks the peak of the 900 μm continuum emission imaged by Liu et al. (2020). The circle in the bottom left represents the synthesised beam. (b) Same as the top panel, but for the map obtained on September 3, 2021. Contour levels range from 0.29 to 16.24 in steps of 1.45 mJy beam−1.

In the text

Fig. 2 Maps of the line emission obtained with ALMA. (a) Map of the emission averaged over the H13CN(1−0) line (white contours) overlaid on the continuum image at 3 mm obtained on June 6, 2019. Contour levels range from 1.2 to 4.56 in steps of 0.48 mJy beam−1. The circle in the bottom left represents the synthesised beam. (b) Same as the top panel, but for the SO(22−11) line. Contour levels range from 1.65 to 4.62 in steps of 0.33 mJy beam−1.

In the text

	Fig. 3 Maps of the 3.6 cm continuum emission obtained with the VLA on July 10, 2016 (red) and of the 3 mm continuum emission (blue), also shown in Fig. 1b. The dashed lines denote the symmetry axes of the bipolar structures. Red contour levels range from 0.1 to 0.6 in steps of 0.1 mJy beam−1. Blue contour levels range from 0.29 to 8.99 in steps of 2.9 mJy beam−1.
In the text

Fig. 4 Continuum spectra of the variable, compact component in S255IR NIRS 3. The date of the observations is given in the top left of each panel. The vertical error bars assume a calibration error of 10% in all bands, while the horizontal bars indicate the bandwidth covered at each frequency. The upper limits denote that the flux density is also contributed by the emission from the large-scale lobes, due to the limited angular resolution. Red symbols mark data in the L (1.5 GHz) and S (3 GHz) bands. The curves are the best fits obtained with the model of a thermal jet described in Sect. 4.3.2.

In the text

Fig. 5 Distribution of the peaks of the 3 mm and 1.3 cm continuum maps obtained, respectively, with ALMA and with the A and B configurations of the VLA. The circles indicate our VLA data, the triangle the data of OLHKP, and the two squares our ALMA data. The dashed line connects the points in order of time (the first point in the bottom-right corner corresponds to the observation on July 10, 2016). The third point from the top left is obtained from the OLHKP data. The cross marks the position of the (proto)star, and the solid red line is the axis of the radio jet (PA = 48°). The dotted ellipses are the synthesised beams scaled down by a factor of 5.

In the text

Fig. 6 Flux density of the variable, compact source in S255IR NIRS 3 as a function of time from the beginning of the radio burst (July 10, 2016) in all bands observed with the VLA. The two rightmost points are the measurements at 6 and 22 GHz by OLHKP. Curves are colour-coded by frequency (in gigahertz). Arrows denote points with upper limits (see Sect. 2.1), which are connected by dashed lines. The two black arrows pointing to the right indicate the flux densities at 92 GHz of the NE lobe (see Table 1) measured 1061 and 1881 days after the radio burst. Typical error bars are shown in the bottom right of the figure. The vertical dotted lines mark the beginning and end of an array configuration, as labelled above the figure.

In the text

Fig. 7 Maps of the 3 mm continuum emission (contours) overlaid on the corresponding map of the spectral index at 3 mm. (a) Data obtained with ALMA on June 6, 2019. Black contour levels range from 10% to 90% in steps of 10% of the peak emission. White contour levels correspond to a spectral index equal to 0. The cross marks the position of the (proto)star. The circle in the bottom left represents the synthesised beam. (b) Same as the top panel, but for the map obtained on September 3, 2021.

In the text

	Fig. 8 Map of the 3 mm continuum emission obtained with ALMA on September 3, 2021, overlaid on the 1.3 cm image obtained by OLHKP on May 5, 2018. Contour levels are the same as in Fig. 1b. The cross marks the position of the (proto)star. The ellipse in the bottom left represents the synthesised beam at 1.3 cm.
In the text

	Fig. 9 Water masers spots (solid circles) observed by Goddi et al. (2007; cyan), Burns et al. (2016; green), and Hirota et al. (2021; magenta) overlaid on the 3 mm continuum map (blue contours) of Fig. 1b. The dashed lines indicate the jet axes in the two tracers. The cross marks the position of the (proto)star. The circle in the bottom left represents the synthesised beam at 3 mm.
In the text

Fig. 10 Comparison between the observed maps of the NE lobe of the radio jet (left panels) with the best-fit maps obtained with the model (right panels). The top panels correspond to the first ALMA epoch, bottom panels to the second one. The horizontal axis is parallel to the jet axis. Contour levels range from 10% to 90% in steps of 10% of the peak of each image. The circle in the bottom left represents the synthesised beam.

In the text

Fig. 11 Comparison of the projection on the plane of the sky of the outer radius of the jet, ym (blue), with the position of the brightness temperature peak, ypeak (black), and the border between the optically thin and optically thick parts of the jet, y1 (red). All quantities are plotted as a function of time from the onset of the radio outburst (July 10, 2016). The black circles indicate our VLA data, the triangle the OLHKP data, and the two squares our ALMA data. We remark that ypeak is measured from Fig. 5 as the separation between the star and the projection of the peak along the jet axis.

In the text

	Fig. 12 Parameters of the best fits to the continuum spectra, as a function of time from the onset of the radio outburst (July 10, 2016). The plotted parameter is indicated in the bottom right of each panel.
In the text

	Fig. 13 Lower limit on the opening angle estimated from Eq. (10) for all the epochs of our monitoring. The bars indicate the standard deviation of the values obtained in the different bands. The dashed line connects the values of θ0 obtained from the best fits to the spectra (see Table 6).
In the text

	Fig. 14 Ionised jet mass (solid black curve) and inner radius (dashed red curve) as a function of time from the onset of the radio outburst (July 10, 2016). The grey area separates the period when the jet is still powered by the outburst from that when the feeding mechanism is quenched.
In the text

	Fig. B.1 Sketch of the jet model. The x- and y-axes lie in the plane of the sky, and the ɀ-axis is parallel to the line of sight. The y-axis coincides with the projection of the jet axis on the plane of the sky. The dotted line labelled ‘l.o.s.’ denotes a generic line of sight. The star powering the jet lies at the origin of the coordinate system.
In the text