Radio outburst from a massive (proto)star II. A portrait in space and time of the expanding radio jet from S255IR NIRS 3 ⋆

Context. Growing observational evidence indicates that the accretion process leading to star formation may occur in an episodic way, through accretion outbursts revealed in various tracers. This phenomenon has also now been detected in association with a few young massive (proto)stars ( > 8 M ⊙ ), where an increase in the emission has been observed from


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/submillimeter Array (ALMA; see e.g. the review by Beltrán & de Wit 2016).This strongly suggests that diskmediated 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(Stecklum , 2021;;Caratti o Garatti et al. 2017;Hunter et al. 2017Hunter et al. , 2021;;Burns et al. 2020Burns et al. , 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.12  −0.11 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.

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.

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 phasecalibrators were J0632+1022 in the L and S bands, J0559+2353 in the C and X bands, and J0539+1433 in the K and Q 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 CASA1 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 ∼15000 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 m 53 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 −1 , 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 software2 developed by Sánchez-Monge et al. (2018).In this way we also obtained cubes of the continuumsubtracted 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 2 -1 1 ) at 86093.983MHz and H 13 CN(1-0) at 86338.7 MHz.Note: The row labelled with 'days' gives the number of days after the beginning of the radio burst.The error on the flux densities is estimated to be ∼10% (see main text).

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 Fig. 1a and 1b reveals that, while the NE peak is mov-ing 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 2 -1 1 ) and H 13 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 contin-  71 37, 2.3 342, 7.7  K -22.2  16, 0.094 19, 0.36 36, 1.1 111, 3.4  Q -45.5  81, 0.061 42, 0.20 51, 0.57 102, 1.9  this paper H 2 O masers 0 ′′ .5 (900 au) 44 • Hirota et al. (2021)   uum.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 submillimetre 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.

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 −1 , much less than the mean velocity from the beginning of the radio burst, namely ∼472 au/1881 days≃436 km s −1 , 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 −1 from their IR data.We further analyse the expansion of the jet in Sect.4.3.1.four fluxes in a log S ν -log ν 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 ∝ ν γ , with γ = 2-4, because the dust absorption coefficient is believed to vary as ν β 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 ∝ ν −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 ν = S d ν + S ff ν , where S ff ν ∝ ν −0.1 is the optically thin free-free flux and S d ν ∝ ν γ , with γ = 2-4, is the dust flux.As previously explained, the spectral index between ν 1 = 85.2GHz and ν 2 = 99.2GHz was estimated assuming S ν ∝ ν α , and hence we have  b Both dust and free-free emission.
After some algebra, we obtain the ratio 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 Fig. 1a to Fig. 1b.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.

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 observations3 of the Maser Monitoring Organization (M2O)4 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 2 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 andFedriani 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.

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.

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)+ 0 t, with 0 =900 km s −1 .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 −2 , with r the distance from the star.It is possible to demonstrate (see Appendix A) that applying momentum conservation one obtains 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 0 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 0 =900 km s −1 , 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 1 in Reynolds' notation).Beyond this point, the emission is optically thin and the brightness temperature scales with the opacity τ ∝ r −3 , 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 where y = r cos ψ is the projection of r on the plane of the sky and y 0 = r 0 cos ψ with r 0 the inner radius of the jet.For a given set of y 0 , y m , and θ 0 (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 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 θ 0 =9 • , y 0 =0 ′′ .206=367 au, and y m =0 ′′ .29=516 au, for the first map, and θ 0 =2 • .4,y 0 =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 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 1 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 1 than to y m , as previously mentioned.

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 θ 0 , the opening angle of the jet 5 .However, in Paper I we find that θ 0 ≃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 θ 0 < 90 • .All observed spectra have been fitted with Eq. (B.9) using only the measurements with ν≥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 Fig. 11.Comparison of the projection on the plane of the sky of the outer radius of the jet, y m (blue), with the position of the brightness temperature peak, y peak (black), and the border between the optically thin and optically thick parts of the jet, y 1 (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 y peak is measured from Fig. 5 as the separation between the star and the projection of the peak along the jet axis. 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 0 , the inner radius, r 0 , the projection of the outer radius on the plane of the sky, y m , the opening angle, θ 0 , and the parameter Λ = x 0 Ṁ/ 0 , where x 0 is the fraction of ionised gas, 0 the expansion velocity, and Ṁ the total mass loss rate (neutral plus ionised).The quantities T 0 , 0 , and x 0 are assumed to be constant along the jet, while T 0 is also constant in time.
To simplify the fitting procedure as much as possible, we fixed T 0 = 10 4 K and ψ = 10 • (see Paper I), and computed y m from Eq. ( 7).Unlike in Paper I, we decided to leave r 0 free, because a priori the inner radius could change while the jet is expanding.So, we are left with three free parameters: r 0 , θ 0 , and Λ.The best fit to each spectrum has been obtained by minimising the χ 2 given in Eq. ( 10) of Paper I, after varying the parameters over the ranges θ 0 =3 • -50 • , r 0 =50-300 au, and Λ=10 −10 -10 −7 M ⊙ yr −1 /(km s −1 ).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   6).estimate of θ 0 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 where T B is the intrinsic brightness temperature of the source and Θ S and Θ B are, respectively, the full widths at half power of a obtained from Eq. ( 7) the source and synthesised beam.Consequently, This expression can be used to calculate a lower limit on the source diameter, Θ min S , assuming that the free-free emission is optically thick, namely for T B =T 0 =10 4 K. Correspondingly, a lower limit on the opening angle is obtained from where ∆ is the separation between the peak of T SB and the position of the (proto)star.We estimated θ min 0 for all of our maps, and for each epoch we computed the mean θ min 0 obtained from the four bands.This is plotted in Fig. 13 as a function of time.For the sake of comparison we also plot θ 0 obtained from our model fits (dashed line).We conclude that, despite the crude approximations adopted to derive Eq. ( 10), values of θ 0 of a few times 10 • seem plausible and consistent with our model fit results.It is worth noting that θ 0 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 0 , 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 0 starts increasing, supporting the idea that no more material is added to Fig. 14.Plot of the 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.
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 0 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 0 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 0 M jet in Fig. 14  by definition x 0 ≤ 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).

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 ν ≥ 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.

Fig. 1 .
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.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.

Fig. 2 .
Fig. 2. Maps of the line emission obtained with ALMA.a. Map of the emission averaged over the H 13 CN(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.The circle in the bottom left represents the synthesised beam.b.Same as the top panel, but for the SO(2 2 -1 1 ) line.Contour levels range from 1.65 to 4.62 in steps of 0.33 mJy/beam.

Fig. 3 .
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.Blue contour levels range from 0.29 to 8.99 in steps of 2.9 mJy/beam.

Fig. 4 .
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 .
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 .
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 .
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 .
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 circle in the bottom left represents the synthesised beam at 1.3 cm.

Fig. 9 .
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.

Fig. 10 .
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.

Fig. 12 .
Fig. 12. Plot of the 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 .
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 Table6).

Table 1 .
Flux densities of S255IR NIRS 3 observed at different epochs and frequencies.The value in parentheses is the flux density of the NE lobe only (see Sect. 4.1). f

Table 2 .
Typical 1σ RMS noise and synthesised beam for different bands and VLA configurations.

Table 3 .
Main parameters of the ALMA observations.

Table 5 .
Flux densities in millijanskys measured at 92.5 GHz with ALMA towards the NE and SW components of S255IR NIRS 3.
a Free-free emission.

Table 6 .
Parameters of the best fits to the spectra in Fig.4.