© The Authors 2024

1. Introduction

Circumstellar discs and jets are believed to play a crucial role in the formation of stars of all masses. While many objects of this type are known for solar-type stars, disc–jet systems associated with young early-type stars have only been revealed and studied in detail for a limited number of cases (Beltrán & de Wit 2016). The advent of the Atacama Large Millimeter and submillimeter Array (ALMA) as well as the upgrade of the Karl Jansky Very Large Array (VLA) have contributed to improving our knowledge in this field. As a consequence of such a substantial technical improvement, it is not only possible to perform large surveys of many disc–jet candidates, but also to monitor their emission and reveal changes in both the flux density and morphology due to the expansion and evolution of the jet.

With all of the above in mind, we recently conducted an observational campaign in a variety of tracers (Stecklum et al. 2016, 2021; Caratti o Garatti et al. 2017; Moscadelli et al. 2017; Hirota et al. 2021) of the young stellar object (YSO) S255IR NIRS 3 (distance 1.78 kpc; Burns et al. 2016), which is one of the few cases of an accretion outburst detected towards a young massive star (Tapia et al. 2015; Hunter et al. 2017, 2021; Burns et al. 2020, 2023; Chen et al. 2021). As part of this campaign, we made multi-epoch observations of the source at centimetre (with the VLA) and millimetre (with the NOEMA and ALMA interferometers) wavelengths (Cesaroni et al. 2018, hereafter Paper I; Cesaroni et al. 2023, hereafter Paper II), showing that the emission from the ionised jet also underwent a burst, although with a delay of about 1 year with respect to the infrared burst detected by Caratti o Garatti et al. (2017). The flux-density variation is explained with a model that takes into account the expansion of a conical ionised jet. This model can reproduce both the continuum spectrum and the morphology of the jet.

One of the results obtained in our studies is that although the jet initially consists of only the NE lobe (see also Fedriani et al. 2023), a SW lobe appears between approximately 22 and 35 months after the onset of the radio outburst, in July 2016. We propose that this new lobe is powered by the same accretion outburst (which began in mid-June 2015; see Caratti o Garatti et al. 2017) as the NE lobe, despite the significant time lag between the two. In Paper II, we speculated that the delay is due to an inhomogeneous distribution of the material around the YSO, with the gas density to the SW being greater than that to the NE of the star. To shed light on these issues, we performed new VLA subarcsecond-resolution observations of S255IR NIRS 3 to investigate the outbreak of the SW lobe, as well as its intensity variation with time. In the following, we report the results obtained.

2. Observations

We used the A-array configuration of the VLA on July 21, 2023 (project code: 23A-021), to image the continuum emission of the radio jet from S255IR NIRS 3 at three frequencies: 6, 10, and 22.2 GHz, which correspond to the C, X, and K bands, respectively.

The Wideband Interferometric Digital Architecture (WIDAR) correlator was used in dual polarisation mode. The total observing bandwidth (per polarisation) was 4 GHz in the C and X bands and 8 GHz in the K band. The primary flux calibrator was 3C48 and the phase calibrators were J0559+2353 in the C and X bands, and J0539+1433 in the K band.

We made use of the calibrated data set provided by the NRAO pipeline and subsequent inspection of the data and imaging were performed with the CASA¹ package, version 5.6.2-2. For the continuum images, we adopted natural weighting in the X and K bands and ‘Briggs’ weighting with ‘robust=0.5’ in the C band, which provides a good compromise between angular resolution and sensitivity to extended structures. The half-power widths and position angles of the synthesised beams are 0$\stackrel{″}{.}$30 × 0$\stackrel{″}{.}$28, 0$\stackrel{°}{.}$3 at 6 GHz, 0$\stackrel{″}{.}$23 × 0$\stackrel{″}{.}$22, −1$\stackrel{°}{.}$3 at 10 GHz, and 0$\stackrel{″}{.}$10 × 0$\stackrel{″}{.}$093, 0° at 22.2 GHz. The typical noise is ∼0.02 mJy beam⁻¹ in all bands for an on-source time of 26 min in the C and X bands and 20 min in the K band. The calibration uncertainty is estimated to be 15% in all bands.

3. Results and analysis

