Massive Protostars in a Protocluster—A Multi-scale ALMA View of G35.20-0.74N

Massive stars dominate the radiative, mechanical, and chemical feedback to the interstellar medium, by which they regulate the evolution of galaxies. Their feedback also affects nearby low-mass forming stars and their protoplanetary disks. In spite of their importance, many aspects of massive star formation remain poorly understood (see, e.g., Tan et al. 2014; Motte et al. 2018).

One of the challenges to form a massive star is its strong feedback on its natal gas core. In particular, protostellar extreme-ultraviolet radiation can ionize the accretion disk and infall envelope to produce photoevaporative outflows driven by the thermal pressure of the ∼10⁴ K ionized gas (Hollenbach et al. 1994). Detection of a H ii region via free–free continuum emission and/or hydrogen recombination lines (HRLs) has traditionally been considered as an indication of the end of the accretion phase (e.g., Churchwell 2002). However, theoretical calculations have suggested that such feedback may not be strong enough to stop accretion (e.g., Keto 2002; Peters et al. 2010; Tanaka et al. 2017; Kuiper & Hosokawa 2018), and recent high-resolution observations have revealed cases of ongoing accretion after the onset of photoionization (e.g., Zhang et al. 2019c; Maud et al. 2019; Guzman et al. 2020; Moscadelli et al. 2021). Understanding the accretion status after the onset of photoionization is therefore crucial to constrain the growing efficiency under feedback during the later stages of the birth of massive stars.

Another important question is how accretion proceeds in massive star formation. There are a growing number of examples of accretion disks around massive protostars (e.g., Sánchez-Monge et al. 2013; Johnston et al. 2015; Ilee et al. 2016; Ginsburg et al. 2018; Maud et al. 2018; Zhang et al. 2019b; Tanaka et al. 2020; see also Beltrán & de Wit 2016 and references therein), suggesting that disk-mediated accretion is common during massive star formation. On the other hand, there are also some examples of massive young stars in which the accretion appears to occur without mediation from a large, stable disk, but rather from multiple, more chaotic flows (e.g., Goddi et al. 2020). Disks around massive protostars may be prone to gravitational instability and thus relatively likely to fragment and form binaries or close multiple systems. However, direct observations of substructures in massive protostellar disks are rare (e.g., Motogi et al. 2019; Johnston et al. 2020). Massive protostellar binary systems with a few hundreds of au separations have been speculated to be the result of such disk fragmentation, with continued accretion mediated by a circumbinary accretion disk or pseudo-disk (e.g., Beltrán et al. 2016; Kraus et al. 2017; Zhang et al. 2019a), although the morphology of such a circumbinary accretion structure may have evolved significantly after binary formation (e.g., Tanaka et al. 2020). Therefore, detailed, multi-scale observations of the hierarchical structures of the accretion flows around massive protostars are important to understand the global accretion process involved in massive star formation.

On larger scales (≳0.01 pc), the fragmentation of the molecular clump leads to the formation of a cluster with both low- and high-mass stars (see, e.g., Motte et al. 2018). The degree of fragmentation is found to be diverse among different massive star-forming regions (e.g., Palau et al. 2013; Beuther et al. 2018), which is expected to be affected by various factors such as the initial density distribution, magnetic field strength, degree of turbulence, and radiative heating of the cloud. The distributions of the massive stars within the cluster and the spatial separations of the members can provide important information on the fragmentation mechanism and help distinguish if massive stars form via core accretion (e.g., McKee & Tan 2003) or competitive accretion (e.g., Bonnell et al. 2001; Wang et al. 2010). For example, Bonnell & Davies (1998) concluded massive stars formed preferentially in the center of the Orion Nebular Cluster. On the other hand, Moser et al. (2020) did not find strong evidence for primordial mass segregation in the massive protostar population of the massive Infrared Dark Cloud G028.37+00.07. Gravitational interactions and/or radiative and mechanical feedback from massive forming stars in crowded protocluster environments may also affect neighboring low-mass protostars and their disks. It is therefore important to perform detailed studies of massive star formation in the context of star cluster formation.

G35.20-0.74N (a.k.a. IRAS 18566+0136; hereafter G35.2) is an ideal target to study all these processes to understand the formation of massive stars. It is a well-known massive star-forming region at a distance of 2.2 kpc (Zhang et al. 2009; Wu et al. 2014). Previous observations of the Atacama Large Millimeter/submillimeter Array (ALMA) and the Submillimeter Array (SMA) have revealed a filamentary structure with a string of embedded cores (Qiu et al. 2013; Sánchez-Monge et al. 2013, 2014). CO observations have revealed a wide outflow structure extending to >1' from the central source (Gibb et al. 2003; Birks et al. 2006) in the direction of the northeast–southwest (NE–SW), which is near perpendicular to the filament. Consistent with this wide CO outflow, near-infrared (NIR) H₂ observations also show a wide outflow structure in the NE–SW direction, which may contain two outflows (Caratti o Garatii et al. 2015). On the other hand, centimeter radio observations of the Very Large Array (VLA) have revealed a collimated bipolar ionized jet along the north–south (N–S) direction (Heaton & Little 1988; Gibb et al. 2003; Beltrán et al. 2016), with its driving source being one of the sub-mm cores of the filament (core B, Sánchez-Monge et al. 2013, 2014; Beltrán et al. 2016). Fedriani et al. (2019) revealed NIR emission spatially coincident with the radio jet in the central region. This N–S outflow is also prominent in mid-infrared (MIR; De Buizer 2006). At these wavelengths, the emission is elongated in the N–S direction but peaked to the north of the identified radio source and continuum core and is thought to trace the blueshifted outflow cavity. At longer wavelengths of 30–40 μm, the SOFIA-FORCAST observations have revealed the southern, far-facing outflow cavity (Zhang et al. 2013; De Buizer et al. 2017). The driving source of the N–S outflow (core B) is reported to have a Keplerian accretion disk (Sánchez-Monge et al. 2013), with the dynamical mass estimated to be about 18 M_⊙. Beltrán et al. (2016) further identified a binary system in core B (VLA sources 8a and 8b). By fitting the infrared spectral energy distribution (SED) assuming a single protostar, core B is estimated to be a protostar with a mass of 12–24M_⊙ and a total luminosity of (4–8) × 10⁴ L_⊙ (Zhang et al. 2013; De Buizer et al. 2017; Zhang & Tan 2018). In this paper, we present a detailed study of G35.20-0.74N using multiwavelength continuum, molecular line, and HRL data obtained by ALMA and VLA.

The observations were carried out with ALMA in Band 6 (1.3 mm) on 2016 April 24 with the C36-3 (hereafter C3) configuration, on 2016 September 8 with the C36-6 (hereafter C6) configuration, and on 2017 November 4 with the C43-9 (hereafter C9) configuration (ALMA project ID: 2015.1.01454.S and 2017.1.00181.S). The total integration times are 3.5, 6.5, and 16.8 minutes in the three configurations. Forty-one antennas were used with baselines ranging from 15–463 m in the C3 configuration, 38 antennas were used with baselines ranging from 15 m to 2.5 km in the C6 configuration, and 49 antennas were used with baselines ranging from 113 m to 13.9 km in the C9 configuration. J1751+0939, J1924-2914, and J2000-1748 were used for bandpass calibration and flux calibration. J1830+0619 and J1851+0035 were used as phase calibrators. The science target was observed with single pointings, and the primary beam size (half power beamwidth) was 229.

The data were calibrated and imaged in CASA (McMullin et al. 2007). After pipeline calibration, we performed self-calibration using the continuum data obtained from a 2 GHz-wide spectral window with line-free bandwidth of 1.2 GHz. We first performed two phase-only self-calibration iterations with solution intervals of 30 and 6 s, and then one iteration of amplitude self-calibration with the solution interval equal to the scan intervals(ranging from ∼1 to ∼5 minutes). We applied the self-calibration phase and amplitude solutions to the other spectral windows. Such self-calibration was performed for the data of three configurations separately before any combining was performed. The CASA tclean task was used to image the data, using Briggs weighting with the robust parameter set to 0.5. In order to better show structures on different scales, below we present images made of different configuration combinations. The resolution of the C3+C6 configuration is about 025, and the resolution of the image of just the C9 configuration is about 003, the highest achieved so far for this region at millimeter wavelengths. Detailed synthesized beam sizes of images of different configuration combinations, as well as their rms noise levels, are listed in Table 1.

Table 1. Properties of the Presented Continuum Images

Telescope	Band	Wavelength	Configuration	Synthesized Beam	Noise rms	Noise rms
		(mm)			mJy beam−1	K
ALMA	Band 7 a	0.85	Extended (cycle 0)	047 × 042 (P. A. = 514)	0.73	0.039
ALMA	Band 6	1.3	C3+C6	025 × 025 (P. A. = −363)	0.48	0.17
ALMA	Band 6	1.3	C3+C6+C9	0045 × 0033 (P. A. = −640)	0.085	1.3
ALMA	Band 6	1.3	C6+C9	0043 × 0031 (P. A. = −633)	0.088	1.5
ALMA	Band 6	1.3	C9	0037 × 0027 (P. A. = −660)	0.061	1.4
VLA	Q band b	7	A	0046 × 0036 (P. A. = −520)	0.020	7.7
VLA	Q band b	7	B	014 × 012 (P. A. = −445)	0.030	1.1
VLA	Ku band b	13	A	0079 × 0076 (P. A. = −463)	0.013	5.2
VLA	Ku band b	13	B	025 × 024 (P. A. = 395)	0.020	0.8

Notes.

^a Reimaged using the data published by Sánchez-Monge et al. (2013). ^b Beltrán et al. (2016).

Download table as:  ASCIITypeset image

We also utilize the ALMA Band 7 (0.85 mm) continuum data published by Sánchez-Monge et al. (2013). Please refer to that paper for details of the observation. Instead of directly using the published 0.85 mm continuum image, we reproduced the continuum image using the calibrated data obtained from the ALMA archive. From the four 2 GHz wide spectral windows, we imaged the continuum emission using line-free channels, which span a total of 0.86 GHz. We performed two phase-only self-calibration iterations with solution intervals of 30 and 6 s, and then one iteration of amplitude self-calibration with the solution interval equal to the scan intervals (ranging from ∼1 to ∼5 minutes). The CASA tclean task was used to image the data, using Briggs weighting with the robust parameter set to 0.5.

In addition, we also utilize the VLA K-band (1.3 cm) and Q-band (7 mm) continuum images published by Beltrán et al. (2016). See that paper for details of the VLA observations. The synthesized beams and noise levels of these images are also listed in Table 1.

3.1. Lower-resolution 1.3 mm Continuum

Panels (a) and (b) of Figure 1 show the ALMA 1.3 mm continuum image, combining the data of the C3 and C6 configurations, and the reproduced 0.85 mm continuum image. In general, the 1.3 mm continuum morphology is consistent with that of the 0.85 mm continuum, i.e., showing multiple cores along the filamentary structure. Compared to that presented by Sánchez-Monge et al. (2013), the new 0.85 mm continuum image shows a compact source to the west of core A (labeled as source 3 in the figure; see Section 3.2). This source is seen as a compact continuum source in the 1.3 mm continuum image, as well as all the VLA images, but was not seen in the previously presented 0.85 mm image. Detection of this source in the reproduced 0.85 mm continuum image is probably due to the improved sensitivity brought by the self-calibration. The new 0.85 mm image also shows more diffuse emission in the southern end of the filament, which is also seen in the 1.3 mm continuum image. Panel (c) shows a comparison between the 1.3 and 0.85 mm continuum. The overall morphology and the locations of the bright sources/regions coincide very well in the two images.

