Measuring the ionisation fraction in a jet from a massive protostar

It is important to determine if massive stars form via disc accretion, like their low-mass counterparts. Theory and observation indicate that protostellar jets are a natural consequence of accretion discs and are likely to be crucial for removing angular momentum during the collapse. However, massive protostars are typically rarer, more distant and more dust enshrouded, making observational studies of their jets more challenging. A fundamental question is whether the degree of ionisation in jets is similar across the mass spectrum. Here we determine an ionisation fraction of ~5–12% in the jet from the massive protostar G35.20-0.74N, based on spatially coincident infrared and radio emission. This is similar to the values found in jets from lower-mass young stars, implying a unified mechanism of shock ionisation applies in jets across most of the protostellar mass spectrum, up to at least ~10 solar masses.


Introduction
Although massive stars are rare, they play a fundamental role in the Universe, synthesising most of the chemical elements and providing a major feedback into the molecular clouds where stars are born (see ref. 1 and references therein for a review). However, massive star formation and evolution are still a matter of debate. An excellent tool to investigate massive star formation is provided by protostellar jets. The detection of several jets driven by massive protostars [2][3][4][5] and the discovery of dusty molecular discs around high-mass young stellar objects (HMYSOs, M * > 8 M , L bol > 5 × 10 3 L ) through near-infrared (NIR) interferometry 6 , strongly support the idea that HMYSOs form in a similar way to low-mass young stars [6][7][8][9][10] . Collimated jets associated with discs of few 100 au have also been observed at mm wavelengths, supporting this scenario 11,12 . The jets themselves are thought to be launched centrifugally along magnetic field lines 13 and thus the jet's ionisation fraction (χ e ) is a key parameter as determines the strength of the coupling of the magnetic field to the ionised gas. Often, when irradiated by a nearby OB star, the jet is considered to be fully ionised or alternatively a somewhat arbitrary value is used [14][15][16] . In a number of cases for low-mass young stars, χ e has been measured using emission line diagnostics in the optical/NIR regime 17 . Determining χ e is important if we are to correctly deduce such parameters as the total (neutral plus ionised) mass-loss rate in the jet (Ṁ ejec ), which, in turn, is often compared with the mass-accretion rate (Ṁ acc ). Since dynamical quantities are fundamental inputs in massive star formation models, a good constraint on the ionisation fraction along massive jets is critical. The radio continuum emission from protostellar jets is generally interpreted as thermal bremsstrahlung, and, since the radio emission does not suffer from extinction due to dust, it is an excellent way of observing the ionised component 15,18,19 . However, no velocity information can be inferred from these observations (unless we rely on proper motion studies and we know the jet geometry) nor can the ionisation fraction be directly derived. On the other hand, NIR spectroscopy directly provides us with physical parameters and dynamical information for the (molecular and atomic) jet through the study of emission lines 2, 3 . Therefore, radio and NIR observations of both line and continuum emission allow for a complementary analysis of protostellar jets.
The well-known high-mass star-forming region G35.20-0.74 is located at 2.2 kpc 20 in the tail of the Aquila constellation. G35.20-0.74N (hereafter G35.2N) is a main formation site of Btype stars, has a bolometric luminosity of 3 × 10 4 L and hosts two main cores, core A and core B. Both cores display discs in Keplerian rotation 21,22 . Core B is a binary system which consists of two B-type protostars with masses of 11 and 6 M 22 , sources 8a and 8b, respectively. The rotation axis of the disc, i.e. the jet axis, of source 8a has an inclination angle of i ∼ 19 ± 1 • with respect to the plane of the sky 23 . Perpendicular to this disc a radio jet close to the central engine has been detected 22,24 as well as a molecular hydrogen (H 2 ) outflow 2,25 . Source 8a is thought to drive one of two parsec-scale bipolar outflows 2,22 in this region with an initial north (blue-shifted lobe) south (red-shifted lobe) orientation (see Fig. 3 of ref. 2 for a complete view of both parsecscale jets). Here, we focus on the protostellar jet driven by source 8a (see Fig. 1).
Here, we show a unique example of spatially coincident NIR and radio jet emission. We use multi-wavelength observations, i.e. data from the Hubble Space Telescope (HST), Karl G. Jansky Very Large Array (VLA), and Very Large Telescope (VLT), of the outflow from G35.20-0.74N to determine the ionisation fraction, χ e , in a jet from a massive young star. The values found, ∼ 5−12%, are similar to those found for solar mass young stars. Our observations confirm that the ionising mechanism giving rise to the radio emission originates from shocks seen in the NIR jet 18 .