Figure 1 shows the maps obtained at the three frequencies. It is quite evident that the jet structure is bipolar with a prominent lobe to the NE. The ratio between the intensities of the two lobes appears to change with frequency, which suggests different opacities, and therefore column densities, of the two lobes. To better investigate this aspect, we plot in Fig. 2 the maps of the spectral index computed from the ratio between the 6 GHz and 22 GHz maps after cleaning the latter with the same beam as the former. These maps, obtained from consecutive observations, are only compared in the higher intensity contours to minimise possible systematic uncertainties (see Sanna et al. 2018). One sees that the spectral index is steeper to the SW than to the NE, which is consistent with the free-free emission being thicker. This result agrees with the idea proposed in Paper II that the expansion of the SW lobe could be curbed by a density enhancement in that direction.

Fig. 1. Maps of the continuum emission from S255IR NIRS 3 obtained with the VLA. (a) Data acquired at 22.2 GHz on July 21, 2023. Contour levels are quantified by the colour scale to the right of the panel. The ellipse in the bottom right is the half-power width of the synthesised beam. The cross marks the peak of the 3 mm continuum emission from Paper II, which we assume to be the position of the star. (b) Same as top panel, but for the 10 GHz emission. c. Same as top panel, but for the 6 GHz emission.

Fig. 2. Map of the spectral index obtained from the ratio between the 6 GHz and 22.2 GHz images. The contours are a map of the 22.2 GHz continuum emission obtained with the same clean beam as the 6 GHz image. Contour levels range from 0.5 to 5.5 in steps of 1 mJy beam−1. The ellipse in the bottom right is the half-power width of the synthesised beam.

3.1. Expansion of the jet lobes

While in Paper II we established that the NE lobe is expanding with velocity decreasing with time, the data obtained in that study were insufficient to gain insight into the expansion of the SW lobe. With the new VLA image at 22.2 GHz, it is now possible to shed light on this aspect in particular. For this purpose, it is necessary to compare the ALMA maps acquired at 3 mm in 2019 and 2021 with the VLA 1.3 cm map obtained in 2023.

A caveat with this approach is that the 3 mm continuum is heavily contaminated by emission from dust grains, especially in the SW lobe, whereas the 1.3 cm continuum is mostly due to pure free-free emission. To circumvent this problem, we applied Eq. (3) of Paper II, which allows us to compute the expected fraction of the flux density contributed by free-free emission based on the measured spectral index at 3 mm. In practice, we corrected the flux in each pixel of the 3 mm continuum maps (see Fig. 1 of Paper II) using the corresponding spectral index (see Fig. 7 of Paper II). In the calculation, we adopted the dust spectral index of 2.08 derived in Eq. (4) of Paper I. We note that for some pixels, especially in the NE lobe, the spectral index lies below the minimum value of –0.1 expected for free-free or dust emission. In these cases, we used –0.1 for the calculation.

The ‘pure’ free-free maps obtained with this method are shown in Fig. 3, where a comparison is also made with our new 1.3 cm image. Although it is not possible to identify a precise peak of emission in the 2021 map, there is no doubt that the SW lobe becomes more elongated with time. A mean projected expansion speed of ∼285 km s⁻¹ is estimated from the separation of ∼0$\stackrel{″}{.}$14 = 248 au between the SW peak of the 2019 image and that of the 1.3 cm map, and the corresponding time interval of 1506 days. We note that this value is close to the true expansion speed because the inclination of the jet on the plane of the sky is ∼10° (see Paper II). Such a speed is about half that estimated for the NE peak in Paper II, of namely ∼440 km s⁻¹, which suggests that the expansion has indeed been slowed down more on the SW side than on the NE side. This evidence supports the idea that the SW lobe was also powered by the same outburst as the NE lobe, but was impeded at the beginning of its expansion.

Fig. 3. ALMA maps of the 92.2 GHz emission (colour scale, saturated to emphasise the emission from the SW lobe), after subtracting the dust contribution (see text). The contours are the map of the 22.2 GHz continuum emission obtained in 2023. Contour levels range from 0.34 to 2.72 in steps of 0.34 mJy beam−1. The ellipses in the bottom left are the half-power widths of the synthesised beams. a. Obtained from the ALMA data acquired on June 6, 2019 (see Paper II). b. Same as top panel, but for the data obtained on September 3, 2021.