Zoom In Zoom Out Reset image size

Figure 1(d) shows the spectral index map between ALMA 1.3 and 0.85 mm continuum. Here, the spectral index is defined as ${\alpha }_{\nu }=\mathrm{log}({I}_{{\nu }_{1}}/{I}_{{\nu }_{2}})/\mathrm{log}({\nu }_{1}/{\nu }_{2})$, where I_ν is the intensity at the frequency ν, ν₁ = 343 GHz (0.85 mm) and ν₂ = 234 GHz (1.3 mm). The spectral index map is obtained with the 0.85 mm data and 1.3 mm C3+C6 data, using multifrequency synthesis (mtmfs) in the CASA tclean task with nterms = 2. In order to derive the spectral index map more accurately, we only used the 1.3 mm data with baselines between 20 and 440 kλ to match the uv range of the 0.85 mm data.

Typical dust continuum emission has spectral indices of α_ν > 2, with α_ν = 2 in the optically thick case and 3 ≲ α_ν ≲ 4 in the optically thin case (assuming a typical dust emissivity spectral index of 1–2). Free–free continuum emission has spectral indices of −0.1 ≲ α_ν < 2, with α_ν = 2 in the optically thick case and α_ν ≈ −0.1 in the optically thin case. Cores A and B have spectral indices below 2, which suggests a contribution of free–free emission, expected for forming massive stars. However, the observed spectral indices are not far from 2. Optically thick dust emission or a combination of free–free and dust emission can also lead to such spectral indices. The regions with spectral indices above 2, which are dominated by dust emission, are mostly around cores A and B along the filament, suggesting dusty envelopes around the small H ii regions created by the forming massive stars. Cores C, D, and F all have spectral indices close to 3, dominated by dust emissions, while core E has a spectral index slightly above 2, indicating a higher free–free fraction or high optical depth. Source 3 has a very flat SED between 1.3 and 0.85 mm, suggesting it is dominated by free–free emission.

3.2. Long-baseline 1.3 mm Continuum

Figure 2 shows the high-resolution view of the 1.3 mm continuum emission. From the long-baseline continuum data (C9 configuration), we are able to identify a total of 22 compact sources in the whole field of view with a primary beam response >0.1 (20 sources within the radius of primary beam response >0.5). The candidates of compact sources are first identified with maximum intensities >5σ_noise (1σ_noise = 0.06 mJy beam⁻¹) and size estimated for regions with intensities >4σ_noise larger than that of one beam size. We further removed the emission peaks caused by the extended filament emissions or the bright patches of the residual pattern, and made sure that the remaining compact sources clearly stand out from the residual pattern fluctuation. Note that the residual pattern close to the central region has a fluctuation level of 1σ_residual = 0.1 mJy beam⁻¹. And the confirmed compact sources all have peaks ≳6σ_residual (≳10σ_noise). Figure 3 shows the continuum image of each source. Table 2 lists the identified compact sources, as well as their associated cores and corresponding VLA compact sources. The positions of these sources are also marked in panels (a) and (b) of Figure 1 to show their association with the larger-scale structures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Identified Compact Continuum Sources in ALMA 1.3 mm Observations

Source	VLA	Associated	α(ICRS)	δ(ICRS)	Ipeak,1.3 mm c , d	S1.3 mm c , e	Mass g	Radius h	Driving Outflow i
	Source a	Core b	hh:mm:ss	° ' ''	(mJy beam−1)	(mJy)	(M⊙)	(au)
1	3	A	18:58:12.953	1:40:37.359	22.99 ± 0.06	42.19 ± 0.15 (42) f	0.70	35	Outflow 2
2	8a	B	18:58:13.037	1:40:35.931	27.80 ± 0.06	35.20 ± 0.15 (25) f	0.42	17	Outflow 1
3	1		18:58:12.815	1:40:36.604	13.03 ± 0.06	14.52 ± 0.16 (9.7) f	0.16	<12
4	8b	B	18:58:13.015	1:40:36.110	1.37 ± 0.06	1.49 ± 0.15 (1.3) f	0.021	<36
5			18:58:12.991	1:40:36.200	0.74 ± 0.06	1.05 ± 0.15	0.067	<39
6			18:58:12.992	1:40:36.338	0.75 ± 0.06	0.95 ± 0.15	0.060	17
7			18:58:13.070	1:40:35.768	1.07 ± 0.06	1.29 ± 0.15	0.083	15
8			18:58:13.111	1:40:35.831	0.65 ± 0.06	0.65 ± 0.15	0.040	<39
9			18:58:13.070	1:40:37.597	1.89 ± 0.06	3.65 ± 0.15	0.23	34
10			18:58:12.921	1:40:39.988	0.86 ± 0.06	1.57 ± 0.16	0.10	33
11		D	18:58:12.815	1:40:40.000	0.89 ± 0.07	3.77 ± 0.17	0.24	257
12			18:58:12.836	1:40:34.715	1.21 ± 0.06	2.42 ± 0.16	0.15	37
13			18:58:12.725	1:40:35.606	0.58 ± 0.07	0.74 ± 0.17	0.047	11
14		C	18:58:13.131	1:40:33.233	2.98 ± 0.06	6.60 ± 0.16	0.42	44	Outflow 3
15			18:58:13.301	1:40:34.313	0.74 ± 0.07	1.88 ± 0.17	0.12	52
16			18:58:12.957	1:40:30.656	2.12 ± 0.07	4.57 ± 0.17	0.29	38
17		E	18:58:13.190	1:40:30.678	0.99 ± 0.07	0.99 ± 0.18	0.043	<30	Outflow 4
18			18:58:12.883	1:40:31.778	1.00 ± 0.07	1.48 ± 0.17	0.094	27
19			18:58:12.858	1:40:26.403	0.77 ± 0.09	1.51 ± 0.23	0.096	34
20			18:58:12.290	1:40:36.928	0.80 ± 0.10	1.55 ± 0.25	0.099	32
21			18:58:13.955	1:40:27.608	1.64 ± 0.20	3.14 ± 0.48	0.20	44
22			18:58:12.002	1:40:46.487	3.05 ± 0.30	7.83 ± 0.74	0.50	57

Notes.

^a Compact sources identified from VLA observations by Beltrán et al. (2016). ^b Core structures identified from ALMA 0.85 mm observations by Sánchez-Monge et al. (2014). A source is considered to be associated with one of the identified cores only if it is close to the continuum peak of that core. ^c Primary beam response has been corrected for the peak intensities, integrated flux densities and their uncertainties. ^d Peak intensity measured from the C9 configuration image. ^e Flux density integrated over a circle with a radius of 005. The flux density uncertainties are estimated by calculating flux densities with apertures of the same size (005) over random off-source positions many times on the primary beam uncorrected map, and then scaled by the primary beam response at that source position. ^f Numbers in the bracket show the dust continuum flux densities, after subtraction of the free–free emission. The free–free flux densities are estimated from the model fitting to the continuum SED and H30α flux. ^g Gas mass calculated from dust continuum flux densities. A gas-to-dust mass ratio of 100 is assumed. For sources 1−4 (massive sources with photoionized regions), a dust temperature of 100 K is assumed, and for sources 5−22, a dust temperature of 30 K is assumed. ^h Defined as $\sqrt{{D}_{\mathrm{maj}}{D}_{\min }}/2$, where D_maj and D_min are the major axes FWHM and minor axis FWHM of the Gaussian component of the continuum image after deconvolution of the synthesized beam. The 2D Gaussian fit is performed using the CASA task imfit. Upper limits are given if the source is unresolved. ⁱ Identified from ¹²CO maps (see Section 6.2).

Download table as:  ASCIITypeset image

Most of the compact continuum sources are along the filament, especially around cores A and B. Among the six cores previously identified from the lower-resolution 0.85 mm continuum image, five (except core F) have associated compact continuum sources. ALMA source 3, which is a bright and compact source in both the long-baseline and lower-resolution 1.3 mm images is now also shown in our new 0.85 mm continuum image (Figure 1(b)). ALMA source 16 is also seen to be associated with a dense core structure in the new 0.85 mm continuum, which was not reported previously.

Figure 4 overlays the VLA 1.3 and 7 mm continuum emission images in the B and A configurations with the 1.3 mm low- and high-resolution continuum images. As previously pointed out, the large-scale VLA continuum emission mostly traces a radio jet originating from core B (source 2 in our list). As can be seen in this figure, ALMA sources 1−4 have associated emissions in the VLA bands. ALMA sources 2 and 4, which are associated with VLA sources 8a and 8b, respectively, are likely members of a binary system, as suggested by Beltrán et al. (2016).

Zoom In Zoom Out Reset image size

ALMA sources 1−4 are the only four 1.3 mm compact sources associated with VLA compact emission (VLA sources 3, 8a, 1, and 8b). The other emission knots identified in the VLA continuum maps do not have any counterpart compact source at 1.3 mm. These sources mostly have flat or negative spectral indices in VLA bands (see Beltrán et al. 2016). Extrapolating from the VLA bands to ALMA 1.3 mm, their fluxes should be below the continuum detection limits of ALMA unless there are additional dust contributions. Therefore, the non-detection of these VLA sources in ALMA 1.3 mm supports the idea that they are jet knots caused by shock ionization and are not protostars (Beltrán et al. 2016). Such a scenario is also supported by the fact that Fedriani et al. (2019) detected atomic ([Fe ii] and Brγ) and molecular (H₂) emissions associated with shocks toward the north side of core B, i.e., the blueshifted outflow (also see Section 6.2). On the other hand, the other ALMA compact sources that are not detected in VLA bands (ALMA sources 5−22) should be lower-mass protostars (or relatively massive but young sources which have not yet reached the main sequence to produce detectable photoionized regions) forming along with the photoionizing massive stars.

For the ALMA sources with corresponding centimeter emission (sources 1−4), Figure 5 shows their continuum SEDs from 6 cm to 1.3 mm. The VLA flux densities are obtained from Beltrán et al. (2016). The 1.3 mm flux densities are obtained from the long-baseline (C9) data. These sources have spectral indices from 6 cm to 7 mm (VLA bands) either near flat (source 4) or slightly higher than 1 (sources 1−3), consistent with free–free emission. Free–free emission tends to have its spectral index decrease as the frequency becomes higher. However, in all these sources, the spectral indices tend to increase at 1.3 mm (especially in sources 1 and 4), indicating dust contributions. As the free–free emission in these sources is highly compact, the different baseline ranges in these data have minimum effects on the SED slope. For the 1.3 mm data, we also only focus on the most compact emission structures here, so filtering out of the extended dust emission in the C9 configuration data does not affect our discussion. The compact dust emission confirms that they are forming stars rather than jet knots. For sources 1−3, the existence of bright free–free and dust emission is consistent with forming massive stars creating small ionized regions via photoionization. This was speculated by Beltrán et al. (2016) and is now further supported by the ALMA observations. For source 4, while such a scenario is also possible, it is 1.3 mm flux is on a similar level as the other low-mass sources in the region. If source 4 is not yet massive enough to produce enough ionizing photons, the free–free emission could be caused by shock ionization of an unresolved jet, which is seen in protostellar sources with a wide range of masses (e.g., Anglada et al. 2018; Tychoniec et al. 2018). We will discuss these SEDs in detail in Section 5.