Results
NIR imaging and spectroscopy in G35. 2N We performed high-resolution NIR imaging and long-slit spectroscopy of G35.2N using the infrared spectrometer and array camera (ISAAC) at the VLT. We performed imaging to study in detail the NIR jet and spectroscopy to measure its kinematics and derive its physical and dynamic properties. We combine our ob-servations with VLA and HST data to form the most complete view of a jet close to its source to date (see Fig. 1, panel b). We observed the atomic jet in the form of [Fe II] and Brγ (hydrogen recombination line) in the NIR as well as the ionised jet in the form of H II (ionised hydrogen) in the radio. G35.2N, thus, represents a unique example because both the NIR and radio jet are visible and we see that both atomic and ionised emission are spatially coincident. Therefore, we can combine the information from both regimes and infer the ionisation fraction in a HMYSO jet. White contours are ISAAC H 2 − K emission (corresponding to the H 2 emission at 2.121 µm), which mostly delineate the outflow cavity walls 2 . Green contours are VLA C-band emission at 6 cm (5.8 GHz), that comes from the ionised radio jet. Notably, the radio data, not affected by extinction, trace both the blueand red-shifted lobes, whereas the NIR emission mostly comes from the blue-shifted jet. The position of the radio source driving the jet is labelled as Core B and the main knots of study are labelled as K1, K2, K3, and K4. The cyan contours show the Atacama Large Millimeter Array (ALMA) 870 µm emission from dust, which reveals the locii where dust and gas condense into stars 23 .
The spectral images of the protostellar jet close to the source are shown in the position-velocity (PV) diagram (Fig. 2). Our spectroscopic slit, at a position angle (PA) of −14.7 • , encompasses knots K2, K3, and K4. The velocities are measured with respect to the Local Standard of Rest (LSR) and subsequently corrected by the velocity of the parent cloud (v LSR ∼ 33 km s −1 , ref. 21,23 ). H 2 emission, namely the 1 − 0 S(1) transition at 2.121 µm and the 2 − 1 S(2) transition at 2.154 µm, is observed towards the outflow (see Fig. 2, panel f for the spectral image of the 1 − 0 S(1) line). The radial velocity (v rad ) of the H 2 lines peaks at red-shifted velocities of ≈ 5 km s −1 , corresponding to a total velocity of v tot ≈ 15 km s −1 (v tot = v rad / sin i, where i is the inclination of the jet axis with respect to the sky), which might be an indication of slow oblique shocks against the cavity walls. [Fe II] emission is also detected along the jet (see Fig. 2, panels b, c, and d). Blue-shifted radial velocities range from few km s −1 to −200 km s −1 , which correspond to a total velocity of v tot ≈ −600 km s −1 . We also observe Brγ emission at 2.166µm (hydrogen recombination line transition from n up = 7 to n low = 4). This emission is observed to be extended up to ∼ 8 from the central source (i.e. ∼ 17 600 au at a distance of 2.2 kpc; see Fig. 2, panel e). This emission covers the same velocity range of the [Fe II] emission, consistent with proper motions given by ref. 22 . Emission from [Fe II] and Brγ is, thus, indicative of shocked material 3,27 . This atomic emission is spatially coincident with the radio jet emission towards the north (blue-shifted lobe) where the visual extinction (A V ≈ 25 mag 28     tection of NIR lines. Conversely, towards the south (red-shifted lobe), the visual extinction reaches up to ∼ 41 mag (increasing up to ∼ 170 mag towards the driving source and its immediate surroundings) 28 hindering any NIR jet-detection, although the redshifted radio jet is still visible (Fig. 1, panel b). To investigate the possibility that some of the Brγ emission is due to scattered light, we have thoroughly checked the PV diagram of Figure 2. Brγ emission close to 0 km s −1 is likely scattered emission from the central source. But, there is also material moving radially at more than −150 km s −1 (which corresponds to v tot ≥ −500 km s −1 ). Moreover, the full width at zero intensity (FWZI), which gives an estimate of the jet shock velocity 29 , is larger than 500 km s −1 . This value is consistent with the total velocities derived above, supporting the jet geometry. Finally, the [Fe II] emission spatially coincides with that from Brγ and their velocities are similar. This evidence suggests that the Brγ line is likely tracing shocked material as [Fe II] does. This is also confirmed by the fact that all Brγ emission is blue-shifted which is consistent with the lobe being directed towards us. Thus, the Brγ emission likely comes from a combination of scattered light and from directly visible shocked material. Notably, the radio continuum, which traces ionised gas, emits co-spatially with the atomic emitting region implying that the ejected material is partially ionised.

Determination of ionisation fraction and dynamic properties
We can determine the degree of ionisation of the jet from the massive protostar by directly comparing the total number density (n tot ) with the electron number density (n e ). The ionisation fraction is defined as χ e = n e /n tot . On the one hand, the n tot can be derived from the properties of the [Fe II] emission line at 1.644 µm (among the brightest transitions in the NIR). The product of the total number density times the volume emitting region (n tot V ) can be written as the ratio between the observed luminosity of the line and the emissivity per particle for the line calculated theoretically 17  On the other hand, the n e can be derived from the ratio of the [Fe II] lines 26,30 or from the properties of the radio emission 31 (see Methods). Interestingly, both independent methods provide us with the same electron density estimates within the errors. As the radio regime here provides us with smaller errors, we then adopt these electron density estimates. We are able to give a direct measurement of the ionisation fraction in a HMYSO jet for knots K1 (12±6%), K3 (7±1%), and K4 (5±1%). Additionally, in the case of K2 (< 17%) we are able to set an upper limit due to the non-detection of radio emission. Figure 3 shows a plot of the variation of the ionisation fraction with distance where there is not a strong indication that the χ e changes with distance as the inferred values are the same within the error bars. It is worth noting that the uncertainty in the knot sizes does not significantly affect the ionisation fraction measurement of knots K3 and K4, whereas in the case of K1, the uncertainty in the knot size represents half of the error of the ionisation fraction. Finally, we would like to stress that the determination of the ionisation fraction rests under a few reasonable assumptions, namely i) all Fe is ionised, ii) the Fe abundance is solar, and iii) there is no Fe dust depletion (see Methods for a detailed discussion of these assumptions).
Now, using the mass of each knot, its velocity, and its length, we can also determine the total mass-loss rate using the [Fe II] emission line at 1.644 µm (see Methods). We obtain a mass-loss rate of 4.7 ± 0.3 × 10 −5 M yr −1 and 4.9 ± 0.4 × 10 −5 M yr −1 for knots K3 and K4, respectively. In the case of knot K1 there is no velocity information, because our spectroscopic slit does not encompass it, and in the case of knot K2 just an upper limit on n e could be inferred, therefore no dynamic properties were derived. The momentum rate is 0.018 ± 0.001 M yr −1 km s −1 and 0.011 ± 0.001 M yr −1 km s −1 for knots K3 and K4, respectively.
The total mass-loss rate can also be calculated using the properties of radio emission for these knots in combination with the ionisation fraction derived above. In an analogous way, we combine the size of the radio knots, the electron density, the tangential velocity from our spectroscopic observations (i.e. v tan = v rad / tan(i)), and the ionisation fraction, to calculate the total mass-loss rate for knots K3 and K4 (see Methods). It results inṀ ejec = 1.6 ± 1.0 × 10 −5 M yr −1 for knot K3 anḋ M ejec = 2.0 ± 1.0 × 10 −5 M yr −1 for knot K4. Remarkably, these values match very well, within a factor of 2 to 4, with those obtained in the NIR regime. Such a small variation is caused by the slight difference in the size of the NIR and radio emitting regions, with the latter being smaller. This is not unexpected as the jet is supposed to have an onion-like structure (see e.g., ref. 32,33 ), and the more ionised component is expected to be confined to the innermost region (i.e., closest to the jet axis). The ionised mass-loss rate (Ṁ ionised ) of the considered radio knots is ∼ 10 −6 M yr −1 being a factor of ∼ 10 smaller than what we get in the NIR, as expected from the ionisation fraction values obtained here. This is in excellent agreement with the ionised mass-loss rate calculated on source (core B) following Reynolds' formulation 34 , which provides a value of 1.8 ± 0.3 × 10 −6 M yr −1 , consistent with ref. 22 (see Methods).
The mass-accretion rate can also be estimated assuming that the ratioṀ ejec /Ṁ acc is between ∼ 10% and ∼ 30% 35,36 . Considering a mass-loss rate of ∼ 4 − 6 × 10 −5 M yr −1 , we obtain a mass-accretion rate of 1.3 × 10 −4 Ṁ acc 6 × 10 −4 M yr −1 consistent with ref. 37,38 . A typical jet phase timescale in HMYSOs is ∼ 2 × 10 4 yr 3, 39 , which gives a lower limit of the age of the star, and, in turn, an estimate of the mass of the central source being 2.6 M * 12 M , which agrees with previous estimates using other methods 22,23,37 .

