Characterizing the radio continuum nature of sources in the massive star-forming region W75N (B)

The massive star-forming region W75N~(B) is thought to host a cluster of massive protostars (VLA~1, VLA~2, and VLA~3) undergoing different evolutionary stages. In this work, we present radio continuum data with the highest sensitivity and angular resolution obtained to date in this region, using the VLA-A and covering a wide range of frequencies (4-48~GHz), which allowed us to study the morphology and the nature of the emission of the different radio continuum sources. We also performed complementary studies with multi-epoch VLA data and ALMA archive data at 1.3 mm wavelength. We find that VLA~1 is driving a thermal radio jet at scales of $\approx$0.1 arcsec ($\approx$130 au), but also shows signs of an incipient hyper-compact HII region at scales of $\lesssim$ 1 arcsec ($\lesssim$ 1300~au). VLA~3 is also driving a thermal radio jet at scales of a few tenths of arcsec (few hundred of au). We conclude that this jet is shock-exciting the radio continuum sources Bc and VLA~4 (obscured HH objects), which show proper motions moving outward from VLA~3 at velocities of $\approx$112--118~km/s. We have also detected three new weak radio continuum sources, two of them associated with millimeter continuum cores observed with ALMA, suggesting that these two sources are also embedded YSOs in this massive star-forming region.