Zoom In Zoom Out Reset image size

H30α HRL emission is detected toward sources 2 and 3. Figure 6 shows the integrated H30α emission maps of these two sources and their H30α spectra at the continuum peak positions. The H30α peaks coincide well with the continuum peaks, which supports the existence of small ionized regions at the centers of sources 2 and 3 that contribute both to the H30α emission and the free–free continuum emission. Even though the H30α emission in source 2 is only detected with low signal-to-noise ratio (S/N) (6σ), the fact that the broad spectral feature is spatially coincident with the continuum emission suggests that the H30α detection is likely real. The H30α lines are relatively well fitted with Gaussian profiles. The H30α central velocities are derived to be 29.4 ± 2.2 km s⁻¹ and 30.2 ± 0.2 km s⁻¹ for sources 2 and 3, respectively, which are similar to each other within the uncertainties, and also consistent with the systemic velocity of source 2 (30.0 ± 0.3 km s⁻¹) estimated by fitting the Keplerian disk kinematics traced in different molecular lines (Sánchez-Monge et al. 2013). The line widths are derived to be 86.1 ± 5.2 km s⁻¹ and 45.0 ± 0.6 km s⁻¹ for sources 2 and 3, respectively. The line width of 86 km s⁻¹ is on the upper edge of observed recombination line widths in hypercompact H ii (HCH ii) regions (e.g., Hoare et al. 2007). Although, larger line widths have been reported in some HCH ii regions (e.g., NGC 7538 IRS 1; Sewilo et al. 2004); a line width of ≲50 km s⁻¹ is more typical for such regions.

Zoom In Zoom Out Reset image size

Panels (e) and (f) of Figure 6 show the velocity structures of the H30α lines in sources 2 and 3. The integrated maps of the blue- and redshifted emission show a shift of emission with velocity in these two sources. For source 3, the high S/N further allows us to determine the shift of the emission centroid with channel velocity (Figure 7). The centroid positions are determined by fitting Gaussian ellipses to the H30α emission in the C3+C6+C9 configurations at channels with peak intensities >10σ. Following the method used by Zhang et al. (2019c), the uncertainties of the centroid positions are determined from the data S/N as well as the phase noise in the passband calibrator. The centroid distributions and kinematics are consistent with a rotating structure with a radius of about 5 mas (11 au). The elliptical distribution of the centroids, and the fact that the centroids with the highest velocities are closer to the center, are similar to what has been seen in molecular lines around massive protostars (e.g., Sánchez-Monge et al. 2013; Ilee et al. 2016), which indicates faster rotation velocities near the central source, as expected for a Keplerian disk.

Zoom In Zoom Out Reset image size

Following the methods used by Sánchez-Monge et al. (2013) and Ilee et al. (2016), we fitted the H30α centroid distribution with a model of a Keplerian disk. The model has six free parameters: mass of the central object, the R.A. and decl. offsets of the center, disk position angle, disk inclination, and systemic velocity. The fitting minimizes the differences between the channel velocity and the line-of-sight velocity of disk rotation at the peak position of each channel. However, in our case, different combinations of disk inclination and central mass can generate fitting with similarly good qualities. Therefore, we cannot well constrain these two properties solely by such fittings. To break the degeneracy, we set the central mass to 16 M_⊙, which is estimated from the free–free continuum (see Section 5.2). Figure 7(b) shows the best-fit model. This model has a systemic velocity of V_sys = + 30 km s⁻¹, R.A. offset of Δα = 1.2 mas, decl. offset of Δδ = 0.2 mas, position angle of the disk major axis of P.A. = 28°, and disk inclination of i = 26° between the disk plane and plane of sky. Although such a simple model fitting cannot rule out other possibilities for the observed H30α kinematic pattern, it shows that disk rotation around a forming massive star can be a valid explanation.

While the velocity gradient of the ionized gas in source 3 could trace a rotating disk, it is more difficult to tell whether the velocity gradient seen in source 2 (Figure 6(e)) is real or not, or what its possible origin is if it is real, due to the low S/N. The direction of the tentative H30α velocity gradient in source 2 is at a position angle of P.A. ≈ 25°. The disk around source 2 has a position angle of P.A. ≈ − 25° (see Section 6; also Sánchez-Monge et al. 2013), and the ionized jet originating from source 2 has a position angle of P.A. ≈ 20° close to the source, as the jet has an S-shape morphology (Beltrán et al. 2016). Therefore the direction of the velocity gradient of the H30α line is consistent with the jet direction. However, the line-of-sight velocity directions seen in the H30α line (redshifted in the north and blueshifted in the south) are opposite of the outflows. Molecular line observations have shown that the blueshifted molecular outflow is toward the north and NE, while the redshifted outflow is toward the south and southwest (Gibb et al. 2003; Birks et al. 2006; Qiu et al. 2013; see Section 6.2). The atomic line observations show that the jet toward the north is blueshifted (Fedriani et al. 2019). The jet-associated H₂O masers have also been reported to be redshifted (V_lsr >30 km s⁻¹) toward the south. Therefore, despite the fact that the velocity gradient orientation is similar to the jet direction, it is difficult to confirm that the velocity gradient is mainly tracing outflow motion. One possibility is that the observed H30α velocity gradient is affected by the disk's rotation (redshifted on the northern side). In fact, in this source, the Brγ, an NIR HRL, does show one high-velocity component associated with the jet and a slow component associated with the disk (Fedriani et al. 2019), supporting that the ionized gas in this source is in both outflow and disk. Another possibility is that the jet is close to edge-on and therefore the inclination of the innermost part of the outflow is different from that of the outflows at larger distances.

As there are more and more 100-au-scale observations toward massive star-forming regions, kinematics of photoionized gas at disk scales around forming massive stars traced by recombination lines have started to be reported. This includes ionized disks (some are surrounded by neutral disks traced by molecular lines) and ionized outflows just launched from the disk (some are corotating with the disk) (e.g., Zhang et al. 2017; Zhang et al. 2019b, 2019c, Guzman et al. 2020; Jiménez-Serra et al. 2020; Moscadelli et al. 2021). These observations show the effectiveness of using recombination lines as a tracer of the ionized gas in the innermost region around forming massive stars to study accretion and photoionization feedback processes in massive star formation.

For sources 1−4, Figure 5 shows that most SEDs in the VLA bands have spectral indices between 0 and 2, indicating partially optically thick free–free emission. The optical depth of the free–free emission is dependent on the frequency (see the Appendix). For an ionized region with a uniform electron density, a spectral index of α ≈ 1 is only seen in a relatively narrow transition frequency range. However, in our detected sources, the intermediate spectral indices are seen over a wide range of frequencies, which requires a density gradient in the ionized gas (e.g., Keto 2003).

We construct a simple model to explain the observed free–free continuum. The total H30α fluxes are also used to constrain the model. The model details are introduced in the Appendix. In the model, we assume the emitting region is a circular disk on the plane of sky, and the emission measure ($\mathrm{EM}=\int {n}_{e}^{2}{dl}$) follows a power-law dependence on the radius ($\mathrm{EM}\propto {r}^{{\rho }_{\mathrm{EM}}}$). The three free parameters of the model are the electron temperature of the ionized gas (T_e ), the emission measure at the radius of 10 au (EM_10au), and the power-law index of the emission measure distribution ρ_EM. The fitting results are listed in Table 3. Below we discuss the fitting results of each source in detail.

Table 3. Parameters of the Best-fit Models of Continuum SED and H30α Flux

Source	Te	EM10au	ρEM	ne,10au	αdust	${\dot{N}}_{{\rm{i}}}$ a	m*,ZAMS	Spectral Type	L*,ZAMS
	(104 K)	(109 pc cm−6)		(107 cm−3)		(1046 s−1)	(M⊙)		(103 L⊙)
Source 1	1.0	0.2	2	0.2	2.3	0.009–1	7.3–12.0	B2–B0.5	1.8–9.0
Source 2	0.9	20	3	2		0.8–48	11.7–18.4	B0.5–O9.5	8.4–31
Source 3	0.6	6.3	3	1.1		0.3–30	10.7–17.3	B1–B0	6.2–26
Source 4	1.2	0.2	2	0.2		0.013	7.6	B2	2.1

Note.

^a As ${\dot{N}}_{{\rm{i}}}$ is sensitive to ρ_EM, we show a range of ${\dot{N}}_{{\rm{i}}}$ if ρ_EM is not well constrained by the fitting (see text for details).

Download table as:  ASCIITypeset image

5.1. ALMA Source 2

For source 2, we simultaneously fit the VLA continuum fluxes (assuming they are free–free dominated) and H30α flux. Figure 8 shows the best-fit model and the χ² distribution in the parameter space. The 1.3 mm continuum flux is not used in the fitting as it contains significant dust contribution. The ionized gas temperature and emission measure are both well constrained by the fitting, with the best-fit model having T_e = 9000 K and EM_10au = 2 × 10¹⁰ pc cm⁻⁶. The power-law index of the emission measure gradient is not well constrained. The best-fit model has ρ_EM = 3, but all ρ_EM values between 2 and 3 can produce similarly good fits. Note that ρ_EM = 2 corresponds to a power-law index of ${\rho }_{{n}_{e}}=1.5$ for the electron density gradient, which is consistent with a photoionized disk, while ρ_EM = 3 corresponds to a power-law index of ${\rho }_{{n}_{e}}=2$ in the electron density gradient, which is consistent with an expansion flow (see the Appendix). The best-fit model with ρ_EM = 3 is therefore consistent with the wide recombination line (see Section 4). The derived EM_10au = 2 × 10¹⁰ pc cm⁻⁶ also leads to an electron density of n_e ≃ 2.0 × 10⁷ (l/10au)^−1/2 cm⁻³ at 10 au scale (l is the line-of-sight scale of the ionized region). The compact size (<100 au, seen from the image) of the ionized region, the derived emission measure and electron density, as well as the H30α line width of 80 km s⁻¹, are all consistent with a HCH ii region (e.g., Hoare et al. 2007).

Zoom In Zoom Out Reset image size

From the fitted parameters T_e , EM_10au, and ρ_EM, we further estimate the ionizing photon rate ${\dot{N}}_{{\rm{i}}}$ (see the Appendix). Since ${\dot{N}}_{{\rm{i}}}$ is sensitive to ρ_EM (see the last term in Equation (A19)) and ρ_EM is not well constrained, the models with similarly good fit give a range of ${\dot{N}}_{{\rm{i}}}=(0.8-48)\times {10}^{46}\,{{\rm{s}}}^{-1}$, with ${\dot{N}}_{{\rm{i}}}=48\times {10}^{46}\,{{\rm{s}}}^{-1}$ being the best model (Figure 8(e)). For zero-age main sequence (ZAMS) stars, this corresponds to a stellar mass of 11.7–18.4 M_⊙ (Davies et al. 2011), i.e., the spectral type B0.5–O9.5 (Mottram et al. 2011). Our estimates of the emission measure, ionizing photon rate, and stellar mass, are higher than that estimated by Beltrán et al. (2016) (EM −2.7 × 10⁸ pc cm⁻⁶, ${\dot{N}}_{{\rm{i}}}=1\times {10}^{45}\,{{\rm{s}}}^{-1}$, and B1-type ZAMS star). This is because their estimation is based on the assumption of optically thin free–free emission, and a uniform density distribution. In our model, we consider a gradient in the emission measure (and therefore the optical depth of free–free emission), which leads to higher values of ionizing photon rate and stellar mass. We note that the derived ZAMS mass should be treated as a reference mass because a forming massive star that is still in a bloated phase before reaching the main sequence requires more mass than a ZAMS star to produce the same amount of ionizing photons (e.g., Hosokawa et al. 2010). Detailed protostellar evolution calculations taking into account various accretion histories of massive star formation show that accreting massive young stars with a wide range of mass have similar levels of ionizing photon rates (e.g., Zhang et al. 2014; Tanaka et al. 2016; Zhang & Tan 2018), which we will discuss in more details in Section 5.5.