Discussion
In the case of jets from low-mass young stars, the ionisation fraction is typically found to be around 10% or less 27,40 . Recent studies of two jets from HMYSOs, tentatively derive similar values in an indirect manner and under different assumptions than used here (∼ 10% in IRAS13481-6124 3 , and 14% in S255IR NIRS3 10 ). In the case of IRAS13481-6124, the authors assume that the mass-loss rate calculated along the parsec-scale jet, traced by the H 2 , has remained roughly constant for the formation history of the star and, thus, should be similar very close to the star where the actual mass-loss rate could not be calculated. Then, they compare this value with the ionised mass-loss rate inferred from the radio regime very close to the star. With this comparison of mass-loss rates, they give a rough estimate of the ionisation fraction. In the case of S255 NIRS3, the authors estimate the mass-accretion rate from the NIR regime and assume that at least 10% of the accreted material should be ejected. This provides an upper limit to the ionisation fraction by comparing to the ionised mass-loss rate calculated from the radio continuum emission. In any case, these previous estimates were made without a direct evidence of spatially coincident atomic and ionised jet emission. In any event, our result is consistent with these indirect findings, i.e., that the ionisation fraction is low and similar to what is found in the low-mass regime. Additional studies of this kind however are warranted to determine whether this is the norm amongst outflows from HMYSOs.
In the regions of the jet of G35.2N probed by our observations -i.e. knots K1, K2, K3 and K4-we see a similar ionisation degree as found in low-mass YSOs, suggesting a mechanism of shock ionisation, rather than photo-ionisation, for setting the ionisation degree. However, towards the later stages of massive star formation the protostar is expected to contract close to a zero age main sequence (ZAMS) structure and emit copious UV radiation, which would then set a higher ionisation degree, closer to unity 41 . This stage might not have been reached yet by G35.2N and at current high accretion rates, dust would quench the full ionisation of the disc surface, therefore keeping its wind mostly neutral. Another possible scenario to explain the low χ e is that if the central engine is accreting at high rates, then the protostar is still swollen and too cool to photo-ionise much. Both scenarios are in agreement with the fact that the ionised mass-loss rate on source is very similar to the ionised mass-loss rate on the studied knots. This implies that the χ e should not be much higher very close to the star or that the total mass-loss rate has dramatically changed (1 order of magnitude) in the past few 100 years with respect to the knots K3 and K4. This reasoning also strengths the idea that the UC/HC H II region, if present, must be very small and confined within 100 au 22 . Anyway, in extreme cases, some jets from HMYSOs are even observed to emerge from hypercompact H II regions 42 . Our results do indicate that ionising photons from the protostar are confined to a small region, which constrains models of both evolutionary state and feedback of this source. Thus, application of these methods to a larger sample can help probe these processes more generally during massive star formation.
Finally the wider implications of our findings should be mentioned. The high collimation, velocities and momentum injection rate observed in HMYSO jets suggest that they must be (magneto-)centrifugally launched from very deep in the star's gravitational well. Coupling of magnetic fields to the gas requires a sufficient level of ionisation: the values we observe are high enough for efficient coupling 13 and are an important constraint for theoretical models of these jets and outflows 43,44 , which are likely to be the dominant feedback mechanism even in massive protostellar cores 45 . The measurements that we present here are relevant to conditions at locations that are away from the disc and jet launching region. These results are important for the dynamics of outflows at distances of few 1000 to few 10 000 au from the protostar, where processes of jet collimation and outflowcore interaction are likely to be occurring. However, observations at higher spatial and spectral resolution would be required to probe the immediate jet environment close to the young star, where the jet launching mechanism could be revealed. Simulations at high resolution have been able to differentiate between a magneto-centrifugally launched highly collimated jet and a slow wide angle magnetic-pressure driven tower flow 44 , but our observations cannot discriminate between them. If we accept that indeed jets from HMYSOs are magneto-centrifugally driven, then it should be noted that earlier magnetohydrodynamical (MHD) simulations estimate lower jet velocities than those observed here, at small scales [46][47][48] and at large scales 49 . However, more recent MHD simulations obtained velocities of several 100 km s −1 for their magneto-centrifugal jets 44 , consistent with our results. In terms of mass ejection and momentum rates, MHD simulations reproduce a wide range of values depending on the magnetic and rotational energies 44,48 , being difficult to discriminate between the different simulations. Ref. 43 provides dynamical estimates remarkably close to our observations, in particular, the authors obtain mass-loss and momentum rates of the order of ∼ 10 −5 M yr −1 and ∼ 10 −2 M yr −1 km s −1 , respectively, for an 8 M protostar. Nevertheless, simulated outflow properties are highly dependent on resolution (see ref. 44 for an overview of convergence aspects and numerical resolution problems in previous simulation works).