3.2. Evolution of the jet emission

In Table 1, we give the flux densities of the free-free continuum emission from each jet lobe obtained from the new VLA data. We also list the ‘pure’ free-free fluxes at 92.2 GHz computed as previously explained from the ALMA data presented in Paper II. We used the flux densities in this table as well as those already obtained in Paper II to produce the light curves of both lobes, which are shown in Fig. 4. At all frequencies, we see that the free-free emission from the NE lobe (circles connected by solid lines) has decreased by a factor of about 3 at the last epoch with respect to the previous measurement. For the SW lobe (triangles and dashed line), we have single-epoch data at centimetre wavelengths, but the previous two measurements at 3 mm indicate a slight increase in the flux density. To compare the emission from the SW lobe in 2021 with that in 2023 is not trivial because the two measurements available are taken at different frequencies. However, it is very reasonable to assume that the 92.2 GHz flux is close to that measured at 22.2 GHz, because at the last epochs of the monitoring in Paper II (Fig. 4) the spectrum appears to be approximately flat above 22.2 GHz and the turn-off frequency ν_m (see Fig. A.1 of Paper I) appears to decrease with time. It is therefore likely that the flux from the SW lobe has dropped by a factor of about 2 after year 2021.

Fig. 4. Plot of the free-free continuum emission from the radio jet at 6, 10, and 22.2 GHz (coloured curves) as a function of time measured from the onset of the radio outburst on July 10, 2016 (see Paper I). Circles and solid lines indicate the emission from the NE lobe, while triangles and the dashed line denote the emission from the SW lobe. The error bar in the bottom right corresponds to an uncertainty of 15%. The black curve with error bars is the light curve of the accretion outburst obtained at K band in the near-infrared (see Fedriani et al. 2023).

We conclude that the jet activity has been weakening since 2021, as expected given that the energy input provided by the accretion event has come to an end. This can be appreciated in Fig. 4 by comparing the radio light curves with the IR light curve of the accretion outburst obtained from Table B.1 of Fedriani et al. (2023). From this comparison, the time lag of ∼600 days between the peaks of the two curves clearly stands out, and is roughly consistent with the delay of ∼400 days (see Paper I) between the onset of the NIR burst and that of the radio burst. Also, both curves are characterised by a rapid increase and a slow decline, but the timescales of the IR light curve are shorter, as expected for radiative processes (the IR burst) as opposed to hydrodynamical processes (the jet expansion).

3.3. Modelling the spatial and spectral emission from the jet

With the new VLA data, for the first time since the onset of the radio outburst in S255IR NIRS 3, we can both resolve the jet structure and measure the continuum flux at three frequencies. We can therefore fit both the jet morphology and the spectrum at the same time. For this purpose, we adopted the same model described in Paper II, namely a slightly modified version of the Reynolds (1986) model. In summary, this consists of a conical jet with opening angle θ₀ and vertex coinciding with the position of the star. The ionised gas temperature, T₀, ionisation fraction, x₀, and internal velocity, v₀, do not depend on the distance from the star, r, whose projection on the plane of the sky is y = r cosψ, with ψ = 10° being the inclination of the jet on the plane of the sky (see Paper II). The gas number density is n = n₀(r₀/r)². The ionised gas is confined between an inner radius r₀ and an outer radius r_m. All parameters but T₀ may change with time.

The model fit to the 22.2 GHz map of the jet was performed assuming that the whole jet is optically thin. This implies that the brightness temperature is T_B = T₀τ(y)=T₀τ₀(r₀/y)³ (see Eq. (B.3) of Paper II) and one can fit the map normalised with respect to the peak intensity ${\stackrel{\sim }{T}}_{\mathrm{B}}\equiv T_{\mathrm{B}}/T_{\mathrm{B}}({\mathit{y}}_0)={({\mathit{y}}_0/\mathit{y})}^3$. In this way, only three free parameters are left, that is, θ₀, y₀, and y_m. The best fit is obtained, after convolving the model to the synthesised beam of the observations, by minimising the expression ${\displaystyle \sum _i}{({\stackrel{\sim }{T}}_{\mathrm{B}}^i(\mathrm{model})-{\stackrel{\sim }{T}}_{\mathrm{B}}^i(\mathrm{data}))}^2$, where i is a generic pixel of the image. We then fix θ₀, y₀, and y_m to the best-fit values and vary τ₀ to minimise Eq. (6) of Paper II.