The best-fit model gives a free–free flux of 10 mJy at 1.3 mm, leaving about 25 mJy to be the dust emission at 1.3 mm. Such a 1.3 mm dust emission flux leads to about 0.9 mJy of dust emission at 7 mm assuming an optically thick dust spectral index of 2. The dust component will be much lower if an optically thin dust spectral index of 3 is used. This confirms that the continuum at 7 mm (and at lower frequencies) is indeed dominated by the free–free emission, and that the method of using a free–free continuum model to fit the VLA band continuum is valid.

5.2. ALMA Source 3

As shown in Section 4, source 3 has strong H30α emission. The H30α flux integrated over area and velocity is 1.7 Jy km s⁻¹. Under the local thermal equilibrium (LTE) condition, such an H30α flux gives an optically thin limit of 1.3 mm free–free continuum ranging from 13 mJy with T_e = 6000 K to 51 mJy for T_e = 20,000 K (Equation (A11)). Such a free–free flux is higher than the observed total continuum flux at 1.3 mm, except for the low end of T_e . If the continuum is not optically thin, the intensity ratio between the recombination line (after subtracting the continuum) and the free–free continuum becomes even lower than that estimated in Equation (A11), which leads to an even higher free–free continuum flux. Therefore, the LTE condition for the H30α of this source may be in question, i.e., the observed bright HRL emission in this source might be enhanced by maser emission. To date, millimeter HRL masers have been reported in a modest number of sources (MWC349A, Martín-Pintado et al. 1989; Jiménez-Serra et al. 2013; MonR2-IRS2, Jiménez-Serra et al. 2013; and G45.47+0.05, Zhang et al. 2019c). Here, we consider this to be another case of a possible millimeter HRL maser.

As we do not have a non-LTE effect in our simple model, we use the observed H30α flux as an upper limit to constrain the model. Figure 9 shows the best-fit model and the χ² distribution in the parameter space. As neither the 1.3 mm continuum flux nor the H30α flux is used in fitting, we cannot constrain well the free–free emission at millimeter wavelengths. The emission measure is relatively well constrained to EM = 6.3 × 10⁹ pc cm⁻⁶ by the fitting, which leads to an electron density of n_e ≃ 1.1 × 10⁷(l/10au)^−1/2cm⁻³ at the 10 au scale. The compact size, the derived emission measure and electron density, as well as the H30α line width of 40 km s⁻¹ are also consistent with a HCH ii region. The best-fit model has T_e = 6000 K and ρ_EM = 3, but the dependence of χ² on T_e and ρ_EM are relatively weak. As Figure 9(e) shows, the models with similarly good fit give a range of ${\dot{N}}_{{\rm{i}}}=(0.3-30)\times {10}^{46}\,{{\rm{s}}}^{-1}$, which corresponds to a ZAMS stellar mass of 10.7–17.3 M_⊙ (Davies et al. 2011), i.e., the spectral type B1−B0 (Mottram et al. 2011). The best model has ${\dot{N}}_{{\rm{i}}}=2\times {10}^{47}\,{{\rm{s}}}^{-1}$, corresponding to a ZAMS stellar mass of 16.4 M_⊙.

Zoom In Zoom Out Reset image size

5.3. ALMA Source 1

No H30α emission is detected in the ALMA source 1 (core A or VLA source 2). The noise level of the integrated H30α emission map is 0.06 Jy beam⁻¹ km s⁻¹. We use the 3σ level (0.18 Jy km s⁻¹) as the upper limit for the H30α flux for this source. However, our free–free continuum and H30α model failed to reproduce the VLA band continuum fluxes and the H30α upper limit. The fact that the H30α emission is weak and the spectral index between 1.3 and 7 mm is greater than 2 (see Figure 5), indicates that the dust emission may be significant even at 7 mm. Therefore, we include an additional dust component in the SED fitting (Equation (A16)). We assume that the ALMA 1.3 mm continuum is dominated by dust emission and use the dust spectral index α_dust as an additional free parameter. The best-fit model and the χ² distribution in the parameter space are shown in Figure 10. Note that since we only have five data points and one upper limit to constrain four free parameters, the constraints are relatively weak.

Zoom In Zoom Out Reset image size

The emission measure and the dust spectral index are relatively well constrained. The best-fit model has EM_10au = 2 × 10⁸ pc cm⁻⁶ (corresponding to an electron density of n_e ≃ 2.0 × 10⁶ (l/10au)^−1/2 cm⁻³) and α_dust = 2.3. On the other hand, the electron temperature T_e and the emission measure gradient power-law index ρ_EM are not constrained. The models with similarly good fit give an ionizing photon rate range of ${\dot{N}}_{{\rm{i}}}=(0.09-10)\times {10}^{45}\,{{\rm{s}}}^{-1}$, which corresponds to a ZAMS stellar mass of 7.3–12.0 M_⊙ and spectral type B2–B0.5. Interestingly, the derived emission measure and electron density are not very consistent with typical HCH ii regions. However, the compact size seen at radio wavelengths, the low free–free continuum flux, and the derived low central stellar mass, suggest that this is an extremely young photoionization region.

5.4. ALMA Source 4

For source 4, the SED in the VLA frequency range is almost flat, indicating nearly optically thin free–free emission. Adopting the same method used for sources 2 and 3, we fit the three VLA bands with constraints for the upper limit of the H30α integrated flux (3σ = 0.18 Jy km s⁻¹). The results are shown in Figure 11. Note that we do not include the 6 cm continuum in the fitting, as it is significantly higher than the other VLA continuum fluxes (see Figure 5(d)), which may be caused by nonthermal emission, which has negative indices. The ALMA 1.3 mm continuum is also not used because it has a significant contribution from dust emission. The best-fit model has T_e = 12,000 K, EM_10au = 2 × 10⁸ pc cm⁻⁶, and ρ_EM = 2. These parameters lead to an ionizing photon rate of ${\dot{N}}_{{\rm{i}}}=1.3\times {10}^{44}\,{{\rm{s}}}^{-1}$, which corresponds to a ZAMS stellar mass of 7.6 M_⊙ and spectral type B2.

Zoom In Zoom Out Reset image size

However, as discussed in Section 3.2, while source 4 could be a forming massive star creating a small ionized region via photoionization like the fitting model assumes, its low 1.3 mm continuum emission (on the same level as other low-mass sources in the region) suggests that it could also be a low-mass protostar with free–free emissions caused by shock-induced ionization from an unresolved jet. In fact, the negative spectral index in the low-frequency range could indicate the contribution of nonthermal emission caused by the jet shock. Another supporting evidence is that if the free–free emission is caused by photoionization, the mass ratio between sources 2 and 4 is 2–3 as estimated above. If sources 2 and 4 are indeed a binary system, the existence of such a massive companion should have a significant impact on the circumbinary disk structure. However, as we will show in Section 6.1, the disk structure is highly symmetric with respect to source 2 indicating little to no impact from source 4.

5.5. Consistency of Stellar Mass Estimation

Our observations include various molecular lines tracing different gas structures in the region. Here, we present results for the SO₂(22_2,20 − 22_1,21), H₂CO(3_2,1 − 2_2,0)), CH₃OH(4_2,2 − 3_1,2; E), H₂S(2_2,0 − 2_1,1), and CH₃OCH₃(13_0,13 − 12_1,12) lines, which trace dense gas, and ¹²CO(2-1) and SiO(5-4), which trace outflows. The observing parameters for these lines are listed in Table 4. These lines are the main targeted lines of our spectral window setups. While some other weaker lines are also detected in the same spectral windows, we defer the analysis of the full molecular line data to a future study.

Table 4. Parameters of the Observed Lines ^a

Species	Transition	Frequency	Eu /k	S μ2	Velocity Resolution	Synthesized Beam	Channel rms
		(GHz)	(K)	(D2)	(km s−1)		(mJy beam−1)
H	30α	231.90090			0.63	0039 × 0031 (P. A. = − 680)	1.7
SO2	222,20 − 221,21	216.6433035	248	35.3	0.68	030 × 029 (P. A. = − 35)	3.3
						0052 × 0040 (P. A. = − 639)	1.9
CH3OH	42,2 − 31,2; E	218.4400630	45.5	13.9	0.68	029 × 028 (P. A. = − 59)	5.6
						0052 × 0040 (P. A. = − 643)	2.0
H2CO	32,1 − 22,0	218.7600660	68.1	9.06	0.67	029 × 028 (P. A. = − 37)	5.4
H2S	22,0 − 21,1	216.7104365	84.0	2.06	0.68	030 × 029 (P. A. = − 49)	3.4
SiO	5 − 4	217.1049190	31.3	48.0	0.67	030 × 029 (P. A. = − 86)	2.9
CH3OCH3	130,13 − 121,12; AA	231.9877830	69.8	170	0.63	028 × 027 (P. A. = − 365)	2.9
	130,13 − 121,12; EE	231.9878580	69.8	272
	130,13 − 121,12; AE	231.9879320	69.8	102
	130,13 − 121,12; EA	231.9879330	69.8	68.0
12CO	2 − 1	230.5380000	16.6	0.0242	0.63	028 × 027 (P. A. = − 141)	4.0

Note.

^a Molecular line information taken from the CDMS database (Müller et al. 2005).

Download table as:  ASCIITypeset image

6.1. Multi-scale Accretion Structures

Figure 12 shows the integrated emission maps of the SO₂, H₂CO, CH₃OH, SiO, H₂S, and CH₃OCH₃ lines. Here, we focus on the molecular lines associated with the central sources, especially source 1 (core A), and sources 2 and 4 (core B). For source 1, the line emission is mostly concentrated near the source center (except for the SiO emission whose peak is more offset). However, these lines show different spatial distributions around source 2. Note that sources 2 and 4 are proposed to form a binary system by Beltrán et al. (2016) based on their closeness and mass estimates. However, it is still unclear whether they are a true bound binary system, or how they relate to the nearby other sources identified in 1.3 mm long-baseline data.

Zoom In Zoom Out Reset image size

The SO₂ emission around source 2 is concentrated toward the center, and shows a morphology consistent with that of two spiral structures that are symmetric with respect to the center (Figure 12(a)). In contrast to the SO₂ line, the center is devoid of H₂CO and CH₃OH emissions. While the spiral structures are still distinguishable in the H₂CO and CH₃OH emission, they mainly show an elliptical ring around sources 2 and 4. In addition, the H₂CO and CH₃OH lines further show a stream of emission to the north of source 2 that connects to source 1, forming a larger elliptical shape surrounding sources 1 and 2. The aspect ratio of this larger elliptical structure is about 2:1. If we assume that this structure is circular, then the inclination between its midplane and the line of sight is about 30°. This is close to the derived inclination of the disk around source 2 estimated from previous observations (Sánchez-Monge et al. 2013). Therefore, despite the fact that material is distributed in the filament-like structure on a larger scale, there may be a near-circular structure with a radius of about 4000 au around the two main forming massive stars (sources 1 and 2), which then feeds into two separate disks feeding each source. In addition, the disk feeding source 2 (and source 4) has formed spiral structures.