Methods
NIR VLT/ISAAC imaging and spectroscopy NIR imaging was obtained on 2013 July 7 in the K and H 2 bands with the ISAAC instrument on VLT (ESO, Chile). This observational strategy allows us to subtract the continuum emission of the nebulosity to obtain the H 2 emission line of the system (see Fig. 1 panel b, white contours). We correct the narrow band image (H 2 ) by subtracting the continuum in the broad band image (K), after first having normalised using the counts of the field stars. Additionally, high-resolution long-slit spectra were obtained with VLT/ISAAC on 2013 September 14 in the H and K bands. The 0.3 × 120 slit was positioned at the brightest point in the NIR with a position angle of −14.47 • (see Fig. 2 panel a, green lines). The spatial sampling was 146 mas pixel −1 and the spectral resolution was R ∼ 10 000 and 8 900, corresponding to a velocity resolution of 30 − 35 km s −1 , for the H and K bands, respectively. Total integration time was 1440 s for each band. The spatial resolution was seeinglimited at ∼ 0.8 − 1 . The standard ABBA nodding technique was applied and the data were reduced in the standard way using the Image Reduction and Analysis Facility software (IRAF, http://ast.noao.edu/data/software).  Figure 1 has been generated in a similar way to the H 2 -K image explained in the previous section.

Radio VLA imaging
VLA observations were carried out on 2012 December 29, at 5.8 GHz, with the VLA in the A-configuration (θ beam ∼ 0.3 ) and utilising a total bandwidth of 2 GHz. Flux/bandpass and gain calibrators used were 3C286 and J1824+1044, respectively. Data reduction was performed in the usual way using the CASA software package (https://casa.nrao.edu/) together with the CASA pipeline (version 4.7.2). After several rounds of phase-only self-calibration, an RMS noise level in the image of σ = 6.4 µJy beam −1 was achieved. Flux densities can be found in Table 1, which are taken from 50 . Observations were conducted under program ID 12B-140, PI: M. Hoare.

Electron density
The electron density can be estimated through the ratio of the [Fe II] lines using a non-local thermodynamic equilibrium model (NLTE) that considers the 16 fine-structure levels 26,30 . In our H band spectroscopic ISAAC observations we detect three [Fe II] lines, the so-called a 4 D 7/2 − a 4 F 9/2 at 1.644 µm, a 4 D 1/2 − a 4 F 5/2 at 1.664 µm, and a 4 D 5/2 − a 4 F 7/2 at 1.677 µm (see Ta-ble 2 for the fluxes of the various knots). The ratio 1.664/1.644 provides a good estimate of the electron density whereas the ratio between 1.644/1.677 provides a less accurate estimate because this ratio is less sensitive to the electron density at values higher than few 10 4 cm −3 . The electron density can also be computed using the radio properties of the knots following Equation (7) of ref. 31 : where u 1 = 0.857 is the density model II (cylinder) conversion factor for computing electron density, T e is the electron temperature, S 5 GHz is the flux density at 5 GHz, D is the distance to the source, and θ G = (θ min · θ maj ) 0.5 is the deconvolved Gaussian width of the knot. This equation is valid if the emission is optically thin, which is the case of the knots under consideration (see previous section). Using Equation 1 and considering T e = 10 4 K, D = 2.2 kpc, the flux density (from Table 1, assuming α = −0.1), and the size of each knot (from Table 5); values of the order of 10 4 cm −3 were obtained for the various knots (see Table 3, Column 2). As we obtain spectral indices of ∼ −0.3 in our knots, there might be some mixture between thermal and non-thermal emission. To remove any non-thermal contribution, we assume α = −0.1 and extrapolate the flux density from the 23 GHz data. Anyhow it is worth noting that the flux dependence is to the power of 0.5 and, therefore, there is no significant change (less than a factor of 1.5) in the electron number estimates.
It is worth mentioning that in the case of knot K3, where we are able to detect the [Fe II] lines with enough S/N ratio (see Table 2) to estimate a reliable electron density, the estimates for the electron density from the NIR and the radio coincide within the errors (see Table 3). The same applies for knot K4, but here the uncertainties are much larger because of the low S/N ratio (see Table 2). Therefore, the two independent methods provide us with the same result. In the case of knot K4, however, the S/N ratio in the [Fe II] lines is quite low to estimate a reliable electron density from the NIR as is evidenced by the larger errors in this knot. We could still calculate a more precise electron density with smaller errors from the radio regime. For knot K1 no information about the electron density could be retrieved from the NIR because our slit did not encompass this knot and only radio emission could be used to calculate the electron density. Finally, for knot K2, only upper limits for the electron density could be given because the fluxes of the lines 1.664 µm and 1.677 µm are below 3σ (see Table 2). For consistency, we then adopt the electron density derived from the radio for the four knots.

Ionisation fraction
In the following the derivation of the ionisation fraction is outlined (the justification for our assumptions is given below). The ionisation fraction is defined as Equation 2 : where n tot is the total number density and n e is the electron number density. On the one hand, the number electron density is calculated following the previous section. On the other hand, the total number density can be derived following the next steps. Firstly, the product n tot V , where V is the volume emitting region, can be expressed as the ratio between the observed luminosity of the line and the emissivity per particle for the line calculated theoretically 17 . That is, where A i and f i are the radiative rates and fractional population of the upper level of the considered transition and F e + F e is the ionisation fraction of the iron having a total abundance with respect to hydrogen of F e H 27, 51 . Here, we have used the [Fe II] line at 1.644 µm. We have assumed that all iron is ionised and solar abundance of 2.8 × 10 −5 under the hypothesis of no dust depletion 52 . Secondly, both the radius and the length of the knots are resolved in our HST and VLT/ISAAC, therefore we can calculate the precise value of the volume. Considering that the volume of the knot is a cylinder with radius r jet = D · tan(θ min /2) D · θ min /2, where θ min is the deconvolved diameter of the knot (see Fig. 4 and Table 4 Column 6), and with length l ⊥ = D ·tan(θ maj ) D ·θ maj where θ maj is the deconvolved length of the knot (see Table 4 Column 5), we can write V = πr 2 jet l ⊥ , being D the distance to the source. Finally, solving for n tot , we can derive the total number density of the region given by Equation 3: From a theoretical point of view, iron is easily ionised (E ion = 7.9 eV) in mild J-shocks (v shock ∼ 25 − 50 km s −1 , see e.g., ref. 53,54 ); in the case of G35.2N we are dealing with much faster shocks (v shock ≥ 500 km s −1 ). The assumption that all iron is ionised is also supported by our observations where both [Fe II] and ionised hydrogen (E ion = 13.6 eV) emissions are co-spatial. Moreover, no [Fe I] emission has ever been detected in protostellar jets driven by HMYSOs. Indeed, [Fe I] has seldom been observed in protostellar jets driven by low-mass YSOs. For example, ref. 55  is observed. The assumption of no dust depletion is sustained by the fact that the faster the shock velocity, the higher the fraction of the dust-forming elements transferred into the gas-phase 56 . Indeed, shocks with velocities ≥ 400 km s −1 completely destroy the dust 57 and release the Fe grains into gas-phase. We certainly measure jet velocities > 500 km s −1 .
On the other hand, it should be noted that ref. 17 found that Fe depletion can reach up to ∼ 90% for the HH 34 jet, driven by a low-mass YSO. However, the jet we are analysing here is driven by a HMYSO, and therefore more energetic shocks are measured. Indeed, the line profile of less energetic jets as HH 34 is very narrow up to large distances from the source (see e.g. Fig. 3 of ref. 27 ). In the case of our Figure 2 panel b, it is clear that the line profile is different, resembling that of a bow-shock. In addition, later studies on HH 34 and other low-mass jets (see ref. 58 and ref. 59 ) show that [Fe II] is mostly detected in the low-velocity component close to the source whereas the high-velocity component show very low or no dust depletion. In any case, if the scenario of high depletion were in place in our study, this would imply that n tot is underestimated by an order of magnitude and, therefore, that the ionisation fraction is also overestimated by an order of magnitude.
For knots K3 and K4 we are able to provide a good estimate of the ionisation fraction. However, for knot K2 only an upper limit for the electron density could be obtained and therefore an upper limit for the ionisation fraction is given. Finally, for knot K1 the uncertainty is quite large due to that the knot geometry is not very clear in our NIR images and we have considered larger errors.

Mass-loss and momentum rate
The mass-loss rate was computed following the equation: where M is the mass of the knot, v ⊥ the tangential velocity of the knot, and l ⊥ the length of the knot in the plane of the sky. The mass of the knot is not a direct measurement, but it can be written as M = µm H n tot V where µ = 1.24 is the mean atomic weight, m H = 1.67 × 10 −27 kg is the mass of the hydrogen, n tot is the total density of the hydrogen (both neutral and ionised) of the flow, and V is the volume of the emitting region (assumed to be a cylinder). We derived above that the product n tot V can be expressed as the ratio between the line luminosity and the line emissivity of the [Fe II] line at 1.644µm. Therefore, the Equation 4 can be written as Equation 5: The luminosity of the line is expressed as L = 4πD 2 F , where D is the distance to the source from Earth and F the dereddened flux using the Rieke & Lebofsky extinction law 60 considering A V = 25 mag 28 . The input parameters values and the mass-loss rates can be found in Table 4.
The mass-loss rate can also be calculated using the radio emission properties of the various knots. In this case, the mass of the knot is calculated as M = V · n e · µ/χ e , where V is the volume of the knot (considered to be a cylinder as in the case of the NIR), n e is the number electron density calculated in the previous section, µ is average atomic weight, and χ e is the ionisation fraction (see previous section). The ejection time, which is the period of time that the jet lobe was ejected over, is equal to τ ejec = D · tan(θ maj )/v ⊥ . See Table 5 Columns 2 and 3 for the deconvolved dimensions of the emission lobes. Finally, the massloss rate is then given byṀ ejec = M/τ ejec , see Table 5 for input parameters and mass-loss rates estimates (Column 5 for ionised mass-loss rate where χ e was simply considered equal to 1, Column 6 for total mass-loss rate considering the ionisation fraction calculated above).
The mass-loss rates for knots K1 and K2 are quite uncertain. In the case of K1, no reliable velocity could be measure; and in the case of K2, the electron density could not be retrieved, just an upper limit.
Finally, the momentum rate (or thrust) is given byṖ = M ejec v tot (see Column 9 of Table 4 for the momentum rate of NIR jet, and Column 7 of Table 5 for the momentum rate of the ionised radio jet,Ṗ ionised ).

Ionised mass-loss rate on source
The ionised mass-loss rate can be calculated following the Reynolds' formulation 34 (see also ref. 18 ). For a canonical jet (α = 0.6) the ionised mass-loss rate is given by Equation 6: where S ν is the flux density, ν is the frequency, v j is the jet velocity, θ 0 is the jet opening angle, i is the inclination of the jet, D is the distance, and T the temperature. From our observations we obtain the following parameters for core B: S 5.8 = 0.794 ± 0.03 mJy, ν = 5.8 GHz, v j = 600 ± 100 km s −1 , θ 0 = 52.3 • ± 4.4 • , D = 2.2 kpc, T = 10000 K, the resulting ionised mass-loss rate is 1.81 ± 0.33 × 10 −6 M yr −1 , consistent with ref. 22 . This value is very similar to those found along the knots (see Table 5, Column 5). This indicates that has not been a great variation in the ionised mass-loss rate between the knots located up to ∼ 18 000 au and the jet very close to the massive protostar. Table 4: Input parameters and mass-loss rates for the [Fe II] 1.644µm line. The source of uncertainties in Columns 4, 5, and 6 is coming from measurement uncertainties, whereas in Columns 2, 3, 7, 8, and 9 comes from performing the error propagation in the specific equations. 9.12 ± 0.33 3.1 ± 0.2 205 ± 10 1.8 ± 0.1 0.52 ± 0.01 4.5 ± 0.2 4.9 ± 0.3 1.1 ± 0.1 * Assuming a tangential velocity of 150 − 300 km s −1 from the proper motions of ref. 22 . Table 5: Input parameters and mass-loss rates for the Radio C-band (6 cm). The source of uncertainties in Columns 2 and 3 is coming from measurement uncertainties and in Columns 4, 5, and 6 performing the error propagation in the specific equations. · · · · · · · · · · · · · · · · · · K3 0.61 ± 0.03 0.23 ± 0.01 2.0 ± 0.5 1.1 ± 0.3 1.6 ± 0.4 3.9 ± 0.9 K4 0.93 ± 0.07 0.33 ± 0.03 4.8 ± 1.7 1.0 ± 0.3 2.0 ± 0.7 2.3 ± 0.7 * Assuming a tangential velocity of 150 − 300 km s −1 from the proper motions of ref. 22 .  Figure 4: Gaussian fitting to the jet width using the HST [Fe II] 1.644 µm continuum-subtracted image. Black data points are the data, red data points are the best Gaussian fit, and the blue dots are residuals from the fit. Best fit values for the jet width are indicated in the figure legend for the K2, K3, and K4 knots in panels a, b, and c, respectively.

Data availability
All data used in this study are public and can be accessed through the different data archives of the various instruments using the Program ID and/or PI name given in the section above. ESO archive: http://archive.eso.org/eso/ eso_archive_main.html, HST archive: http://hst. esac.esa.int/ehst/#home, VLA archive: https:// archive.nrao.edu/archive/advquery.jsp. Upon reasonable request the authors will provide all data supporting this study.