The fit to the spectrum depends on four parameters, θ₀, y₀, y_m, and Λ, where the latter is defined as Λ = x₀ Ṁ/v₀ with Ṁ being the mass-loss rate of the jet. It is worth noting that Λ can be expressed as a function of τ₀ through Eqs. (B.6) and (B.11) of Paper II. The best fit is obtained by minimising Eq. (10) of Paper I. As the two lobes evolve in a different way, it is convenient to discuss the fit to each of them separately.

3.3.1. NE lobe

As a first step, we fitted the morphology of the NE lobe as explained above. Figure 5a shows the 22.2 GHz map of the NE lobe rotated by 42° counterclockwise for the sake of comparison with the model fit in Fig. 5b. The corresponding best-fit parameters are given in Table 2.

Fig. 5. Comparison between observed and model maps of the 22.2 GHz continuum emission from the jet lobes. (a) Map of the NE lobe of the jet, rotated by 42° counterclockwise. Contour levels range from 10% to 90% in steps of 10% of the peak emission. (b) Best-fit model of the NE lobe. Contour levels are the same as for the data. (c) Same as panel a, but for the SW lobe, rotated by 48° counterclockwise. (d) Same as panel b, but for the SW lobe.

The second step is the fit to the spectrum. As only three data points are available, we prefer to reduce the number of free parameters by fixing y₀ to the value obtained from the fit to the map, 329 au, while y_m is computed from Eq. (7) of Paper II. The best fit is obtained for the values of θ₀ and Λ in Table 2 and is shown in Fig. 6a. We note that the values of θ₀ and Λ are in good agreement with those derived from the fit to the map, which lends support to the results obtained. We stress that assuming y_m = 760 au, as obtained from the fit to the map, would change θ₀ and Λ by only ∼10%.

Fig. 6. Spectra of the continuum emission from the jet lobes. The vertical bars correspond to a 15% error on the flux density, whereas the horizontal bars indicate the bandwidth covered at each frequency. The curves are the best fits obtained with the model described in the text. a. NE lobe. b. SW lobe.

Figure 7 shows the best-fit values of θ₀ and r₀, and the ionised jet mass computed from Eq. (B.13) of Paper II as a function of time for all the data already presented in Paper II plus the new data analysed in the present study. The conclusion is that the inner radius keeps expanding, while the jet recollimates, as indicated by the dramatic decrease in θ₀. Furthermore, the ionised gas recombines, because the decrease in x₀M_jet cannot be due to a decrease in M_jet, which remains constant after the inner radius begins to grow and the jet feeding stops (see also Paper II).

Fig. 7. Parameters obtained from the best fit to the continuum spectrum of the NE lobe as a function of time. (a) Ionised mass of the jet. (b) Inner radius of the NE lobe. The red point was obtained by fitting the map of the 22.2 GHz emission (Fig. 5b). (c) Opening angle of the jet.

3.3.2. SW lobe

The model fit to the map and spectrum of the SW lobe was performed in the same way as for the NE lobe, with the only difference being that we do not have an expression analogous to Eq. (7) of Paper II to calculate y_m as a function of time. For this reason, when fitting the spectrum, we assumed y_m equal to the value obtained from the fit to the 22.2 GHz map. The best-fit parameters are given in Table 2. Also in this case, as for the NE lobe, there is good agreement between the values of θ₀ and Λ estimated from the map and those derived from the spectrum.

While the opening angle and Λ are similar to those of the NE lobe, the inner and outer radii are significantly smaller, a result consistent with the SW lobe expansion being hindered by the environment. The ionised jet mass is $\sim (6_{-2}^{+3})\times {10}^{-6}$M_⊙, which is only marginally less than that of the other lobe. In summary, the scenario depicted by these results is that of a lobe whose expansion is delayed, but with physical properties that are basically the same as those characterising the NE lobe. We conclude that both lobes are powered by the same accretion event despite the different size, expansion speed, and morphology. Such a difference can be explained either by an asymmetry in the way the central engine feeds the two lobes, or by an interaction with the surrounding material, which may be denser on the SW side (as suggested in Paper II). On the basis of the data available to us, it is impossible to discriminate between these two hypotheses, but further high-resolution monitoring of the jet might be helpful.