The kinematics of these lines are consistent with such a scenario. Figures 13 and 14 show the channel maps of the SO₂ and CH₃OH lines around these sources. The spiral structures around source 2 seen in the integrated maps have highly symmetric and ordered velocity structures, with the southeastern arm blueshifted and the northwest (NW) arm redshifted. Source 1 (core A) is slightly more redshifted than source 2 (core B), and the stream connecting these two sources (Figure 14) is redshifted with respect to source 2 and joins the source 1 disk at a velocity slightly more redshifted than the source 1 systemic velocity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Compared to the H₂CO and CH₃OH lines, the H₂S and CH₃OCH₃ emissions are more concentrated around the two main sources, but their morphologies are similar to the H₂CO and CH₃OH lines. SiO emission around source 2 appears to follow the spiral arm structures seen in the SO₂ line. However, SiO emission is widely spread out in this region with multiple components tracing outflows (originated from source 1) and emission along the filament, as we will discuss in more detail in Section 6.2. Therefore, it is difficult to conclude that the SiO emission around source 2 actually traces the disk spiral structure.

6.1.1. Velocity Structure around Source 2

Figure 15 shows the position–velocity (PV) diagrams along a position angle of P.A. = −25° across source 2 to show the kinematics of its proposed disk in these molecular lines. This position angle is consistent with the position angle of the disk derived from the previous 0.85 mm observation at a lower angular resolution (Sánchez-Monge et al. 2013). The spiral structures are also roughly aligned along this direction. Velocity gradients consistent with rotation are seen in all molecular lines. However, only in the SO₂ line, a high-velocity component reaching a rotation velocity of about 10 km s⁻¹ is seen within about 02 (440 au) from the center. This high-velocity component is missing in all other molecular lines shown here. On the other hand, the low-velocity component of SO₂ is highly consistent with the other lines. The high-velocity inner component seen in the SO₂ line may be caused by an inner disk rotating faster.

Zoom In Zoom Out Reset image size

There are two reasons why the SO₂ line can better trace the inner region and highlight the spiral structure. First, the observed SO₂ transition has a higher upper energy level(E_u /k = 248 K) than the other lines shown here (see Table 4). Therefore it is more sensitive to the warmer inner region while the other molecular lines are more dominated by the outer cooler area. Second, as a typical shock tracer, the SO₂ emission may be enhanced by shocks associated with the spiral arms, and therefore more concentrated along the spiral arms than in the other species. Such behavior is very similar to what has been reported in several other massive protostellar disks, such as G339.88-1.26 (Zhang et al. 2019b) and IRAS16547-4247 (Tanaka et al. 2020). In those cases, SO₂ lines are also found to trace the inner region in the disk/envelope system with faster rotation than molecular lines such as CH₃OH and H₂CO. However, in those sources, SiO emission is found to trace the innermost part of the disk with the highest rotation velocities, which is not seen in the case of G35.2.

6.1.2. Fitting the Spiral Structures around Source 2

To better understand the properties of the spiral structure seen in the SO₂ line, we fit them with two spiral shapes; one is a logarithmic spiral with the form of R = R₀ e^{b θ} , and another is an Archimedean spiral with the form of R = R₀ + b θ (e.g., Huang et al. 2018; Kurtovic et al. 2018).

Figure 16(a) shows the integrated SO₂ emission in terms of distance from source 2 and polar angle. We determine the polar angle of the two spiral arms by their emission peaks at each distance on the emission polar map. Only the radius range between 02 and 08 is used as the spiral structure is most clear in this range. Figure 16(b) shows the fit to the determined spiral shapes. The parameters of the best-fit shapes are listed in Table 5. The pitch angles of the spirals $\phi ={\tan }^{-1}(\tfrac{{dR}}{{Rd}\theta })$ are also listed. It appears that the logarithmic shape fits the spiral arms better than the Archimedean shape, reproducing the curved shape (also see Table 5). The pitch angle is 69° for the NW arm and 62° for the southeast (SE) arm. Note that the two spiral arms lie along the major axis of the disk, the large pitch angles are partly due to the high inclination (closer to edge-on view) and the low spatial resolution.

Zoom In Zoom Out Reset image size

Table 5. Parameters of the Best-fit Spirals

Spiral Arm	Logarithmic (R = R0 eb θ )				Archimedean (R = R0 + b θ)
	R0	b	Pitch Angle a	χ2	R0	b	Pitch Angle	χ2
	(arcsec)	(deg−1)	(deg)		(arcsec)	(arcsec deg−1)	(deg)
NW	0.042	0.044	69°	0.77	−0.69	0.022	81°–57°	57
SE	0.00032	0.032	62°	3.9	−2.96	0.015	78°–48°	48

Note.

^a Pitch angle from 02–08.

Download table as:  ASCIITypeset image

There are commonly two mechanisms to form the spiral structures in a disk, gravitational instability and perturbation by a companion or a planet (e.g., Forgan et al. 2018; Huang et al. 2018). A pair of symmetric logarithmic spiral arms are expected if the spiral structures are formed due to gravitational instability in the disk (e.g., Forgan et al. 2018). Also, the disk around source 2 is estimated to be about 3 M_⊙ (Sánchez-Monge et al. 2013), which is relatively massive with respect to the mass of the central sources (sources 2 and 4 have a dynamical mass of about 18 M_⊙ in total; see Section 5.5). Such a massive disk is prone to gravitational instability. Another possibility is that the spiral structures are caused by the perturbation of the binary companion, i.e., source 4. There are also several lines of evidence supporting such a scenario. First, the NW spiral arm shows a local minimum close to the location of source 4 (see Figure 16(a)), which may indicate that the binary companion has opened a gap in the spiral. Furthermore, despite that a single Archimedean spiral cannot well fit the NW arm, this can be fit with two Archimedean shapes with a discontinuity at r ≈ 03, which leads to a local maximum of pitch angle. In fact, simulations have shown that the spirals formed by perturbation of a companion have an increasing pitch angle toward the position of the companion (e.g., Zhu et al. 2015; Bae & Zhu 2018), which is similar to what we observed in the NW arm. Therefore, with the current data, it is difficult to determine unambiguously which is the formation mechanism of the spiral structures in the source 2 disk.

However, it is also possible that the rotation structure around source 2 on the 2000 au scale is not truly a disk. The extent of the spiral structures (∼2000 au) is much larger than the size of typical Keplerian disks around massive protostars seen in high-resolution observations, which have typical radii ≲ a few × 10² au (e.g., Ginsburg et al. 2018; Ilee et al. 2018; Zhang et al. 2019a; Maud et al. 2019). Furthermore, although the kinematics of the rotation structure around source 2 was explained by Keplerian rotation with lower-resolution data (Sánchez-Monge et al. 2013), its kinematics can also be explained by rotating-infalling material in the outer part feeding an inner small disk (see Figure 15). The high-velocity component shown in the SO₂ PV diagram (Figure 15(a)) has a radius of ∼500 au, which is more consistent with other Keplerian disks in massive star formation. If this is the case, the spiral structures can be the infalling streamers feeding the inner disk. Such streamers on several × 10³ au scales were reported in other sources (e.g., Pineda et al. 2020), and can be explained by accretion of unsmooth turbulent material or even cloudlets (e.g., Kuffmeier et al. 2019; Hanawa et al. 2022). However, it is unclear whether such a scenario can explain the high symmetry seen in the pair of spiral structures.

6.1.3. Velocity Structure around Source 1

Rotation velocity gradients are also seen around source 1. In the low-resolution images, the emission around source 1 appears to be blueshifted in the NW and redshifted in the SE, opposite from that around source 2 (see Figure 17). Similar velocity gradient directions were reported by Sánchez-Monge et al. (2014) with ALMA 0.85 mm observations. However, with higher resolution, the innermost region (<01, ∼200 au) around source 1 shows a different velocity gradient direction, as shown by panels (c) and (d) in Figure 17. The innermost SO₂ emission is blueshifted in the south and redshifted in the north. The CH₃OH line has blueshifted emission in the SE. There is some redshifted emission in the NW of source 1 in the innermost region, but is affected by absorption. The velocity gradients in the innermost region around source 1 revealed by these two molecules are opposite from that in the outer region but generally consistent with that seen around source 2. This rotation direction is also consistent with the relative velocities of sources 1 and 2, as source 1 (the NW source) is redshifted with respect to source 2.

Zoom In Zoom Out Reset image size

Figure 18 shows PV diagrams of SO₂ and CH₃OH lines, along a direction with a maximum velocity gradient across source 1 seen in the SO₂ moment 1 map (P.A. = 0°). In the SO₂ PV diagram, there are higher-velocity components with clear gradients in the innermost region that can be explained by rotation. However, in the case of pure rotation (e.g., a Keplerian disk), the blueshifted emission would be limited to one side while the redshifted emission would be limited to the other side, which is not compatible with the observed pattern. Instead, a combined motion of infall and rotation can reproduce the observed kinematics. Compared to the SO₂ emission, the CH₃OH emission shows a similar increase in velocity gradient toward the center. However, it does not show the central high-velocity component as SO₂, indicating that the CH₃OH traces slightly outer regions than SO₂, similar to the case of source 2 (Section 6.1.1). The CH₃OH line shows strong absorption against the bright continuum emission in redshifted velocities (see Figure 18(b)). This indicates infalling material in front of the central source along the line of sight, consistent with a combined motion of infall and rotation. Furthermore, if the accreting material is not evenly distributed (e.g., there are unsmooth accretion flows), the combined motion of infall and rotation can cause the change of velocity gradient directions from the outer region to the inner region (e.g., Liu 2017).

Zoom In Zoom Out Reset image size

To demonstrate the possibility of the combined motion of infall and rotation, we construct a simple model to fit the observed PV diagrams. Following the method used by Sakai et al. (2014a) and Zhang et al. (2019a), we fit the observed PV diagrams with a simple model of infalling rotation with the rotation velocity v_φ and infall velocity v_r described as

$$\begin{eqnarray}&&{v}_{\varphi }(r)={v}_{\mathrm{CB}}\displaystyle \frac{{r}_{\mathrm{CB}}}{r},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{v}_{r}(r)=-{v}_{\mathrm{CB}}\displaystyle \frac{\sqrt{{r}_{\mathrm{CB}}(r-{r}_{\mathrm{CB}})}}{r}.\end{eqnarray} \tag{ 2 }$$

Such motion conserves both angular momentum and mechanical energy. r_CB is the innermost radius that such infalling gas can reach with the angular momentum conserved (i.e., the centrifugal barrier), where v_φ = v_CB and v_r =0.v_CB is the rotational velocity at the centrifugal barrier. In such a model, the central mass is ${m}_{* }={r}_{\mathrm{CB}}{v}_{\mathrm{CB}}^{2}/(2G)$. Since we only focus on the kinematics pattern here, we adopt a simple geometry and density distribution in the model. We assume the disk has a height of h(r) = 0.2r on each side of the midplane, and the density distribution follows ρ(r) ∝ r^−1.5 (e.g., Oya et al. 2016). For simplicity, we also assume the emissions are optically thin and the excitation conditions are uniform across the region. Note that the density distribution, excitation conditions, and disk height profiles have a combined effect on the detailed emission pattern, but do not affect the general kinematics pattern, which we aim to explore here. We also fix the outer radius to be 02 (4.4 × 10² au), as we only focus on the kinematics of the innermost region. Therefore, we have in total three free parameters: the central mass m_*, the radius of the centrifugal barrier r_CB, and the inclination angle i.