INTRODUCTION
Although it is well-known that the most massive stars have a great impact on the galactic environment, many aspects related to their early evolutionary stages still remain unknown. For instance, massive protostars are deeply embedded in dense molecular gas, located at typical distances of few thousand parsecs. Thus, detailed stud-Contact e-mail: adriana.rodriguez@unc.edu.ar ies of these objects require observations with very high sensitivity and angular resolution. One of the best known massive star-forming regions is W75N (B), located in the Cygnus X complex at a distance of 1.3 kpc (Rygl et al. 2012), comprising dense molecular clouds (Dickel, Dickel & Wilson 1978;Persi, Tapia & Smith 2006) and showing strong maser emission of different molecular species (e.g., Baart et al. 1986;Hunter et al. 1994;Torrelles et al. 1997;Surcis et al. 2009;Krasnov et al. 2015;Colom et al. 2018). This region constitutes an excellent laboratory to study early stages of massive star formation, since it hosts a cluster of massive protostars (e.g., Shepherd, Testi & Stark 2003), probably undergoing different evolutionary phases (Torrelles et al. 1997).
Since its discovery, W75N(B) has been widely studied, revealing the presence of five radio continuum sources (named VLA 1, VLA 2, VLA 3, VLA 4, and Bc; e.g., Hunter et al. 1994;Torrelles et al. 1997;Carrasco-González et al. 2015) and a large-scale high-velocity molecular outflow (e.g., Davis et al. 1998;Shepherd, Testi & Stark 2003). Among the five radio continuum sources, VLA 1 was proposed to be an evolved young stellar object (YSO), whereas VLA 2 is probably the least evolved YSO in the region (e.g., Torrelles et al. 1997). These two sources are the only ones in the region that are associated with 22 GHz water (e.g., Torrelles et al. 1997;Surcis et al. 2009Surcis et al. , 2011Surcis et al. , 2014Kim et al. 2013) and 6.7 GHz methanol maser emission, which was actually detected from a location in between them (e.g., Minier, Booth & Conway 2000;Surcis et al. 2009). Furthermore, polarimetric maser observations show the presence of a magnetic field oriented in the direction of the molecular outflow (e.g., Hutawarakorn, Cohen & Brebner 2002;Surcis et al. 2009Surcis et al. , 2011Surcis et al. , 2014. However, despite the deep studies conducted so far towards W75N(B), the nature of some of the radio continuum sources in the region is not well known yet.
In this work, we analyze radio continuum data obtained with the Karl Jansky Very Large Array (VLA) over a wide range of frequencies (4 to 48 GHz), which provide images with the highest sensitivity (rms = 8 µJy/beam) and angular resolution (0. 12 × 0. 09, PA= -69 • ) obtained to date in this region. Part of these data were presented by Carrasco-González et al. (2015), who focused their attention on the remarkable source VLA 2, reporting through radio continuum and H 2 O maser observations the transition from an uncollimated outflow to a collimated outflow over a period of only 18 years. In this work, we focus on the remaining sources in the field: VLA 1, VLA 3, VLA 4, and Bc (see Fig. 1a). These observations allow us to perform a deep multifrequency study of the morphology of the sources and of their nature. We also analyze Atacama Large Millimeter Array (ALMA) 1.3 mm continuum and spectral line archive data obtained toward this region.

VLA
The star-forming region W75N (B) was observed with the VLA of the National Radio Astronomy Observatory (NRAO) 1 in its Aconfiguration at C (6 cm), Ku (2 cm), K (1.3 cm), and Q (7 mm) bands (project code 14A-007). A detailed description of the observations and calibration procedures can be found in Carrasco-González et al. (2015). Deconvolved images were obtained with the task CLEAN of the Common Astronomy Software Applications (CASA 2 , version 4.1.0) data reduction package, using multifrequency synthesis (parameter nterms = 2) and multi-scale cleaning (Rau & Cornwell 2011). Primary beam corrections were applied. We split each band data set to build images of narrower bandwidth (1 and 2 GHz), using different weightings, i.e., natural, uniform, and Briggs (Briggs 1995) to achieve the best compromise between sensitivity and angular resolution, depending on the analysis performed (e.g., spectral energy distributions, angular size vs. frequency). We 1 NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. 2 https://science.nrao.edu/facilities/vla/data-processing also made a single image combining all four bands (C, Ku, K, and Q; Fig. 1), as well as individual images integrating the full bandwidth of each frequency band (Fig. 2). Moreover, the multifrequency synthesis cleaning technique allows us to obtain a spectral index map covering the entire range of the observed frequencies.
All the radio continuum images presented in this paper, as well as the spectral energy distribution analysis of the different sources are based on the VLA project code 14A-007 (epoch 2014.29). However, in order to study the kinematics of some of the sources in the region (VLA 4 and Bc) we also reanalyzed previously reported, multi-epoch VLA archive data (project codes AT141, AF381, and AS831; see Carrasco-González et al. 2010, for details on the observations). This, along with our new K-band observations, allow us to compute proper motions in a period spanning 22 years, from 1992 to 2014. Calibration of these archive data was undertaken following standard VLA procedures, using the Astronomical Image Processing System (AIPS 3 ) data reduction package.
Parameters of the data sets and images are summarized in Table  1 and Table A1, respectively.

ALMA
W75N (B) was observed with ALMA at 1.3 mm during three sessions, on May 6th, 7th, and 11th 2018 (archive ALMA data, project code: 2017.1.01593.S). In total, approximately 16 minutes were spent on source. During the session on May 7th, only one minute of useful data on W75N (B) was obtained. The phase center for the W75N (B) observations was RA(J2000) = 20:38:37.0 and Dec(J2000) = +42:37:51.0, which is ∼18 arcsec north of the VLA 1-VLA 2-VLA 3 sources. As a result, the mm continuum sources discussed in this work are detected towards the edge of the ALMA primary beam (FWHM 27 arcsec). The ALMA images presented in this paper (Section 3.2) have not been corrected by primary beam but they are of good enough quality for the identification of different mm continuum sources. The observations were performed using four spectral windows (spws). Two spws had 1.875 GHz bandwidth and were centered on 217.117 and 230.552 GHz. Two further spws had 117.188 MHz bandwidth and were centered at 216.124 and 231.334 GHz. All spws had 1920 spectral channels after hanning smoothing. During the observations, 46 ALMA telescopes participated, with a minimum baseline length of 15 m and a maximum baseline length of 500 m. This resulted in a maximum recoverable scale of ∼7 arcsec.
The observations were initially calibrated using the ALMA pipeline of CASA 5.1.1 (McMullin et al. 2007). Subsequently, after excluding the channels showing line emissions, we performed two interactions of phase self-calibration on the continuum emission by using CASA 5.1.1. This improved the continuum signal-to-noise ratio by a factor of three. Next, the data were imaged and cleaned (task TCLEAN) centered on VLA 2, using Briggs weighting with a robust parameter of 0.5, yielding a synthesized beam size for the 1.3 mm continuum observations of 1. 73 × 0. 86 with a position angle of −4 • (Fig. 3).  -4, 9, 13, 18, 25, 50, 100, and 200 times the rms, 8 µJy/beam. Panels (b) and (c) show a close-up of the northern region containing sources VLA 1, 2, and 3, and the southern region containing sources Bc and VLA 4, respectively. In both cases, intensity contours of panel (a) are shown over the spectral index map (color scale). The pixels shown in spectral index maps are those with S/N > 7 in the continuum image. Synthesized beam = 0. 12×0. 09 (PA = -69 • ).

VLA
The most recent data set (epoch 2014.29, see Table A1) provides radio continuum images of W75N (B) with unprecedented sensitivity. Also, these observations enable us to perform a detailed study of different structures associated with the sources and their emission nature within several ranges of frequencies.
In Fig. 2 we show the images of the radio continuum sources at C, Ku, K, and Q bands, covering the entire bandwidth in each case (see table A1 for the image parameters). Due to spatial filtering effects that worsen at increasing frequency for a given configuration, C and Ku bands are more sensitive to extended emission than K and Q bands, which tend to favor the detection of compact structures. Taking this into account, and given the wide frequency coverage of our observations, both extended and compact emission of the sources can be studied with high sensitivity. We note that all the sources are detected with a S/N ratio much higher than 5 from C to K band (Figs. 2a,b,c). However, at the highest frequency (Q band in  -4, 5, 7, 9, 15, 30, 100, 200, 300, 600 times 10 µJy/beam (22 GHz image, natural weighting); -3, 5, 10, 30, 50, 100, 300, 500 times 20 µJy/beam (44 GHz image, natural weighting). In each panel, the synthesized beam is indicated with a white ellipse at the bottom left.   Table  2 (a detailed discussion of the parameters of VLA 2 is given in Carrasco-González et al. 2015).
In addition to all these sources, three new weak (<100 µJy, see Table 2) compact radio continuum sources are detected in the images at Ku and K bands, as well as in the image obtained by combining all four bands. These new sources are located ∼8 arcsec northeast (VLA [NE]) and ∼6 arcsec southwest (VLA [SW]) from VLA 2, and ∼0.5 arscec northeast (Bd) from VLA 4. Figs. 3 and A1 show the position and contour maps, respectively, of these weak radio sources).

ALMA
The ALMA continuum observations at 1.3 mm show four cores in a region of ∼14 arcsec (MM1, MM2, MM3, MM[N]; Fig. 3). Three of them (MM1, MM2, MM3) have been previously identified by Minh et al. (2010) with the Submillimeter Array (SMA) at 217 and 347 GHz, with angular resolution similar to that in our ALMA images. The fourth millimeter core, MM[N], is located ∼9 arcsec north from MM1 and has not been previously reported. We did not detect with the VLA any radio continuum source toward MM[N], with an upper limit of ∼30 µJy (4σ in the combination of C+Ku+K+Q bands; Table A1). From Fig. 3 we see that the massive protostars VLA 1, VLA 2, and VLA 3 are associated with the brightest millimeter core, MM1, although the limited angular resolution of the ALMA observations (1. 73 × 0. 86, PA = −4 • ) and the north-south distribution of the sources (with a maximum angular separation  cating they are embedded YSOs. However, the fact that Bc, VLA 4, and Bd are not clearly associated with any millimeter continuum core suggests that they are probably not YSOs (see Section 4.3 for the discussion, in particular, on the nature of the Bc and VLA 4 radio continuum sources).
Because these millimeter continuum cores are found at the edge of the ALMA primary beam and they were only observed at a single frequency band, we are not able to derive with accuracy their physical parameters with the present data. For those estimates we refer to Minh et al. (2010).
In Fig. 3 we also show the images of the integrated intensity and velocity field (first order momentum) of the molecular lines CH 3 OH [5(1,4)-4(2,2); rest frequency 216.94552 GHz], SO 2 [22(2,20)-22(1,21); 216.64330 GHz], and SiO [v=0 (5-4); 217.10498 GHz] as observed with the ALMA archive data. A main molecular core centered on VLA 1-VLA 2-VLA 3 is detected through the CH 3 OH and SO 2 lines. This molecular core, of ∼4 arcsec (∼5200 au) size, exhibits a velocity gradient of ∼5 km/s along the northwest-southeast direction, which is fully consistent with the velocity gradient reported by Minh et al. (2010) in H 2 CO with the SMA (∼1 km/s arcsec −1 ). Given their angular resolution, these ALMA observations (and the SMA observations; Minh et al. 2010) cannot resolve the structure and motions of the dust and molecular gas around each of the individual sources VLA 1, VLA2, and VLA 3.
On the other hand, the SiO emission shows an irregular distribution covering a broader velocity range (V LSR ≈ 10-35 km/s) than CH 3 OH and SO 2 (Fig. 3), with the highest velocity emission (V LSR ≈ 30-35 km/s) located ∼6 arcsec northeast from VLA 1-VLA 2-VLA 3. This high-velocity SiO emission could be tracing outflow motions driven by any of the central VLA sources. ALMA observations with higher angular resolution are clearly necessary to identify, isolate, and study the expected different dust and molecular gas components around the individual massive protostars.

VLA 1
It can be seen from Figs. 1 and 4 that VLA 1 exhibits a tail-shaped extended component. This is the first time that this particular structure is observed in VLA 1, due to the high sensitivity of our images. In Fig. 1b we show the spatial distribution of spectral indices (α, defined as S ν ∝ ν α ) along the source, covering the whole range of frequencies. It can be noticed that the central region of VLA 1 presents positive spectral indices (∼ +0.5, measured from the spectral index map, at the continuum emission peak), while the tailshaped structure is dominated by a flat spectrum with α 0. By studying the emission at different bands, we can see that the morphology of the VLA 1 source seems to be composed by an extended and a compact components (see Figure 4). At the lowest frequencies (15 and 22 GHz), the higher optical depth and lower angular resolutions emphasize the extended component. In Figures  4a and 4b, we can clearly see a "cometary" tail curved in the northeast direction. As we go to higher frequencies (44 GHz), the higher angular resolution and lower optical depth allow us to filter out most of the low-brightness extended emission, and the most compact higher brightness emission is resolved into an elongated source in the NE-SW direction (PA = +42 • ± 5 • ; see Figures 4c and 4d).
To study the emission nature of both, the extended and the compact components of VLA 1, we compute the spectral energy distribution (SED) of the source in the whole range of observed frequencies with low angular resolution, and the SED at Q band (where most of the extended emission is filtered) with high angular resolution. In Fig. 5 (top panel) we show the SED over the entire range of frequencies, obtained by measuring flux densities in images with 2 GHz bandwidth, using uniform weighting for C band and natural weighting for Ku, K, and Q bands. All images were convolved to 0. 37, corresponding to the lowest resolution in C band. Flux densities were determined by a Gaussian fit within a circular region of 1.25 arcsec diameter enclosing the source. We note that these data were observed with the telescope array same configuration, and in this case, the largest scales in the images could be more heavily filtered at high frequencies, which could result in spuriously lower values of the spectral index. However, we limited the study of the low angular resolution SED to the core of the emission, which has a size of ∼400 milliarcsec (mas). Emission with this size is fully recovered at all bands. Just for description purposes, we have performed a fit to the observed flux densities from 4 to 47 GHz through an ad hoc function S ν = aν α [1−e −b/ν β ] (Fig. 5, top panel). The fit gives a = 0.62, α = 1.38, b = 13.29, and β = 1.42, with S ν in mJy and ν in GHz. Within the uncertainties in the observations, this SED is consistent with an HII region thermal bremsstrahlung spectrum, opaque at low frequencies ( 10 GHz) and optically thin at high frequencies ( 20 GHz). The size of the extended emission, including the tail, is of the order of 1 arcsec, corresponding to an extension of ∼0.006 pc. This size is significantly smaller than 0.1 pc, suggesting it could be classified as a Hypercompact (HC) HII region, according to Kurtz (2005).
We want to note, however, that we cannot rule out the possibility that some dust contribution from the extended emission in VLA 1 is present at Q-band (Figs. 4a and 4b). If this were true, the real spectral energy distribution of the ionized gas would be flatter at high frequencies than that shown in Fig. 5 (top panel). In any case, as stated in Section 3.2, very high angular resolution observations with ALMA are clearly needed to individually resolve the dust content of each YSO in the region.
In order to study VLA 1 with as less contribution as possible from the extended emission, we measured flux densities with very high angular resolution in eight R0-weighted (Briggs weighting using parameter robust = 0) images of 1 GHz bandwidth within the Q band (see table A1 for image details). All flux densities are measured within a circular region of 0.16 arcsec diameter enclosing the source. The resulting SED is shown at the bottom panel of Fig. 5. Applying a linear fit to this SED, we obtain a spectral index α = +0.5 ± 0.4, consistent with partially optically thick free-free emission from a thermal radio jet, as predicted by models given by Reynolds (1986), and consistent with typical values measured for thermal radio jets (e.g., Anglada, Rodríguez & Carrasco-González 2018). This is also in agreement with previous works by Baart et al. (1986) and Torrelles et al. (2003), who detected the presence of a radio jet traced by the distribution of OH and H 2 O masers, respectively, along the same direction as the radio continuum emission (PA ≈ +43 • , Torrelles et al. 1997Torrelles et al. , 2003. According to Reynolds models, the distance to the driving source where the jet becomes optically thin varies with frequency as a power law. This distance is interpreted as the angular size of the semi-major axis of the jet, θ ν . Therefore, θ ν ∝ ν −0.7/ , where the index is related to the spectral index α as α = 1.3 − 0.7/ in the case of an isothermal jet with constant velocity and ionization fraction. Thus, from the SED in Q band we derive = +0.9 ± 0.5, suggesting a conical jet ( = 1).
All these results indicate that VLA 1 might consist of a thermal radio jet surrounded by an HCHII region. Assuming this is the case, we can estimate the mass-loss rate and some parameters of the HII region. In order to estimate the protostar mass-loss rate M, we follow Equation (3) where S ν is the flux density at the frequency ν, α the spectral index, θ o = 0.88 rad the jet injection opening angle, estimated as θ o = 2 arctan(θ min /θ maj ), where θ min and θ maj are the deconvolved minor and major axes of the Gaussian fit to the source respectively. T e = 10 4 K is the electron temperature, x o the ionization fraction, V jet the jet velocity, i the jet inclination angle, ν t the turn-over frequency, and d = 1.3 kpc the distance to the region W75N (B). The values of ν and S ν correspond to the Q-band image, being ν = 44 GHz the central frequency of the band and S ν = 15.7±0.1 mJy. The spectral index α = +0.5±0.4 is derived from a linear fit to the SED (see Fig. 5, bottom panel), from which we can also infer that the jet emission is partially optically thick up to ν= 47 GHz. Thus, we take this value as a lower limit to the turn-over frequency ν t (above which the entire jet becomes optically thin). Since we cannot specify the jet inclination angle i, we adopt i =45 • , as variations from 45 • to 90 • only change the mass-loss estimate by less than 10%. Moreover, both the ionization fraction and the jet velocity are unknown. Typical values for V jet range from 100 to 1000 km s −1 , while the ionization fraction is usually assumed to be 10% (e.g., Anglada, Rodríguez & Carrasco-González 2018) for low-mass protostars. However, this value is very uncertain, and could probably be higher for high mass protostars. According to this, we estimate lower and upper limits for the mass loss rate of ∼3.5×10 −7 M yr −1 (assuming x o = 1 and V jet = 100 km s −1 ) and ∼3.5×10 −5 M yr −1 (assuming x o = 0.1 and V jet = 1000 km s −1 ), respectively. On the other hand, we can obtain some estimates of the physical parameters of the HII region using data from C to K bands. Derived parameters are presented in Table 3. Within the Rayleigh-Jeans regime, the brightness temperature T B can be written in terms of the flux density S ν at frequency ν, and the solid angle subtended by the source Ω S , as T B = S ν c 2 2kν 2 Ω S . The solid angle Ω S of the elliptical Gaussian fitted to the source brightness profile is calculated as (π/4 ln 2) × FWHM maj × FWHM min . At the central frequency of each band, T B is always small compared with the electronic temperature T e , assumed to be of the order of 10 4 K for an HII region (since this is the temperature at which hydrogen ionizes). Knowing the brightness temperature, the optical depth τ ν can be calculated from T B = T e (1 − e −τ ν ). As the HII region is more optically thin at the K band, we choose this band to estimate some parameters. Assuming the ionization number equals the recombination number, we can estimate the ionizing photon rate N, i.e. the number of ionizing photons λ < 912 Å per unit of time, necessary to account for the emission observed at K band: where η = 3 × 10 −13 cm 3 s −1 the "case B" recombination coefficient (i.e., the number of recombinations per unit time, volume, and electron and ion density) to levels 2 for T e ≈10 4 K; n e and n p are, respectively, the number density of electrons and protons (assumed to be the same), and L is the characteristic size of the region (estimated as the geometric mean of the source major and minor axes). Assuming a homogeneous, spherical HII region of depth L, n e can be expressed in terms of the emission measure (EM) and the geometrical depth as n e = (EM/L) 1/2 (with EM = ∫ L 0 n e n p dl which, in the case of a homogeneous ionized hydrogen medium of depth L approximates to E M n 2 e L). In turn, EM can be derived from the expression of the optical depth as EM[cm −6 pc] = 12.2τ ν [T e /K] 1.35 [ν/GHz] 2.1 . Thus, computing S ν , Ω, T B , τ ν , L, n e , and EM, at the central frequency of band K (Table 3), we finally obtain an ionizing photon rate N ≈ 6×10 44 photons s −1 . This value is much lower than typical estimations for O-B stars (≈10 48 photons s −1 ). Moreover, characteristic values of emission measure and electron number density for HCHII regions are EM 10 10 pc cm −6 and n e 10 6 cm −3 , respectively (Kurtz 2005). In the case of VLA 1, the electron number density we obtain is of the same order of typical values for HCHII regions, while the emission measure is about two orders of magnitude lower (see Table  3).
Mass-loss rates in the range of 10 −10 M yr −1 (for lowmass YSOs) to 10 −5 M yr −1 (for high-mass YSOs) have been determined in the literature (see Anglada, Rodríguez & Carrasco-González 2018, and references therein). Thus, our estimates for the mass-loss rate of VLA 1, together with the physical parameters we obtain for the HII region, lead us to conclude that VLA 1 might be a massive protostar driving a thermal radio jet, which seems to be at the very beginning of the photoionization stage. The presence of a radio jet coexisting with an UCHII has been also re- the SED is computed in the whole range of observed frequencies. Flux densities are obtained from Gaussian fits to 2 GHz bandwidth images within a circular region of 1.25 arcsec diameter enclosing the source. We use uniform weighting at C band and natural weighting at Ku, K, and Q bands. A fit to the measured flux densities is also shown (S ν = aν α [1−e −b/ν β ], with a = 0.62, α = 1.38, b = 13.29, and β = 1.42, with S ν in mJy and ν in GHz). This fitted ad hoc function is only for description purposes of the observed SED (see Section 4.1). Bottom panel: Spectrum at Q band. Flux densities are obtained from Gaussian fits to 1 GHz bandwidth images within a circular region of 0.16 arcsec diameter. We use Briggs weighting (robust 0). The solid line is a linear least-squares fit to the log data, from which we derive a spectral index α = +0.5 ± 0.4 that correspond to = +0.9 ± 0.5, consistent with a thermal radio jet (see Section 4.1). All data points are shown with measurement errors, considering both fitting and calibration uncertainties. Note that these panels trace different components in the source: while the top panel corresponds to the HCHII region, the bottom one traces the compact jet.

VLA 3
VLA 3 was previously proposed to be a partially optically thick compact HII region (e.g., Torrelles et al. 1997;Shepherd, Kurtz & Testi 2004), but it was later suggested to be a thermal radio jet with a spectral index α 3.6−2cm = +0.6 ± 0.1 (Carrasco-González Figure 6. Dependence of the angular size of the jet with frequency (top panel) and spectral energy distribution (bottom panel) of VLA 3. Flux densities and semi-major axes θ ν , are obtained from Gaussian fits to the brightness profile of VLA 3, in 2 GHz bandwidth images. To measure flux densities we use uniform-weighted images (C band) and natural-weighted images (Ku, K, and Q bands), while angular sizes were measured in uniformweighted images (C, Ku, and K bands) and  images (Q band). Solid lines are linear least-squares fits to the log data. Angular size error bars correspond to fitting errors, while flux error bars consider both fitting and calibration errors. The flux density S, r min , r maj , and the characteristic size L correspond to values derived from a by-dimensional Gaussian fit to VLA 1 at K-band (robust 0 weighting). et al. 2010). In our data ( Figs. 1 and 2), VLA 3 appears as an elongated source, with its major axis oriented in the northwestsoutheast direction at all wavelengths, with a position angle PA = -17 • ±2 • (Fig. 1). The ionized emission is characterized by a spectral index α +1 in the central region of VLA 3 (Fig. 1b), consistent with partially optically thick free-free emission. Its elongated shape and its spectral index suggest that VLA 3 could be a thermal radio jet. Thus, according to theoretical models given by Reynolds (1986), we expect to find that both the SED and the jet angular size depend on the observed frequency as power laws (see Section 4.1). Therefore, we measured the flux densities and angular sizes in several images of 2 GHz bandwidth each, covering a frequency range from 7 to 48 GHz. Flux densities are measured within a circular region of 0.5 arcsec diameter in each image, and the angular sizes correspond to the deconvolved major axis of the bi-dimensional Gaussian fit to the emission. Sizes vary in the range 30-200 mas, and therefore, we can rule out that a significant amount of extended emission is filtered out at the highest frequencies. In Fig. 6 we can see that both the angular-size θ ν (top panel) and the SED (bottom panel) do vary as power laws of the frequency. In the ideal case of a conical thermal jet, with constant velocity, temperature, and ionization fraction, values of +0.6 and −0.7 are expected for the spectral index and the slope of the size vs frequency, respectively (Reynolds 1986). In our case, the behavior of the size with frequency is consistent with a conical jet, while the slightly larger value of the spectral index would indicate some deviation from the ideal physical conditions. These results support that VLA 3 is associated with a thermal radio jet as previously proposed by Carrasco-González et al. (2010). Therefore, in order to estimate the protostar mass-loss rate M we follow Equation 1. In this case ν, S ν , and α correspond to the combined image (C+Ku+K+Q bands), being ν = 26 GHz the central frequency of the band, S ν = 9.06±0.07 mJy, and α = +1.27 (computed from the spectral index map at the continuum emission peak). From the SED (bottom panel of Fig. 6) we can see that the emission is partially optically thick up to ν= 47 GHz, thus, we take this value as a lower limit to the turn-over frequency ν t . As in the case of VLA 1, we also adopt a jet inclination angle i =45 • , and estimate lower and upper limits for the mass loss rate considering different approximations for the ionization fraction and the jet velocity, i.e., ∼4×10 −6 M yr −1 (assuming x o = 1 and V jet = 100 km s −1 ) and ∼4×10 −4 M yr −1 (assuming x o = 0.1 and V jet = 1000 km s −1 ), respectively. Such mass-loss rates are significantly higher than those estimated in low-and intermediate-mass YSOs (e.g., Beltrán et al. 2001;Anglada, Rodríguez & Carrasco-González 2018), but similar to the values obtained in high-mass YSOs (e.g., Rodríguez et al. 1994;Guzmán et al. 2012;Añez-López et al. 2020), supporting that VLA 3 is excited by a massive protostar.

Bc and VLA 4
In Fig. 1c we show the radio image with the highest resolution and sensitivity to date of Bc and VLA 4. This allows us to resolve their structure, and study the nature of their emission through the spectral index map. The source Bc is clearly resolved into two components (labeled as Bc [E] and Bc [W] in Fig. 1). We note that Bc [W]-Bc [E] form an elongated structure, with its minor axis aligned with the VLA 3 jet direction as we would expect to observe in a bow-shock produced by the impact of a supersonic jet with the environment gas (e.g., Tafalla et al. 2017;Castellanos-Ramírez, Raga & Rodríguez-González 2018). This supports the scenario proposed by Carrasco-González et al. (2010) who interpreted Bc as an obscured radio  (1992.98, 2001.40, 2006.47, 2014.29): RA(J2000) = 20h 38m 36.47s, DEC(J2000)= 42 • 37 34.15 . The solid lines are least-square fits to the data. Velocities on the plane of sky are derived by assuming a distance to the region of 1.3 kpc (Rygl et al. 2012). Bc and VLA 4 are moving away from the system with PAs of ∼-20 • and ∼-10 • , respectively.
Herbig-Haro (HH) object, possibly excited by the VLA 3 jet. A flattened structure similar to that of Bc is also seen in the frontal region of the shock of the obscured HH 80N object (Rodríguez-Kamenetzky et al. 2019). Carrasco-González et al. (2010) studied the kinematics of these sources by computing proper motions relative to VLA 3 in three epochs (1992.90, 1998.23, and 2006.38) spanning 13.48 years. Adopting a distance to the region of 2 kpc (Dickel, Wendker & Bieritz 1969), they derived for Bc a velocity of 220±70 km s −1 moving on the plane of the sky, and toward the south, approximately along the major axis of the VLA 3 radio jet. Regarding VLA 4, they suggested it could be either an independent star or shock-excited gas produced by a previous ejection from VLA 3. However, even though they noticed a small displacement of VLA 4 to the south with respect to VLA 3 between 2000 and 2006, they were not able to distinguish between these two scenarios. Regarding this, we measured the proper motion of both sources along a time span of 22 years, from 1992 to 2014. Positions of the sources in 2014 were measured in the K band image, since it is the highest angular resolution image where Bc and VLA 4 are detected (see Fig. 2). We compute the proper motions relative to the average position of the system VLA 1-VLA 2-VLA 3 assuming it is stationary, instead of VLA 3 only. In this way, we reduce possible errors associated to variations in the shape of VLA 3 due to observational effects (e.g. differences in the beam size at each epoch) or real changes (e.g., new ejections from VLA 3). Although Bc shows substructure, we measure the displacement of the two-component complex. Fig. 7 shows that both Bc and VLA 4 are moving away from the system in the southeast direction (PA -20 • and -10 • , individually, see Fig.  7) with velocities on the plane of the sky of 1.9×10 −2 arcsec yr −1 (∼118 km s −1 ), and 1.8×10 −2 arcsec yr −1 (∼112 km s −1 ), respectively (assuming the updated distance to the region of 1.3 kpc; Rygl et al. 2012). The velocity of Bc differs from the value given by Carrasco-González et al. (2010), but this is mostly due to the difference in the adopted distance to the region (considering a distance of 2 kpc we obtain velocities of 170 and 180 km s −1 for VLA 4 and Bc, respectively, closer to the results reported by these authors). Thus, the shape and proper motions of Bc and VLA 4 are consistent with both sources tracing shock-excited gas (obscured HH objects). Moreover, from the spectral index map (Fig. 1c) we see that Bc and VLA 4 are dominated by flat spectral indices (α ∼ 0), as it would be expected for optically thin free-free emission produced by shock-ionized material. Thus, flat spectral indices, together with the shape of the sources, their proper motions, and propagation direction constitute solid evidence supporting the shock scenario.
Among the sources VLA 1, VLA 2, and VLA 3, VLA 1 and VLA 2 have associated outflows along the northeast-southwest direction (e.g., Torrelles et al. 1997Torrelles et al. , 2003. Therefore, VLA 3 is the only source elongated in the northwest-southeast direction (PA -17 • ), consistent with the direction of the proper motions of Bc and VLA 4 (PA -20 • and -10 • ). This, along with the results found in section 4.2, suggests that VLA 3 is the driving source of Bc and VLA 4. In addition, the fact that Bc and VLA 4 are not associated with any of the detected millimeter cores (see Section 3.2; Fig. 3) suggests that they are not protostars, further supporting our shock-excited gas interpretation for these sources. This kind of obscured HH objects, exhibiting proper motions higher than 100 km s −1 , have also been observed in radio continuum in other intermediate-and high-mass star-forming regions: e.g., Serpens (Curiel et al. 1993;Rodríguez-Kamenetzky et al. 2016), GGD 27 (Martí et al. 1995(Martí et al. , 1998Masqué et al. 2015;Rodríguez-Kamenetzky et al. 2019), Cepheus A (Curiel et al. 2006).

CONCLUSIONS
We presented an analysis of high-sensitivity, high-resolution multifrequency VLA observations of the massive star-forming region W75N (B), together with complementary studies performed with ALMA and VLA archive data. Our study leads us to the following conclusions: • VLA 1 is detected at all the observed frequencies (4-48 GHz). Its SED over the entire range of frequencies is consistent with thermal free-free emission from an HCHII region ( 1 arcsec, 1300 au), while the high angular resolution spectrum of the most compact component at high frequencies (40-48 GHz) is compatible with a thermal radio jet at scales of ≈0.1 arcsec (≈130 au), with a spectral index α ≈ +0.5 (S ν ∝ ν α ). This suggests that VLA 1 is driving a thermal radio jet, and it is likely at the early stage of the photoionization.
• VLA 3 shows an elongated structure at scales of few tenths of arcsec (few hundred of au), with its major axis oriented in the northwest-southeast direction (PA ≈ -17 • ). Both the SED and the size dependence with frequency indicates that this source is also driving a thermal radio jet.
• We computed proper motions of the radio continuum sources Bc and VLA 4 in a time interval of 22 years. We found both sources are moving away toward the south, in a similar direction as the VLA 3 thermal radio jet, with velocities of ≈112-118 km s −1 (≈1.8-1.9×10 −2 arcsec yr −1 ). From the SED analysis we found these sources are dominated by flat spectral indices, as it is expected for optically thin free-free emission produced by shock-ionized material. These results support the scenario in which Bc and VLA 4 are obscured HH objects tracing shocks of the jet driven by VLA 3.
• Four 1.3 mm continuum cores are observed with ALMA (MM1, MM2,MM3,and MM[N]) in a region of ∼14 arcsec. Three of these millimeter cores, MM1, MM2, and MM3, had previously been identified with the SMA interferometer, while MM[N] had not been previously reported. VLA 1, VLA 2, and VLA 3 are associated with the brightest core MM1. Bc and VLA 4 are not associated with any of the millimeter continuum cores, supporting they are not YSOs but shock-excited gas as concluded from our VLA observations.
• We have detected three new weak compact radio continuum sources (VLA[SW], VLA [NE], and Bd). Two of them, VLA [SW] and VLA [NE] are associated with the millimeter cores MM2 and MM3, respectively, suggesting they are embedded YSOs belonging to the W75N (B) massive star-forming region.
• With our VLA observations we have identified a cluster of at least five YSOs (VLA 1, VLA 2, VLA 3, VLA[SW] and VLA [NE]) in a region of ∼10 arsec (∼13000 au). All of the sources for which information was obtained on their structure and SED, exhibit accretion/outflow activity at different relative stages of their evolution. In this sense, given that VLA 1 has indications that it has already started the photoionization stage, it could be relatively more evolved than VLA 2 and VLA 3. To further characterize this cluster of YSOs, ALMA observations with very high-angular resolution are needed to resolve individually the gas/dust content of each YSO, as well as to study their expected different chemical compositions.  Fig. A1 shows a radio continuum map of the region, with closeups of the three new compact sources of < 100 µJy detected in the field: VLA[NE], VLA [SW], and Bd. This paper has been typeset from a T E X/L A T E X file prepared by the author.