4. Summary and conclusions

We used the VLA to image the ionised jet from S255IR NIRS 3 at three bands with subarcsecond resolution in order to complement our previous multi-epoch observations of the source (Papers I and II). The results confirm the appearance of a SW lobe, delayed with respect to the NE lobe, and allow us to trace its expansion and establish that the whole jet is fading away. The two lobes have different sizes, expansion speeds, and structures, hinting at an asymmetry in either the engine powering them or the surrounding material. We used the model developed in our previous papers to fit both the morphology of the lobes and their continuum spectra between 6 GHz and 22.2 GHz. Our findings confirm a common origin for the two lobes, both powered by the same accretion outburst studied by Caratti o Garatti et al. (2017). We also find that the radio jet is recollimating and recombining.

Acknowledgments

This study is based on observations made under project 23A-021 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. A.C.G. acknowledges support 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 FJC2021-046802-I, PID2020-114461GB-I00 and EX2021-001131-S funded by MCIN/AEI/10.13039/501100011033 and by “European Union NextGenerationEU/PRTR”.

References

All Tables

All Figures

Fig. 1. Maps of the continuum emission from S255IR NIRS 3 obtained with the VLA. (a) Data acquired at 22.2 GHz on July 21, 2023. Contour levels are quantified by the colour scale to the right of the panel. The ellipse in the bottom right is the half-power width of the synthesised beam. The cross marks the peak of the 3 mm continuum emission from Paper II, which we assume to be the position of the star. (b) Same as top panel, but for the 10 GHz emission. c. Same as top panel, but for the 6 GHz emission.

In the text

	Fig. 2. Map of the spectral index obtained from the ratio between the 6 GHz and 22.2 GHz images. The contours are a map of the 22.2 GHz continuum emission obtained with the same clean beam as the 6 GHz image. Contour levels range from 0.5 to 5.5 in steps of 1 mJy beam−1. The ellipse in the bottom right is the half-power width of the synthesised beam.
In the text

Fig. 3. ALMA maps of the 92.2 GHz emission (colour scale, saturated to emphasise the emission from the SW lobe), after subtracting the dust contribution (see text). The contours are the map of the 22.2 GHz continuum emission obtained in 2023. Contour levels range from 0.34 to 2.72 in steps of 0.34 mJy beam−1. The ellipses in the bottom left are the half-power widths of the synthesised beams. a. Obtained from the ALMA data acquired on June 6, 2019 (see Paper II). b. Same as top panel, but for the data obtained on September 3, 2021.

In the text

Fig. 4. Plot of the free-free continuum emission from the radio jet at 6, 10, and 22.2 GHz (coloured curves) as a function of time measured from the onset of the radio outburst on July 10, 2016 (see Paper I). Circles and solid lines indicate the emission from the NE lobe, while triangles and the dashed line denote the emission from the SW lobe. The error bar in the bottom right corresponds to an uncertainty of 15%. The black curve with error bars is the light curve of the accretion outburst obtained at K band in the near-infrared (see Fedriani et al. 2023).

In the text

Fig. 5. Comparison between observed and model maps of the 22.2 GHz continuum emission from the jet lobes. (a) Map of the NE lobe of the jet, rotated by 42° counterclockwise. Contour levels range from 10% to 90% in steps of 10% of the peak emission. (b) Best-fit model of the NE lobe. Contour levels are the same as for the data. (c) Same as panel a, but for the SW lobe, rotated by 48° counterclockwise. (d) Same as panel b, but for the SW lobe.

In the text

	Fig. 6. Spectra of the continuum emission from the jet lobes. The vertical bars correspond to a 15% error on the flux density, whereas the horizontal bars indicate the bandwidth covered at each frequency. The curves are the best fits obtained with the model described in the text. a. NE lobe. b. SW lobe.
In the text

	Fig. 7. Parameters obtained from the best fit to the continuum spectrum of the NE lobe as a function of time. (a) Ionised mass of the jet. (b) Inner radius of the NE lobe. The red point was obtained by fitting the map of the 22.2 GHz emission (Fig. 5b). (c) Opening angle of the jet.
In the text