To obtain the best-fit model, we compare the model PV diagram with the observed PV diagram of the SO₂ line. We focus on the SO₂ line as the CH₃OH line does not trace the innermost region well. We explore the inclination angle i with values ranging from 0° to 40° with an interval of 10°, the angular radius of the centrifugal barrier r_CB/d in a range of 001–005 with an interval of 001, and the central dynamical mass in a range of 4–12 M_⊙ with an interval of 1 M_⊙, and approximately determine the best-fit model by eye. The SO₂ PV diagram can be well reproduced by the model with r_CB/d = 003 (r_CB ≈ 70 au), i = 30° between the line of sight and the disk midplane, m_* = 6 M_⊙ (Figure 18(a)). In Figure 18(b), we show a comparison of this model and the CH₃OH emission. The CH₃OH emission is consistent with the model in the outer part, but does not show the innermost high-velocity component. This confirms the scenario that CH₃OH traces a region with slightly larger radii than that traced by the SO₂ emission.

Our model fitting shows that infalling rotation can be a viable explanation for the gas kinematics in the innermost region around source 1. We note that we cannot rule out all the other possibilities. Especially, our data have limited sensitivities for the extended emissions, and the rotation direction changes from larger scale to the small scale, therefore our simple model cannot model the kinematics of more extended emission. Nevertheless, the model fitting suggests that the rotating structure around source 1 is not rotationally supported. The derived central mass of 6 M_⊙ is close to that estimated from the ionizing photon rate (see Section 5), and higher than the previously derived dynamical mass of 4 M_⊙ (Sánchez-Monge et al. 2014). In the previous estimation, a Keplerian disk was assumed, higher dynamical mass is expected if the disk is not rotationally supported. Despite the non-detection of the rotationally supported disk in our observation, such a disk may exist inside the centrifugal barrier, i.e., within a radius of 003 (∼70 au). The resolution and line coverage of our current observation may not be able to effectively probe such a disk.

Overall, the hierarchical accretion structure around sources 1 and 2, in which both sources are fed by a common rotating structure, and each of them is fed by the individual rotating structure, indicates that the accretion toward these two sources is still highly active and the start of photoionization has not stopped the accretion.

6.2. Molecular Outflows

Figure 19 shows the ¹²CO channel maps of G35.2. The ¹²CO emission shows a complex morphology. Our short observations with limited u − v coverage cannot effectively recover the extended component of the emission, causing significant missing flux, which makes it even more difficult to analyze the outflows in this region.

Zoom In Zoom Out Reset image size

In the blueshifted velocities (with respect to the systemic velocity, which is around V_lsr = + 30 to + 35 km s⁻¹; see Section 6.1), the most distinguishable outflow feature is a wide arc-shaped structure toward the NE. This feature is best seen in velocity channels V_lsr ≲ + 16 km s⁻¹, but can also be seen in less blueshifted velocities down to V_lsr = + 24 km s⁻¹. The wide-angle CO outflow toward the NE has been seen in previous observations, including single-dish observations (e.g., Gibb et al. 2003; Birks et al. 2006). This structure is also seen in the NIR both in molecular and atomic emission (Caratti o Garatii et al. 2015; Fedriani et al. 2019). On the other hand, in the redshifted velocity channels, the emission appears to be more concentrated in a few collimated components, which are better shown in the integrated emission map shown in Figure 20. We identify in total four collimated outflow components, which are labeled by dashed lines with different colors in Figure 19 (in the channel of V_lsr = + 52 km s⁻¹) and Figure 20. Assuming that they are separate outflows, we identify their driving sources by comparing their locations with the positions of the continuum sources.

Zoom In Zoom Out Reset image size

Outflow 1 (marked with the black dashed line) is most clearly seen in the channels of + 40 km s⁻¹ ≤ V_lsr ≤ + 52 km s⁻¹, but it can also be seen in slightly blueshifted channels of V_lsr = + 28 km s⁻¹. At V_lsr = + 48 km s⁻¹, it is clear that there is a string of bright CO blobs associated with this outflow. Despite the fact that this outflow component does not have a well-defined structure in our data, we believe it traces the same N–S outflow component identified in previous CO observations (e.g., Birks et al. 2006), and is consistent with the radio jet (Gibb et al. 2003; Beltrán et al. 2016) and the bright outflow cavity seen in the NIR to MIR (De Buizer 2006; Zhang et al. 2013; Fedriani et al. 2019). This outflow appears to originate from source 2 (core B).

There appears to be an outflow associated with source 1 (core A), which is most clearly seen as two bright blobs to the NE of source 1 at velocities + 56 km s⁻¹ ≤ V_lsr ≤ + 64 km s⁻¹. This outflow has a position angle of about 30°, if the blobs toward the opposite direction are also associated with it (marked with the gray dashed line in Figures 19 and 20). This outflow is also traced by the SiO emission (Figure 20), which is best seen as a high-velocity redshifted emission blob (up to V_lsr > + 50 km s⁻¹) originating from source 1 in the same direction as the ¹²CO outflow. This SiO outflow has been reported by Sánchez-Monge et al. (2014). Its blueshifted emission is more extended and mostly resides along the filamentary continuum structure. It is possible that, in addition to the outflow, the SiO emission also traces the shocks caused by the gas flows forming the filament.

Outflow 3 (marked with the blue dashed line) originates from source 14 (core C). It is clearly seen at velocities V_lsr ≥ + 48 km s⁻¹. To the NE of source 14, there is a high-velocity component in this outflow that can be seen from V_lsr = + 68 km s⁻¹ up to V_lsr > + 90 km s⁻¹. This outflow has a position angle of 70°. The position angles of outflows 2 and 3 are highly consistent with the two larger-scale outflows seen in NIR shocked H₂ emission (Caratti o Garatii et al. 2015). Therefore, they are likely the same outflows seen in different wavelengths and different tracers. If this is the case, the wide CO outflow toward the NE should be a combination of these two outflows rather than a single wide outflow. Furthermore, if these two outflows are driven by sources 1 and 15, they are separated from the N–S outflow (outflow 1) driven by source 2. This suggests that the large NE outflow is not caused by the precession of the N–S outflow/jet, but rather, they are independent outflows driven by different sources.

Finally, we identify a fourth outflow (marked by a red dashed line in Figure 19) at a position angle of 110° originating from source 17 (core E). This outflow is best seen to the east of source 17 at velocities of + 48 km s⁻¹ ≤ V_lsr ≤ + 60 km s⁻¹, but can also be seen at high outflowing velocities (V_lsr = + 80 km s⁻¹, with a slightly different position angle). This outflow has not been reported by previous studies.

Note that outflows 3 and 4 are only seen clearly in the high-velocity redshifted channels. The crowded nature of this region, as well as the fact that the observation is very shallow with a lot of missing flux and side-lobe effects, make it very difficult to identify the different outflow structures in all the velocity channels. Therefore, with the current data, it is hard to conclude if the corresponding blueshifted components of outflows 3 and 4 are really missing or they are not seen because of the observation difficulties.

As presented in Section 3, we identify 22 compact sources from the ALMA long-baseline continuum image. As discussed in Section 5, the emission from 6 cm to 1.3 mm in sources 1−3 suggest that they are massive young stars that have created small photoionized regions around them. A similar scenario may be also true for source 4, but it can also be a low-mass protostar with a shock-ionized region caused by jet activity. Other sources are assumed to be low-mass members in the forming cluster. Following previous similar studies (e.g., Plunkett et al. 2018; Busquet et al. 2019), we use several methods to characterize their distribution.

We assume that the detected long-baseline (C9) 1.3 mm dust continuum emission in these sources comes from their disks and/or inner envelope (referred to as disks for simplicity), and therefore use the continuum fluxes to estimate the masses of such structures of the sources using the following equation:

$$\begin{eqnarray}&&M=\displaystyle \frac{{d}^{2}{F}_{\mathrm{dust},1.3\mathrm{mm}}}{{\kappa }_{1.3\mathrm{mm}}{B}_{1.3\mathrm{mm}}({T}_{\mathrm{dust}})},\end{eqnarray} \tag{ 3 }$$

where d = 2.2 kpc is the distance to the source, F_dust,1.3mm is the estimated dust continuum emission flux, and B_1.3mm(T) is the Planck function. An opacity of κ_1.3mm = 0.00899 cm² g⁻¹ is used (Ossenkopf & Henning 1994; for MRN dust with thin ice mantles and gas density of 10⁶ cm⁻³; a gas-to-dust mass ratio of 100 is assumed). The disk masses are listed in Table 2. For sources 1−4 (massive sources), a dust temperature of 100 K is assumed(suitable for typical massive sources; e.g., Zhang et al. 2014), and for sources 5−22, a dust temperature of 30 K is assumed (e.g., Tobin et al. 2013). Note that for sources 1−4, the free–free component estimated from modeling the SEDs and HRL fluxes (see Section 5) has been subtracted from the continuum fluxes. As the continuum fluxes are estimated with only the long-baseline image, the larger filament or accretion structure is filtered out, i.e., only the individual disk (and/or inner envelope) of each source is included in the mass estimation. Based on the assumed dust opacity, gas-to-dust ratio, and dust temperature of 30 K, the column density sensitivity is 6.3 g cm⁻² at the center of the field of view, which corresponds to a mass sensitivity of 3.8 × 10⁻³ M_⊙ for unresolved sources.

Figure 21(a) shows disk masses as a function of the projected distance from massive stars. Here, for any source, we calculate its distances to sources 1 and 2, and use the smaller value as its distance from massive stars. It appears that the sources within a distance of 3000 au of the massive stars have relatively low disk masses. Such a trend can be explained by the depletion of disks in the regions close to the massive stars. It may be due to the tidal/dynamic interactions between the sources in such a crowded region. In fact, our molecular line observations have revealed complex accretion structures on multiple scales that fit such a scenario (see Section 6.1). Another possibility of disk depletion close to the massive stars is photoevaporation by high UV radiation. However, this is less likely to be the main cause in our case. As discussed in Section 5, the massive sources in this cluster have not yet developed large H ii regions, and are still highly embedded inside dense disk/envelope structures. However, such a trend may be caused by other effects. Many of the sources close to the massive stars are associated with the larger filament and the hierarchical accretion structures (see Section 6.1). The long-baseline observations detect the small over-densities inside the larger structures associated with these sources. Whether such over-densities can represent the individual circumstellar disks is uncertain.

Zoom In Zoom Out Reset image size

Figure 21(b) shows the distribution function of the projected separations between pairs of sources in the cluster. We calculate the distribution function p(s) by counting the number of pairs (N_i ) with separation in the interval of (s, s + Δs), normalized by the total pair number (e.g., Cartwright & Whitworth 2004; Busquet et al. 2019)

$$\begin{eqnarray}&&p(s){\rm{\Delta }}s=\displaystyle \frac{2{N}_{i}}{{N}_{\mathrm{tot}}({N}_{\mathrm{tot}}-1)},\end{eqnarray} \tag{ 4 }$$

where N_tot = 22 is the total source number. Here, the separation s is normalized to the cluster radius ${r}_{\max }=18\buildrel{\prime\prime}\over{.} 1$ (4.0 × 10⁴ au), which is determined by the maximum separation between any source to the average position of all the sources. The distribution function shows a single peak at about 5'' (normalized s ≈ 0.25). As demonstrated by Cartwright & Whitworth (2004), a single peak in the distribution function is consistent with a smooth radial gradient of source density without subclustering. With N_tot = 22 and ${r}_{\max }=4.0\times {10}^{4}\,\mathrm{au}$, this cluster has a number density of protostellar sources of n = 7.2 × 10² pc⁻³. There are 20 sources within 2.5 × 10⁴ au (primary beam response >0.5), which yields a number density of n = 2.7 × 10³ pc⁻³.

We further use the minimum spanning tree (MST) method to characterize the clustering of these sources (Figure 21(c)). It shows the shortest path length connecting all the sources without including closed loops (Kruskal 1956). We construct the MST using the python package NetworkX (Hagberg et al. 2008). Following Cartwright & Whitworth (2004), from the constructed MST, we derive the Q parameter, which is used to distinguish between a smooth overall radial density gradient and multi-scale fractal subclustering. The Q parameter is defined as $\bar{Q}\equiv \bar{m}/\bar{s}$, where $\bar{m}$ is the mean MST path length normalized by the factor $\sqrt{{N}_{\mathrm{tot}}A}/({N}_{\mathrm{tot}}-1)$ with $A=\pi {r}_{\max }^{2}$ being the cluster area, and $\bar{s}$ is the mean separation of all the sources normalized to the cluster radius. From $\bar{m}=0.40$ obtained from MST, and $\bar{s}=0.43$ from the separation distribution (Figure 21(b)), we have $\bar{Q}=\bar{m}/\bar{s}=0.95$. $\bar{Q}=0.8$ is considered as a diagnostic boundary; higher values indicate a smooth overall radial density gradient and lower values indicate a fractal subclustering (Cartwright & Whitworth 2004). The result of $\bar{Q}=\bar{m}/\bar{s}=0.95$ is consistent with a cluster without subclustering, as the observed single-peaked separation distribution function has suggested.

Based on the MST method, we further calculate the mass segregation ratio Λ_MSR following the definition by Allison et al. (2009)

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{MSR}}=\displaystyle \frac{{L}_{\mathrm{norm}}}{{L}_{\mathrm{massive}}}\pm \displaystyle \frac{{\sigma }_{\mathrm{norm}}}{{L}_{\mathrm{massive}}},\end{eqnarray} \tag{ 5 }$$

where L_massive is the path length of the MST of the N most massive sources. L_norm and σ_norm are the average path length and its statistical deviation of the MST of N sources randomly selected from the cluster. To estimate L_norm and σ_norm, we randomly select 500 sets of sources for each N (Allison et al. 2009; Plunkett et al. 2018). Note that we exclude the outermost two sources from this analysis, as the low primary beam response near the edge of the field of view strongly biases against dimmer and therefore less massive sources.

Figure 21(d) shows Λ_MSR as a function of the number of most massive sources in the cluster used in constructing the MST analysis (N_MST). A value of Λ_MSR = 1 indicates a random distribution of sources without mass segregation, and Λ_MSR > 1 suggests that the N most massive sources are more concentrated with respect to the random sample, i.e., a mass segregation. The black data points show Λ_MSR calculated with sources ordered by their disk masses. The mass segregation becomes significant only for the most massive three sources and peaks at N_MST = 3, which corresponds to sources 1 and 2. However, for sources 1−4, the level of dust continuum emission is uncertain because the estimation of free–free contamination depends on the model (see Section 5). Meanwhile, these four sources should be the most massive sources as they have already formed photoionized regions. The red data points show Λ_MSR calculated by setting these four sources to be the most massive ones (ordered by their estimated stellar mass), followed by other sources ordered by the estimated disk masses. In this case, Λ_MSR clearly has a peak at N_MST = 4 indicating that these four massive sources show strong segregation from the other members of the cluster. No strong mass segregation is seen for N_MST ≳ 6, which corresponds to m_disk < 0.16 M_⊙. Considering the short formation timescale of massive stars and the fact that the photoionized regions of most massive sources in this cluster are still highly compact, the mass segregation is unlikely to be caused by dynamical relaxation (as is the case for more evolved clusters). However, instead, it should reflect the primordial mass segregation during the formation of the massive protostellar cluster.

We have presented ALMA 1.3 mm intermediate-resolution and long-baseline observations toward the massive star-forming region G35.20-0.74 (G35.2). Using the continuum emission and emission from different molecular lines, we characterized the massive and low-mass protostars in the region and their accretion and outflow structures. Our main conclusions are as follows:

• 1.  
  The 1.3 mm intermediate-resolution continuum image shows a string of cores along the filamentary structure, similar to what has been seen in the previous 0.85 mm observation. From the 1.3 mm long-baseline observation, we identified 22 compact continuum sources, suggesting a young forming cluster in the region. Among them, four sources have corresponding VLA detections from 6 cm to 7 mm. Three of them (ALMA sources 1–3) are consistent with massive young stars that have initiated photoionization, while the other source (ALMA source 4) could be a low-mass protostar with shock-induced free–free emission. Furthermore, we report the detection of a compact source at the location of source 3 in the 0.85 mm continuum, found as a result of reprocessing previously published 0.85 mm continuum data. The other 18 sources without corresponding VLA detections are considered to be lower-mass members in the cluster.
• 2.  
  Among the massive sources, two (sources 2 and 3) are found to have H30α recombination line emission. The H30α line kinematics shows ordered motions of the ionized gas in the innermost region. For source 3, the H30α kinematics are consistent with disk rotation, while for source 2 it is more consistent with the outflowing motion. We found evidence of potential maser activity in the H30α line emission of source 3, adding another candidate case to the handful of millimeter HRL maser discovered so far.
• 3.  
  We constructed models of photoionized regions to simultaneously fit the multiwavelength free–free fluxes and the H30α total fluxes. The model assumes an isothermal ionized region with a power-law distribution in the emission measure (and electron density), so that it can fit the observed SEDs, which show partially optical thick free–free emission over wide frequency ranges. The derived properties of the ionized regions are consistent with photoionized HCH ii regions.
• 4.  
  Molecular line emission shows multi-scale accretion structures around sources 1 and 2 (cores A and B). We have inferred that sources 1 and 2 are surrounded by a common rotating structure that connects to the individual disks. The source 2 disk is found to have two spiral arms, which may be caused by gravitational instability or perturbation of the binary companion source 4. As another possibility, the spiral structures could also be infalling streamers feeding an inner unresolved disk in source 2. The spiral arms are best seen in the SO₂ emission, which is enhanced by shocks along the spirals. SO₂ also traces the inner faster-rotating part of the disk while other molecular lines such as CH₃OH and H₂CO only trace the outer slower-rotation part. The source 1 disk appears to be more compact than the source 2 disk and not rotationally supported. Its velocity structure is consistent with infalling rotation of unsmooth accretion flows. The existence of such multi-scale accretion structures feeding sources 1 and 2 suggests that the accretion is still going on and is not stopped by the photoionization feedback.
• 5.  
  From the ¹²CO emission, we identified in total four outflows likely driven by four different sources (sources 1, 2, 14, and 17). By comparing with the outflows identified in previous studies with different tracers, we speculate that the N–S CO outflow driven by source 2 is consistent with the radio jet and the bright outflow cavity seen in NIR and MIR. The wide outflow toward the NE is actually composed of two separate outflows driven by source 1 and source 14, and it is not formed by the precession of the N–S jet. Furthermore, we identify another outflow in the east–west direction, driven by source 17.
• 6.  
  We analyzed the distribution of the identified compact continuum sources in the cluster. Assuming the compact dust continuum emission associated with these sources comes from their disks or inner envelopes, we analyzed the relation of their masses with their distances to the massive stars, and found potential disk depletion due to dynamical interactions in the inner region. The spatial distribution of the sources suggests a smooth overall radial density gradient without subclustering, and tentative evidence of primordial mass segregation.

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.01454.S, and ADS/JAO. ALMA#2017.1.00181.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Y.Z. acknowledges support from RIKEN Special Postdoctoral Researcher Program and JSPS KAKENHI grant JP19K14774. K.E.I.T. acknowledges support by JSPS KAKENHI grant Nos. JP19K14760, JP19H05080, JP21H00058, JP21H01145. J.C.T. acknowledges support from ERC project MSTAR, VR grant 2017-04522. Y.-L.Y. acknowledges the support from the Virginia Initiative of Cosmic Origins (VICO) Postdoctoral Fellowship. R.F. acknowledges funding from the European Unions Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 101032092. G.G. acknowledges support by the ANID BASAL project FB210003.

Software: CASA (http://casa.nrao.edu; McMullin et al. 2007), The IDL Astronomy User's Library (https://idlastro.gsfc.nasa.gov), NetworkX (https://networkx.org; Hagberg et al. 2008).

We construct a simple model to simultaneously fit the continuum flux densities observed in multiple VLA bands and ALMA 1.3 mm, as well as the total integrated flux density of the H30α line.

The optical depth of the free–free emission (e.g., Wilson et al. 2013) is

$$\begin{eqnarray}&&{\tau }_{\nu ,\mathrm{ff}}=8.235\times {10}^{-2}{\left(\displaystyle \frac{{T}_{e}}{{\rm{K}}}\right)}^{-1.35}{\left(\displaystyle \frac{\nu }{\mathrm{GHz}}\right)}^{-2.1}\left(\displaystyle \frac{\mathrm{EM}}{\mathrm{pc}\,{\mathrm{cm}}^{-6}}\right).\end{eqnarray} \tag{ A1 }$$

where T_e is the temperature of the ionized gas, and EM is the emission measure defined as $\mathrm{EM}\equiv \int {n}_{e}^{2}{dh}$. Assuming that the ionized gas is isothermal, the radiative transfer of the free–free emission gives the brightness temperature as

$$\begin{eqnarray}&&{T}_{\nu ,\mathrm{ff}}={T}_{e}\left(1-{e}^{-{\tau }_{\nu ,\mathrm{ff}}}\right).\end{eqnarray} \tag{ A2 }$$

The total free–free flux density is then

$$\begin{eqnarray}&&{S}_{\nu ,\mathrm{ff}}=\displaystyle \frac{2k{\nu }^{2}}{{c}^{2}}\int {T}_{\nu ,\mathrm{ff}}d{\rm{\Omega }},\end{eqnarray} \tag{ A3 }$$

where Ω is the solid angle.

The optical depth of the HRL line center under LTE condition (e.g., Wilson et al. 2013) is

$$\begin{eqnarray}&&{\tau }_{\mathrm{HRL},0}=1.92\times {10}^{3}{\left(\displaystyle \frac{{T}_{e}}{{\rm{K}}}\right)}^{-5/2}\left(\displaystyle \frac{\mathrm{EM}}{\mathrm{pc}\,{\mathrm{cm}}^{-6}}\right){\left(\displaystyle \frac{{\rm{\Delta }}\nu }{\mathrm{kHz}}\right)}^{-1}.\end{eqnarray} \tag{ A4 }$$

For H30α line (231.9 GHz), the line center optical depth (Equation (A4)) is

$$\begin{eqnarray}&&{\tau }_{{\rm{H}}30\alpha ,0}=2.48{\left(\displaystyle \frac{{T}_{e}}{{\rm{K}}}\right)}^{-5/2}\left(\displaystyle \frac{\mathrm{EM}}{\mathrm{pc}\,{\mathrm{cm}}^{-6}}\right){\left(\displaystyle \frac{{\rm{\Delta }}v}{\mathrm{km}\,{{\rm{s}}}^{-1}}\right)}^{-1},\end{eqnarray} \tag{ A5 }$$

where Δν and Δv are the FWHM of the line in the frequency and the velocity, respectively.

For a typical electron temperature of 10⁴ K and typical HRL line width of 40 km s⁻¹, H30α is optically thin except for extremely high emission measures (EM > 1.6 × 10¹¹ pc cm⁻⁶). Therefore, the brightness temperature of the free–free continuum subtracted from H30α is

$$\begin{eqnarray}&&{T}_{\mathrm{HRL}}={T}_{e}{e}^{-{\tau }_{\mathrm{ff}}}\left(1-{e}^{-{\tau }_{\mathrm{HRL}}}\right)={T}_{e}{e}^{-{\tau }_{\mathrm{ff}}}{\tau }_{\mathrm{HRL}}.\end{eqnarray} \tag{ A6 }$$

Note that here the factor of ${e}^{-{\tau }_{\mathrm{ff}}}$ includes the effects of the non-optically thin free–free continuum on the line intensity. From Equations (A5) and (A6), the integrated H30α intensity is then

$$\begin{eqnarray}&&\left(\displaystyle \frac{\int {T}_{{\rm{H}}30\alpha }{dv}}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)=2.64{\left(\displaystyle \frac{{T}_{e}}{{\rm{K}}}\right)}^{-3/2}\left(\displaystyle \frac{\mathrm{EM}}{\mathrm{pc}\,{\mathrm{cm}}^{-6}}\right){e}^{-{\tau }_{\mathrm{ff}}},\end{eqnarray} \tag{ A7 }$$

and the total H30α flux density integrated over velocity and space is

$$\begin{eqnarray}&&{S}_{{\rm{H}}30\alpha }=\displaystyle \frac{2k{\nu }^{2}}{{c}^{2}}\iint {T}_{{\rm{H}}30\alpha }{dvd}{\rm{\Omega }}.\end{eqnarray} \tag{ A8 }$$

Note that if the free–free and HRL emissions are optically thin, we have a line-to-continuum ratio of

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\int {T}_{\mathrm{HRL}}{dv}}{{T}_{\mathrm{ff}}}\right)}_{\mathrm{thin}}=\left(\displaystyle \frac{{\tau }_{\mathrm{HRL},0}}{{\tau }_{\mathrm{ff}}}\right)\,{\rm{\Delta }}v\,\left(\displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{\pi }{\mathrm{ln}2}}\right)\end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray}&&\ \ =\ 7.445\times {10}^{3}\,\mathrm{km}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{{T}_{e}}{{\rm{K}}}\right)}^{-1.15}{\left(\displaystyle \frac{\nu }{\mathrm{GHz}}\right)}^{1.1}.\end{eqnarray} \tag{ A10 }$$

For the H30α line, we have

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\int {T}_{{\rm{H}}30\alpha }{dv}}{{T}_{\mathrm{ff}}}\right)}_{\mathrm{thin}}=74.7\,\mathrm{km}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{{T}_{e}}{{10}^{4}\,{\rm{K}}}\right)}^{-1.15}.\end{eqnarray} \tag{ A11 }$$

The true line-to-continuum ratio should be lower than this if the free–free continuum is not optically thin, or higher than this if non-LTE effects become significant (as would be the case for an HRL maser).

For a uniform distribution of the emission measure, the free–free spectral index is α_ff ≈ 2 (optically thick) in lower ν, and α_ff ≈ −0.1 (optically thin) in higher ν, with a relatively narrow transition frequency range. The observed fluxes in VLA bands show consistent spectral indices between 0 and 2 for a quite wide frequency range, which requires a nonuniform distribution of the emission measure (and electron density). For simplicity, we assume the emitting region projected on the plane of sky is a circular disk with an outer radius of r_out, and the emission measure follows a power-law dependence on the radius r,

$$\begin{eqnarray}&&\mathrm{EM}(r)={\mathrm{EM}}_{10\mathrm{au}}{\left(\displaystyle \frac{r}{10\,\mathrm{au}}\right)}^{-{\rho }_{\mathrm{EM}}}.\end{eqnarray} \tag{ A12 }$$

We also assume the disk has an inner radius of r_in to avoid singularity. Then we have the free–free and H30α fluxes as

$$\begin{eqnarray}\begin{array}{rcl}\left(\displaystyle \frac{{S}_{\nu ,\mathrm{ff}}}{\mathrm{Jy}}\right) & = & 4.534\times {10}^{-10}{\left(\displaystyle \frac{\nu }{\mathrm{GHz}}\right)}^{2}{\left(\displaystyle \frac{d}{\mathrm{kpc}}\right)}^{-2}\\ & & \times \left(\displaystyle \frac{{T}_{e}}{{\rm{K}}}\right){\displaystyle \int }_{{R}_{\mathrm{in}}}^{{R}_{\mathrm{out}}}\left(1-{e}^{-{\tau }_{\mathrm{ff}}}\right){RdR},\end{array}\end{eqnarray} \tag{ A13 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}\left(\displaystyle \frac{{S}_{{\rm{H}}30\alpha }}{\mathrm{Jy}\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right) & = & 6.445\times {10}^{-5}{\left(\displaystyle \frac{d}{\mathrm{kpc}}\right)}^{-2}{\left(\displaystyle \frac{{T}_{e}}{{\rm{K}}}\right)}^{-3/2}\\ & & \times \left(\displaystyle \frac{{\mathrm{EM}}_{10\mathrm{au}}}{\mathrm{pc}\,{\mathrm{cm}}^{-6}}\right){\displaystyle \int }_{{R}_{\mathrm{in}}}^{{R}_{\mathrm{out}}}{e}^{-{\tau }_{\mathrm{ff}}}{R}^{1-{\rho }_{\mathrm{EM}}}{dR},\end{array}\end{eqnarray} \tag{ A14 }$$

where R = r/10 au, R_in = r_in/10 au, R_out = r_out/10 au, and τ_ff is given by Equation (A1). These are the observed quantities we fit with the model.

To further reduce the free parameters in the model, we adopt r_in = 10 R_⊙ (0.05 au), and r_out = 50 au. The inner radius of 10 R_⊙ is a typical stellar radius for massive protostars reaching the main sequence (e.g., Zhang & Tan 2018). The outer radius of 50 au is roughly consistent with the observation (see Figure 6). Therefore, we have three free parameters in this model, the electron temperature of the ionized gas (T_e ), the emission measure at the radius of 10 au (EM_10au), and the power-law index of the emission measure distribution ρ_EM. We expect the electron density distribution has a power-law form ${n}_{e}\sim {r}^{-{\rho }_{{n}_{e}}}$ with $1.5\leqslant {\rho }_{{n}_{e}}\leqslant 2$. The power-law index of −1.5 for the electron density corresponds to the case of the surface of a photoionized disk around a massive forming star (Tanaka et al. 2013). The power-law index of −2 of the electron density corresponds to the case of an expansion flow of the ionized gas with constant velocity. The emission measure is then $\mathrm{EM}\sim {n}_{e}^{2}(r)l$, where l is the line-of-sight length scale of the ionized region. As a first-order estimate, we assume l ∝ r and $\mathrm{EM}\sim {n}_{e}^{2}(r)r\sim {r}^{1-2{\rho }_{{n}_{e}}}$. Therefore, we limit the power-law index of the emission measure distribution to the range of 2 ≤ ρ_EM ≤ 3 in the fitting.

To find the best-fit model, we compare the model free–free continuum flux densities (Equation (A13)) to the continuum fluxes observed in VLA bands, and compare the model H30α flux integrated over velocity and space to the observed value from the ALMA observation. For each model, we calculate the test statistics ${\chi }_{\mathrm{ff}}^{2}$ and ${\chi }_{{\rm{H}}30\alpha }^{2}$ as

$$\begin{eqnarray}&&\begin{array}{l}{\chi }_{\mathrm{ff}}^{2}=\sum \displaystyle \frac{{({S}_{\nu ,\mathrm{ff},\mathrm{mod}}-{S}_{\nu ,\mathrm{ff},\mathrm{obs}})}^{2}}{{\sigma }_{\nu ,\mathrm{ff},\mathrm{obs}}^{2}},\\ {\chi }_{{\rm{H}}30\alpha }^{2}=\displaystyle \frac{{({S}_{{\rm{H}}30\alpha ,\mathrm{mod}}-{S}_{{\rm{H}}30\alpha ,\mathrm{obs}})}^{2}}{{\sigma }_{{\rm{H}}30\alpha ,\mathrm{obs}}^{2}},\end{array}\end{eqnarray} \tag{ A15 }$$

where the summation of ${\chi }_{\mathrm{ff}}^{2}$ is for all the VLA bands. Then, we determine the best-fit model by minimizing the value of ${\chi }^{2}=({\chi }_{\mathrm{ff}}^{2}+{\chi }_{{\rm{H}}30\alpha }^{2})/N$, where N is the total number of data points (including continuum and H30α line) used in fitting. Note that the ALMA 1.3 mm continuum is not used for fitting the free–free emission, as it usually includes significant dust emission.

If the H30α line has a non-LTE maser component (such as in source 3), we use the observed H30α flux as an upper limit. If the H30α line is not detected (such as in source 1), we use the 3σ level in the integrated map as an upper limit. If the observed H30α and ALMA 1.3 mm continuum fluxes indicate strong dust continuum contributions even in VLA bands, we also include the dust component following

$$\begin{eqnarray}&&{S}_{\nu ,\mathrm{dust}}={S}_{234\mathrm{GHz},\mathrm{dust}}{\left(\displaystyle \frac{\nu }{234\,\mathrm{GHz}}\right)}^{{\alpha }_{\mathrm{dust}}},\end{eqnarray} \tag{ A16 }$$

and fit all the continuum fluxes (VLA and ALMA) with S_ν = S_ν,dust + S_ν,ff.

From the derived properties, we further estimate the ionizing photon rate by balancing the rates of recombination and photoionization, following the method by Schmiedeke et al. (2016),

$$\begin{eqnarray}&&{\dot{N}}_{{\rm{i}}}=\int {n}_{e}^{2}(\beta -{\beta }_{1}){dV}={d}^{2}\int \mathrm{EM}(\beta -{\beta }_{1})d{\rm{\Omega }},\end{eqnarray} \tag{ A17 }$$

where β and β₁ are the rate coefficients for the recombination to all levels and to the ground state, i.e.,

$$\begin{eqnarray}&&\beta -{\beta }_{1}=4.1\times {10}^{-10}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{{T}_{e}}{{\rm{K}}}\right)}^{-0.8}.\end{eqnarray} \tag{ A18 }$$

From our assumptions on the emission measure distribution, we obtain the ionizing photon rate of

$$\begin{eqnarray}\begin{array}{rcl}{\dot{N}}_{{\rm{i}}} & = & 1.1\times {10}^{45}\,{{\rm{s}}}^{-1}{\left(\displaystyle \frac{{T}_{e}}{{10}^{4}\,{\rm{K}}}\right)}^{-0.8}\\ & & \times \left(\displaystyle \frac{{\mathrm{EM}}_{10\mathrm{au}}}{{10}^{10}\,\mathrm{pc}\,{\mathrm{cm}}^{-6}}\right){\displaystyle \int }_{{R}_{\mathrm{in}}}^{{R}_{\mathrm{out}}}{R}^{1-{\rho }_{\mathrm{EM}}}{dR}.\end{array}\end{eqnarray} \tag{ A19 }$$