© The Authors 2025

1 Introduction

Recent developments in instrumentation and, more specifically, the advent of the Atacama Large Millimeter/submillimeter Array (ALMA), have made it possible to perform population studies of disks around low-mass young stellar objects (YSOs) and T Tauri stars. Studying such disk populations in a statistical manner and linking disk properties to age, exoplanet demographics, and stellar parameters sheds light on the evolution of disks, including their dissipation timescales and mechanisms, as well as (sub)stellar companion formation within them (e.g., Pascucci et al. 2016; Ansdell et al. 2017; Andrews et al. 2013; Tripathi et al. 2017; Pinilla et al. 2020; Manara et al. 2023). Recently, this endeavor has been extended toward higher masses by Stapper et al. (2022), who incorporate Herbig AeBe stars (mostly) in the lower mass range around 2 M_⊙, with afew higher- mass objects up to 10 M_⊙. They conclude that the disks of the studied Herbig stars are skewed toward higher masses and larger sizes than those of T Tauri stars.

It remains challenging to push this endeavor to higher masses as the number of young objects becomes progressively smaller with increasing mass, and the birth environments of these objects are complicated by large distances, crowding, and extinction. As a result, many of the questions regarding the lifetimes and dissipation mechanisms of disks, and the mechanisms of stellar (and possibly planetary) companion formation remain elusive in the higher mass ranges. A fair number of disks, outflows, and disk-like structures have been observed at millimeter and centimeter wavelengths around massive (≥8 M_⊙) young or protostellar objects (e.g., Johnston et al. 2015, 2020; Beltrán & de Wit 2016; Ilee et al. 2016; Caratti o Garatti et al. 2017; Cesaroni et al. 2017; Ilee et al. 2018; Beuther et al. 2018; Maud et al. 2018; Sanna et al. 2019; Burns et al. 2023). However, the observed structures extend to scales from hundreds to thousands of au and, contrary to the disks in the quoted disk population studies, belong to the embedded early stages of formation when the photosphere of the central object is still obscured and its (final) mass is uncertain (but see McLeod et al. 2024).

We searched for candidate massive YSOs in star-forming regions (SFRs) associated with ultra-compact H II regions (e.g., Bik & Thi 2004; Hanson et al. 1996; Lumsden et al. 2013). In this context, a region of particular interest is the M17 H II region (Fig. 1) in the east side of the Omega (M17) giant molecular cloud complex, at a distance of ∼1.7 kpc (Maíz Apellániz et al. 2022; Stoop et al. 2024). Its ionizing source, the young massive cluster NGC 6618, contains more than 100 B stars and at least 16 O stars. At least four of these stars have spectral types earlier than O6 V, the spectral type of θ¹ Ori C, the most luminous star in the Trapezium region of the Orion Nebula Cluster (ONC) (Hanson et al. 1997; Povich et al. 2007, 2009; Hoffmeister et al. 2008). With an estimated total stellar mass of ∼8000–10 000 M_⊙, NGC 6618 is about five to six times as massive as the ONC (Povich et al. 2009; Kroupa et al. 2018). An important property of the candidate massive YSOs in this region is that they have visible photospheres, which allows placing them on the Hertzsprung-Russell diagram (HRD; Ochsendorf et al. 2011).

Ramírez-Tannus et al. (2017, hereafter RT17) characterized the young stellar population in NGC 6618 constraining its age to ⪅1 Myr. Recently, the distance estimate to M17 was adjusted from 1.98 kpc (Xu et al. 2011) to ∼1.7 kpc using Gaia DR3 data (Maíz Apellániz et al. 2022; Stoop et al. 2024). The age estimate was also refined by Stoop et al. (2024) who used the kinematic age of dynamically ejected massive runaway stars yielding an age of 0.65 ± 0.25 Myr. Within NGC 6618, RT17 identified a sample of pre-main-sequence (PMS) stars with masses ∼4–15 M_⊙ that optically reveal their photospheres and, at the same time, show optical and near-infrared (NIR) signatures of a circumstellar disk extending close to the stellar surface (Backs et al. 2023; Poorta et al. 2023). Their stellar parameters were first determined by RT17 and later updated by Backs et al. (2024; in prep.), by fitting the normalized photospheric spectra using the non-LTE stellar atmosphere code FASTWIND (Santolaya-Rey et al. 1997; Puls et al. 2005) and a genetic algorithm based fitting approach KIWI - GA (Brands et al. 2022). In addition to the stellar properties, extensive modeling of the circumstellar disk features in the optical and NIR (Sect. 2.1) suggest that these objects are in the final stages of formation and have small-scale gaseous and dusty disks. However, constraints on the existence of extended disk material and its mass are lacking so far.

The sample of PMS stars in M17 fits the common definition of Herbig Be stars in the literature, but the objects have not been included in statistical studies and catalogs of Herbig stars (e.g., Vioque et al. 2018; Brittain et al. 2023), likely due to their distance, extinction, and the fact that they reside in the complicated environment of a massive SFR. Additionally, Stapper et al. (2022) study a volume-limited sample out to 450 pc, excluding any targets in M17.

In the context of the dearth of disk studies in the higher stellar mass regime, and in massive SFRs in general, this study aims to take the next step by constraining the extended disk properties of four intermediate- to high-mass PMS stars in M17 with masses ∼4–10 M_⊙. With ALMA we detected, for the first time, the millimeter components of these YSOs, adding to the handful of similar objects for which ALMA observations have been obtained, such as MWC 297 (B1.5 Ve, ∼10 M_⊙), which is a rapidly rotating massive YSO at a distance of 418 pc (Sandell & Vacca 2023; Stapper et al. 2022), and IRS2 (early B) in NGC2024 at a distance of 450 pc (Lenorzer et al. 2004; van Terwisga et al. 2020). We adopted similar and refined techniques to constrain the disk mass, radius, and structure of our targets. For the first time we applied a multiwavelength approach to this kind of object, including both photospheric optical to NIR spectra, as well as optical to mid-infrared (MIR) photometry and the newly obtained ALMA millimeter wavelength data.

In the next section (Sect. 2), we present the M17 targets and the ALMA observations. All four targets are detected, and in the covered fields we identify four additional sources, most likely low-mass YSOs. In Sect. 3, we analyze the ALMA data and introduce the methods that we used to measure the disk masses. The millimeter emission is thermal in nature, that is, free-free radiation and/or thermal dust emission. Different estimates of the disk masses are presented in Sect. 4. We discuss our results in Sect. 5 and compare them with those obtained for T Tauri and Herbig stars. In the final section (Sect. 6) we summarize our conclusions.

Fig. 1 Overview of the M17 region in different wavelengths. The left panel shows a NIR color composite image of M17 based on 2MASS data: J (blue), H (green), and K (red), with large-scale 1.3 mm continuum emission contours (archival ALMA Band 6 data, project code 2018.1.01091.S). Contour levels: 5, 15, 30, and 70 × noise rms (respectively 0.018, 0.053, 0.11, and 0.25 Jy). The middle panel shows only the ALMA mosaic, with the beam size in the bottom left corner. In both panels the rectangular region indicates the field in which all four targets in this study (indicated in red) are found. The rightmost panel is a zoomed-in image of this region with all targets and serendipitous discoveries marked in black. To give an impression of the scales, the radius of each object dot is 1″, which corresponds to five times the image size of each panel in Fig. D.1.

2 ALMA observations of M17 targets

2.1 M17 target sample

Four of the M17 PMS stars studied by RT17 were observed with ALMA. These four targets are listed in Table 1 along with their stellar parameters (Backs et al. 2024). Several recent studies were performed to characterize the inner disks around (some of) these objects. Backs et al. (2023) studied the double-peaked hydrogen emission lines in VLT/X-shooter spectra of B243 and B331 using the radiation thermo-chemical code PRODIMO (e.g., Woitke et al. 2016). They find that these lines form very close to the star, such that the disks likely continue (almost) up to the stellar surface, indicating the expected lack of magnetic fields that are typical for low-mass stars (e.g., Alecian et al. 2013). Based on the line-forming region, they estimate inner disk (<20 au) masses on the order of 10⁻³–10⁻⁴ M_⊙. They also find that line luminosity-accretion relations derived for low- to intermediatemass stars (Fairlamb et al. 2017) imply accretion rates of 10⁻⁵ to 10⁻³ M_⊙ yr⁻¹. This would entail a substantial mass reservoir farther out in the disk or in-falling onto the inner disk, which, given the evolutionary state of the stars, seems unlikely.

Derkink et al. (2024) looked at spectroscopic and photometric variability over days to years in multi-epoch data of B268, B275, and B243. They find intrinsic spectroscopic variability for all these objects in the atomic lines originating mostly close to the stellar surface, as well as evidence for asymmetric disk structures in the inner disk, such as spiral arms. They detect no strong photometric variations, (inverse) P Cygni line profiles, or high-velocity components in forbidden emission lines, making the presence of strong accretion (bursts) or jets unlikely.

Finally, Poorta et al. (2023) modeled the first and second overtone CO bandhead emission in all four targets to constrain densities and temperatures in the inner disk. In line with previous studies, they find that this emission is typical for high column densities (N_CO ∼ 10²¹ cm⁻²) and temperatures (2000–5000 K), and originates close to the star (<1 au), though not as close as the previously mentioned hydrogen emission lines. They also fit near- to mid-infrared photometry with a dust disk model, and find that for a canonical dust-to-gas ratio the derived column densities in the inner disk are inconsistent with those derived from the CO emission (see also Sect. 3.5.1). The inner disk structure, though critical for constraining the mass-accretion rate, provides little information about global properties of the disk, such as total mass and size. To investigate these properties, longer wavelength data are needed.

In Fig. 1, we present an overview of the M17 region in 2MASS colors and large-scale millimeter emission in ALMA Band 6 from archival data (project code 2018.1.01091.S). The targets in this study are marked in each panel, and the right panel shows a close-up of the region containing all four targets, as well as the new detections with ALMA that were made in the field of view around three of the four targets (see Sect. 2.3). Two new detections (B243 SW and B275 SE) clearly have 2MASS counterparts. The figure illustrates the high source density as well as the large-scale emission in the region. This emission is largely resolved out in the high-resolution observations taken for this study (Table 2). Figure D.1 shows a zoomed-in image of all detected sources (see also Sect. 2.3).

2.2 ALMA observations and imaging

Under project ID 2019.1.00910.S we obtained ALMA 12m array observations in Band 6 and Band 7 for our four targets in M17, resulting in eight fields to be imaged. The observational setups for the targets were close to identical; an overview of the observations per band is given in Table 2. Specific details per target are provided in Table A.1. The observations were designed to detect the targets in the continuum, but also included the ¹²CO, ¹³CO, and C¹⁸O rotational 2–1 lines in Band 6 and ¹²CO rotational 3–2 line in Band 7. The respective spectral setups are summarized in Tables 3 and 4. The Band 7 data were manually calibrated by the ALMA staff using CASA version 6.1.1 (CASA Team 2022). The Band 6 data were calibrated with the ALMA pipeline version 2021.2.0.128 using CASA version 6.2.1.7. The signal-to-noise ratio of the detected targets was too low for selfcalibration. For all subsequent image deconvolution (imaging) and flux determination we used CASA version 6.2.1.7.

We imaged the continuum using natural weighting to maximize sensitivity, averaging over all spectral windows after including only line-free regions. Because the observed targets lie in a very young SFR, in some cases we had to exclude the data from shorter baselines that are sensitive to extended emission in order to improve the signal-to-noise ratio of the observed targets. In Band 6 we used uv-ranges of >50 kλ, thereby excluding emission from scales larger than ∼4.13″, which at the distance of M17 corresponds to ∼7000 au. In Band 7 the uv-range was already limited to >50 kλ when all baselines were included, and we applied no further limit. An exception was made in both bands for the images surrounding B331. As can be seen in Fig. 1, this target is located in a region with continuum emission on ≳0.4″ (700 au) scales, possibly associated with remnant cloud emission or a chance alignment, and we removed baselines <500 kλ for images in both bands, thereby excluding emission from scales larger than ∼700 au. Zoomed-in images around each detected source (see Sect. 2.3) are shown in Fig. D.1. Larger field images showing the background and multiple sources at once are presented online at Zenodo.

In an attempt to detect the targets in the selected CO lines we inspected image cubes averaged to 1–2 km s⁻¹. We found that all line emission was dominated by large scales, likely from the surrounding cloud. Removing shorter baselines did not yield detections of the continuum-detected sources in ¹²CO or in ¹³CO. In order to estimate the upper limits on the gas mass of possible disks (see Sect. 3.2) we determined the root mean square (rms) of the noise in the line images with uv-ranges of >500 kλ (removing emission from scales larger than ∼700 au). Assuming a full width at half maximum (FWHM) of 10 km s⁻¹, we measured the noise rms in a cube with 10 km s⁻¹ wide channels to be ∼0.8 mJy beam⁻¹, translating to a 3σ upper limit to the integrated line flux of ∼24.0 mJy beam⁻¹ km s⁻¹ for the ¹²CO (3 = 0; J = 2–1) line on scales ≲700 au.

2.3 Detections and determination of fluxes in the continuum

The noise rms was determined from a large emission-free region (radius ∼5″) in the images resulting from the procedure described in Sect. 2.2. Similar different regions were probed in each image and the variation of the rms between them was found to be negligible. A source was considered a detection if its flux was higher than five times the noise rms in one or both bands. Some sources minimally fulfilling this condition on the edge of the primary beam were considered noise peaks after visual inspection. Detections that are not among the original targets are named after the target in whose image they are detected; the upper case letters indicate the position of the detection relative to the original target (e.g., SW for southwest). Three of the four original targets are detected in both bands and four new sources are detected in the fields surrounding the original targets (see Fig. D.1 and the field images online). B268 is only possibly detected in Band 6 (S/N 3), but is included in Fig. D.1 and in the further analysis because it is an original target. This object was imaged with a uv-taper of 0.08″ to increase the signal-to- noise ratio, leading to the large beam size relative to the other images. B331 NE is detected in Band 7 (S/N 6) according to the set criteria; however, it is located in the high-noise region on the edge of the primary beam, is not detected in Band 6, and has no infrared counterpart (see NIR image and ALMA field image). For completeness we present it in Fig. D.1 and list it in Tables 5 and 6, but otherwise exclude it from further analysis. All other detections are in both bands and have optical and/or NIR counterparts. Their exact positions and projected separations are listed in Table 5.

Fluxes were determined using the imfit task in CASA, which fits a two-dimensional Gaussian to the selected image component (full fit results are provided in Appendix B). All resulting major-axis FWHM were at most on the order of the beam size, so we concluded that none of the detected sources is resolved. As expected in that case, the resulting peak fluxes (in units of Jy/beam) are all equal within the errors to the integrated fluxes (in units of Jy). Since the sources are contained in one beam we report only the peak fluxes (in units of Jy). The errors on the flux are taken to be the maximum of the standard flux calibration error of 10%, the noise rms and statistical errors from the fit. The resulting fluxes and their errors are listed in Table 6.

In addition to the ALMA fluxes, additional photometry and Very Large Array (VLA) radio fluxes were collected from catalogs and tables in the literature for all detected sources. For the original targets the found photometric points are listed in Poorta et al. (2023). For the new detections the photometric flux points are listed in Table C.1. VLA fluxes were taken from two extensive surveys in the M17 region with sensitivity down to 0.005 mJy at a beam size of about 0.04 arcsec² (Rodríguez et al. 2012; Yanza et al. 2022). The only detected source is B331, in both surveys. We list a 5σ upper limit for all other ALMA sources. The values are given in Table 6. Finally, all recovered flux points are plotted in Fig. 3.

All the new sources are also detected in acquisition images in the K-band (central wavelength λ_C ∼ 2.2 µm and FWHM ∼0.41 µm) obtained with the Large Binocular Telescope (LBT) on 5 and 6 July 2023 (A. Derkink, 2023, priv. comm.). In these images (see online material), the sources are clearly distinct from the original targets and have not moved since the ALMA observations almost two years earlier (Table 2).

3 Analysis

In order to draw conclusions on the measured fluxes and to constrain the source characteristics, we performed different analyses. First, to determine whether we can attribute the measured continuum fluxes to dust emission from a circumstellar disk or to free-free emission from ionized material we determined the spectral index (Sect. 3.1); we discuss its interpretation in Sect. 4.1. Second, under simplifying assumptions we calculated an upper limit on the gas mass from the nondetection of CO lines (Sect. 3.2). Third, we estimated dust masses from the measured fluxes directly (Sect. 3.3). We also compared the measured fluxes with those obtained in a pre-calculated PRODIMO model grid (Sect. 3.4). Finally, we estimated the disk masses using an adapted version of the disk model developed by Poorta et al. (2023), that was used to fit CO bandhead emission from the same objects (Sect. 3.5).

3.1 Spectral index and origin of detected fluxes

The slope between two or more different radio regimes is a widely used diagnostic for characterizing the source of emission in the Rayleigh-Jeans limit. This slope is quantified as the spectral index which, following Beckwith et al. (1990) and Tychoniec et al. (2020), among others, is defined as $$\alpha =\frac{\log \left(F_{v_1}/F_{v_2}\right)}{\log \left(v_1/v_2\right)},$$(1)

with $F_{v_{1/2}}$ the flux at frequencies ν_1/2, and ν₁ > ν₂ such that a positive spectral index indicates a descending slope toward longer wavelengths. For all detected sources we calculated the spectral index between ALMA Bands 7 and 6. For B331 we also determined the index between ALMA Band 6 and the VLA X- band at 10 GHz. The adopted frequencies for these bands are ν₇ = 341.4183 GHz, ν₆ = 225.4311 GHz, and ν_X = 10 GHz for Band 7, Band 6, and VLA 10 GHz, respectively. For B331 we additionally report a spectral index between the two VLA bands at 4.96 GHz and 8.46 GHz from Rodríguez et al. (2012). All the resulting indices are listed in Table 6, and are discussed in Sect. 4.1.

3.2 Upper limit on the gas mass

Based on the nondetection of the rotational ¹²CO (J = 2–1) line, we estimated an upper limit on the gas mass in the possible disks. Given that the beam sizes and rms values are very similar for the different object images we did this for one average limiting flux value ∼24.0 mJy beam⁻¹ km s⁻¹ (Sect. 2.2). With a beam size of 0.07″ × 0.04″ this limit translates to an antenna temperature of ∼220 K km s⁻¹. We used the non-LTE molecular radiative transfer code RADEX¹ (van der Tak et al. 2007) to estimate line radiation temperatures, using a kinetic temperature T_kin = 100 K, an H₂ density of 10¹⁰ cm⁻³ leading to LTE conditions, and a line width of 10 km s⁻¹ based on Keplerian velocities at a distance of ∼30–50 au from our sources. Under these assumptions RADEX finds our limiting temperature with a CO column density N_CO ≈ 4 × 10¹⁷ cm⁻². The emission at this column density is optically thin (τ ≈ 0.31). Integrating over the beam size and assuming a CO abundance of n(CO)/n(H₂) = 10⁻⁴, this column density results in an upper limit on the gas mass of∼2× 10⁻⁵ M_⊙.

If the disk is much smaller than the 120 × 70 au beam, the emission is likely optically thick and, for an assumed kinetic temperature of 100 K, would have a brightness temperature of ∼100 K. Such a disk would need to take up less than 20% of the beam area to remain undetected, corresponding to a disk radius of <20–30 au, assuming an inclination in the range of 0–60^◦. Under the assumption of such a compact optically thick disk, the CO column density could in principle be arbitrarily high. However, in the studied objects a constraint is provided by the CO bandhead emission, from which Poorta et al. (2023) derive a CO column density at the inner rim of the disk of N_CO ∼ 1– 8 × 10²¹ cm⁻². This density is likely to decline farther out in the disk with typical power law exponents between −1 (e.g., Woitke et al. 2016) and −1.5 (e.g., Poorta et al. 2023). Hence, to derive an upper limit on the gas mass under these conditions we adopted a constant average CO column density of N_CO ∼ 5 × 10²⁰ cm⁻² and a CO abundance of n(CO)/n(H₂) = 10⁻⁴ throughout a disk of 30 au. In this case, we derive an upper limit on the gas mass of∼6× 10⁻³ M_⊙.

3.3 Dust mass estimation: Optically thin approach

Many studies use a simple approximation to estimate dust masses for populations of protoplanetary disks, assuming optically thin dust emission from an isothermal disk (e.g., Ansdell et al. 2016; Cazzoletti et al. 2019; Stapper et al. 2022). Under this approximation the dust mass in the disk is given by (Beckwith et al. 1990) $$M_{\text{dust }}=\frac{F_v{d}^2}{{\kappa }_v{B}_v\left(T_{\text{dust }}\right)},$$(2)

with F_ν the measured flux at frequency ν, d the distance to the source, κ_ν the dust opacity at ν, and B_ν(T_dust) the Planck function at the assumed dust temperature T_dust. For the dust opacity we adopted the power law $${\kappa }_v={\kappa }_0{\left(v/v_0\right)}^{\beta },$$(3)

with κ₀ = 10 cm² /g, ν₀ = 1000 GHz, and β = 1 from Beckwith et al. (1990). The commonly assumed dust temperature for disks around low-mass stars is T_dust = 20 K (Andrews & Williams 2005). However, we anticipated higher temperatures due to the higher luminosities of our stars. Following Andrews et al. (2013) and Stapper et al. (2022) we used $$T_{\text{dust}}=25{\left(\frac{L_{*}}{L_{\odot }}\right)}^{1/4}[\text{K}],$$(4)

resulting in dust temperatures between 110 K and 246 K for the objects in Table 1. For the serendipitous detections, for which we did not know the stellar luminosities, we adopted T_dust = 150 K.

We list the resulting dust mass for all detected sources in Table 6; however, we note that dust emission is not the most likely source of the measured flux for all of these objects (see Sect. 4.1). In Table 8, we list the dust masses under this approximation again, only for the original targets, with the purpose of comparing the results with the other estimates of disk gas and/or dust mass (see Sect. 4.2).

3.4 Disk mass estimation: PRODIMO model grid

To guide observations and to serve as a reference for a first analysis, we computed a grid of models with a version of the radiation thermo-chemical code PRODIMO (Woitke et al. 2009; Kamp et al. 2017) adapted for disks around hot stars (Backs et al. 2023). In these models we used the DIANA-standard model setup with a gas-to-dust mass ratio of 100 (Woitke et al. 2016). The grid consists of 10 × 10 × 5 models with, respectively, total disk masses ranging from 5 × 10⁻⁴ to 5 × 10⁻¹ M_⊙, outer radii from 20 au to 500 au, and inclinations from 15^◦ to 75^◦. The central object for this grid was chosen to be B331 (see Table 1). The inner radius was set to 11.2 au, the dust sublimation radius for this object.

In our current analysis, we used the resulting spectral energy distributions (SEDs) from these models solely as a reference point to compare with our other mass estimations (Table 8) and with the SEDs we calculate in Sect. 3.5. We fitted the model SEDs to the available photometry and ALMA points for each of the four targets in the sample. Before fitting we changed the underlying stellar continuum of each model to match the object in question. The best fit PRODIMO models are all the same, with parameters M_tot = 5 × 10⁻⁴ M_⊙ (the lowest in grid), outer radius 233 au², and inclination 75^◦ (the highest in grid). Because the grid models are unable to fit the near- to mid-infrared photometry simultaneously with the ALMA points, and also provide little constraint on the disk mass, we used another disk model to acquire an improved fit of the data, as described in the next section.

3.5 Dust and gas mass estimation: Thin disk modeling

3.5.1 Model description

It is notoriously difficult to model the inner and outer regions of disks in one self-consistent model. As mentioned in Sect. 2.1, several previous studies constrained the inner disk characteristics of the target sample based on infrared photometry and spectroscopy. In particular, Poorta et al. (2023, hereafter P23) constrained the column density and temperature of the gas in the inner parts of the disk (≲1–2 au) by fitting CO bandhead emission, and characterized the dust emission originating farther out (≲10–30 au) by fitting the near- to mid-infrared SED. For this work, our aim was to fit the SEDs including the infrared photometry and the new ALMA-points, while taking into account the previously obtained constraints on the inner disk. For this we adopted the P23 model, which is a 1D Keplerian disk model that uses power laws for the radial density and temperature structure. It predicts the emission components from the CO gas and the dust assuming constant conversion factors for the H₂ to CO abundance and the dust-to-gas mass ratio. For this work we modeled only the dust emission, with a few adaptations with respect to the original model. For the inner disk column densities we adopted the values resulting from fitting the first and second overtone CO bandheads.

Changes in the dust modeling pertain to the dust opacities and the dust-to-gas mass ratio. Where P23 use single grain-size astronomical silicate opacities from Laor & Draine (1993), we adopted the mass absorption coefficients calculated by Ossenkopf & Henning (1994) for dust in protostellar cores with the Mathis et al. (1977, i.e., MRN) size distribution. The tabulated dust opacities are for varying gas densities (n_H = 10⁵ to 10⁸ cm⁻³) and for no ice mantles, thin ice mantles, or thick ice mantles on the grains. Due to the relatively hot and dense environments expected for the studied objects, we adopted n_H = 10⁸ cm⁻³ and grains without ice mantles.

We also introduced a second change. Instead of using a constant dust-to-gas mass ratio M_d/M_g throughout the disk, we varied this quantity using the following logistic function (see Fig. 2 for an example): $$M_{\text{d}}/M_{\text{g}}={\left(M_{\text{d}}/M_{\text{g}}\right)}_{\text{base }}*{\left[e^{\beta *\left(r-r_{\text{turn }}\right)}\right]}^{-1}+{\left(M_{\text{d}}/M_{\text{g}}\right)}_{\text{lim }}.$$(5)

The function monotonically rises from the minimum limit (M_d/M_g)_base to the maximum limit (M_d/M_g)_lim as a function of disk radius r and free parameters β, the growth rate, and r_turn, the inflection point. As can be seen in Fig. 2, the function approaches a step function for high values of β and a linear function for low values.

The motivation for introducing this function is that the P23 disk models with a constant dust-to-gas ratio cannot account for all the data points at the same time. This also holds for the previously described PRODIMO models, which also have a constant dust-to-gas ratio. In both cases, models that fit the longer wavelength points strongly overpredict the infrared fluxes. This is likely connected to the fact that thermal dust emission in the near- to mid-infrared is not a good probe of the dust mass in the inner disk (e.g., P23, van der Marel 2023). This may be due to settling of the bulk of the dust to the mid-plane in the high-density inner regions of the disk, and only a small amount of grains in the disk surface layer contributing to the observable near- to mid-infrared thermal emission. This also leads to a sensitivity of the NIR emission to geometric effects such as the viewing angle. Since the P23 model lacks vertical structure, these effects would not be accounted for. Alternatively, the dust-to-gas ratio could intrinsically vary, due to dust depletion in the inner disk in the presence of, for example, gaps, traps, pebble formation, or other disk structures and processes (e.g., van der Marel 2023). It is conceivable that a combination of these mechanisms is at play. Regardless of the physical situation, the approach described here is a way to allow the dust distribution in the disk to deviate from the gas column density structure. We briefly discuss the interpretations and implications for the total disk mass in Sect. 5.1.1.

Fig. 2 Illustration of the logistic function (Eq. (5)) used for varying the dust-to-gas mass ratio Md/Mg in the disk. The limiting values are fixed to the same values as in the models used for fitting: (Md/Mg)base = 10−8 and (Md/Mg)lim = 10−2. The value of the free parameter β is varied here to illustrate its effect.

3.5.2 Extrapolation and fitting

With the P23 model and the previous modeling results we estimated the gas and the dust mass of the disk in two ways. First, not taking into account the ALMA flux points and without fitting, we directly extrapolated the column density in the inner disk derived from the CO bandhead fitting by P23, (Σ_i)_CO, using the power law slope from the same fits (q = −1.5). The gas mass is then given by $$M_{\text{gas}}=\frac{n\left({\text{H}}_2\right)}{n(\text{CO})}{m}_{\text{H}}\cdot {\left({\text{Σ}}_i\right)}_{\text{CO}}{\displaystyle {\int }_{{\left(R_i\right)}_{\text{gas}}}^{R_{\text{out }}}2}\pi r{\left(\frac{{\left(R_i\right)}_{\text{gas}}}{r}\right)}^{q}dr.$$(6)

The inner radius for the gas, (R_i)_gas, was again taken from the CO bandhead fits (P23). The outer radius, R_out, was set to half the beam size of the image for each object (Tables 7, B.1, and B.2). Though B268 was imaged with a uv-taper to increase the signal- to-noise ratio (see Sect. 2.3 and Fig. D.1), for the purposes of this analysis we adopted a similar constraint on R_out as for the other sources. Similarly, the dust mass is given by $$M_{\text{dust}}=\frac{M_{\text{d}}}{M_{\text{g}}}\cdot \frac{n\left({\text{H}}_2\right)}{n(\text{CO})}{m}_{\text{H}}\cdot {\left({\text{Σ}}_i\right)}_{\text{CO}}{\displaystyle {\int }_{{\left(R_i\right)}_{\text{dust}}}^{R_{\text{out}}}2}\pi r{\left(\frac{{\left(R_i\right)}_{\text{gas}}}{r}\right)}^{q}dr.$$(7)

The inner radius for the dust, (R_i)_dust, was defined as the point where the dust temperature reaches the dust sublimation temperature 1500 K and was estimated with a radiative equilibrium approach (following Lamers & Cassinelli 1999), using the dust opacities, the stellar temperature, and stellar radius (Table 1). Using the described logistic function (Eq. (5)) for the dust-to-gas mass ratio has a negligible effect on the total disk mass. Therefore, for the purpose of this extrapolation we assumed a constant dust-to-gas mass ratio M_d/M_g = 10⁻².

Second, we fitted an SED to the photometry and ALMA flux points using the P23 model with the described adaptations. The SED has a stellar component that is represented by a Kurucz model matching the stellar parameters (Table 1). The dust emission from the disk is what we focused on. The fixed parameters going into the modeling were the inner and outer dust disk radii (R_i)_dust and R_out (defined as before), the inclination i (taken from P23 CO fit results), and a temperature power law exponent, which we fixed to p = −1 based on exploratory fitting. We fitted the initial column density Σ_i (at (R_i)_dust), power law exponent q, and the logistic function parameters r_turn and β (Eq. (5)). The nonlinear least-squares routine curve_fit from the SciPy module (version 1.9.3; Virtanen et al. 2020) was used to fit the data. The inverse of the errors on the flux were treated as weights, and the errors on the fit parameters were calculated from the covariance matrix returned by the routine. The initial values and bounds of the fit parameters, as well as the fixed parameters and the results are listed in Table 7. The initial value of Σ_i and its bounds were determined by extrapolating (Σ_i)_CO (P23) and its errors to (R_i)_dust, using q = −1.5 (P23). The resulting best fit SEDs are presented together with the data in Fig. 3. In Fig. 4 we show the temperature and dust-to-gas ratio of the best fit models along with their cumulative flux at two different wavelengths, all as a function of the disk radius.

The gas and disk masses resulting from the two approaches described in this section (extrapolation and fitting) are listed in Table 8. These figures and tables are discussed in Sect. 4.2.

4 Results

4.1 Spectral indices and nature of the sources

The most important diagnostics used to distinguish between the different physical origins of radio emission are source size and spectral index. We first briefly review the most common emission sources and their typical scales and spectral indices, and narrow down the possible mechanisms for the detected sources. We then consider the spectral indices for the sources separately.

4.1.1 Origin of radio emission and relation to spectral index

In general, the physical origin of radio emission can be thermal or nonthermal in nature (Yanza et al. 2022). Synchrotron emission by charged particles accelerated in a magnetic field is nonthermal. In a young massive SFR such as M17 this could originate from the hot coronas of young low-mass stars (e.g., Feigelson & Montmerle 1999) or from colliding winds in massive binary systems (e.g., Dougherty & Williams 2000). On much larger spatial scales than considered here, super-shells, for example due to collective stellar winds, may also produce nonthermal emission (e.g., Lozinskaya 1999). Nonthermal emission is most commonly, though not exclusively, associated with a negative spectral index, which we do not observe for any of the sources. For the original targets both mechanisms of compact nonthermal emission are very unlikely. On the one hand, the sources are too massive to harbor a magnetic field associated with a convective envelope. On the other hand, they show no signs of (near) equal mass binarity or strong stellar winds (Ramírez-Tannus et al. 2017; Backs et al. 2024). Furthermore, the only object massive enough to generate a significant stellar wind would be B331 (∼10 M_⊙), but an unresolved companion of similar mass also harboring such a wind would have resulted in a more luminous object. For the new detections a contribution from nonthermal emission cannot be excluded, as these sources could well be low-mass YSOs, but given their spectral indices (Sect. 4.1.2), it would unlikely be the only contributing mechanism.

Thermal radio emission can be due either to free-free radiation or to thermal dust emission. Since free-free radiation is caused by the (negative) acceleration of charged particles, it is always associated with ionized media. The most common of these that are relevant to the observed region are (remnants of) ultra- or hyper-compact H II-regions surrounding young OB stars (e.g., Avalos et al. 2009; Sánchez-Monge et al. 2011), ionized winds of massive stars (e.g., Snell & Bally 1986), externally ionized protoplanetary disks (proplyds) (e.g., Stecklum et al. 1998; Zapata et al. 2004), jets from YSOs (e.g., Anglada et al. 2018), and photo-evaporating disk winds (e.g., Lugo et al. 2004). Thermal dust emission can originate either from a circumstellar disk or from the surrounding molecular cloud. In the latter case, the emission is more diffuse.

The main question for the original target sample is whether the detected emission can be solely or partially attributed to thermal dust emission from a circumstellar disk, or whether free-free emission should be considered. Since all the detected sources are compact (≲150 au in diameter), a large-scale origin of the emission, such as diffuse free-free or dust emission from the surrounding H II-region and molecular cloud, as well as large- scale jets and outflows can be readily excluded as sources of the measured fluxes. Possible contributors to free-free emission originating close to the star include a small-scale H II-region³, the base of a radio-jet, or ionized gas in a stellar wind or a photo- evaporative disk wind. The last are commonly associated with a spectral index α ≃ 0.6 (Panagia & Felli 1975; Wright & Barlow 1975; Lugo et al. 2004). If all the ionized gas close to the star is optically thin, the emission will have an index α ≃ −0.1 in the observed ALMA frequencies.

Finally, the spectral index for optically thick thermal dust emission is that of a blackbody, that is α = 2. For optically thin dust the spectral index will depend on the grain size distribution. Small grains in the interstellar medium (ISM) result in ${\alpha }_{\text{mm}}^{\text{ISM}}=3.7$, whereas lower indices, between 2 and 3, are interpreted as a sign of grain growth in protoplanetary disks around low-mass stars (van der Marel 2023). If the index in such disks is higher, it may indicate a lack of larger grains.

Fig. 3 Best fit SEDs, NIR and MIR photometry points, and ALMA and VLA flux points for our four original targets and four sources serendipitously detected in the covered fields. Each panel shows the flux data for the original target and the Kurucz model for its stellar continuum in blue. The best fit model SED from the thin disk modeling (described in Sect. 3.5) is plotted in orange. Each panel also shows the flux data for the serendipitous detection(s) made in the field of the respective original target (in green and red); these data are not fitted. The downward arrows indicate the upper limit derived for nondetections (mostly VLA 10 GHz). Only B331 (upper right panel) is detected at VLA wavelengths. The dotted lines indicate the slopes between the radio flux points for this object, with a spectral index ∼0.67 for the VLA points at 4.96 GHz and 8.46 GHz (Rodríguez et al. 2012), and ∼0.52 between the ALMA and VLA 10 GHz flux points (Yanza et al. 2022). There is an offset between the data from Rodríguez et al. (2012) and (Yanza et al. 2022).

Fig. 4 Dust temperature (blue) and dust-to-gas mass ratio (black) as a function of disk radius resulting from the thin disk model fit (described in Sect. 3.5) for each object. The dotted black vertical line indicates the fitted inflection point Rturn of the logistic function of the dust-to-gas mass ratio. In red is the normalized cumulative flux distribution to illustrate where in the disk the majority of the thermal dust emission originates for different wavelengths and how this relates to the local dust abundance and temperature. The red dashed line is for the flux at ∼4 µm, the red solid line for the flux at ∼1.3 mm (ALMA Band 6).

4.1.2 Spectral indices for detected sources

The spectral index of each detected source is listed in Table 6. We discuss the targets and serendipitous detections separately.

For two of the original targets, B243 and B275, we measure a spectral index consistent with a flux dominated by thermal dust emission from a circumstellar disk. The slightly lower index for B275 (α = 1.7 ± 0.48), could signify that there is some contribution from free-free emission, though the upper limit on the VLA 10 GHz measurement suggests that this is minimal (see also Fig. 3). For B243 we measure the highest index (α = 2.7 ± 0.61), making it the only target source for which we likely detect (partially) optically thin dust emission, possibly with a particle size distribution depleted in larger grains. The tentative detection of B268 in Band 6 is likely not related to a disk, since the source is undetected in Band 7, giving an upper limit to the spectral index (α < 1.2). B331 is the only object also detected in the VLA studies by Rodríguez et al. (2012) and Yanza et al. (2022). The flux values in those studies show an offset relative to each other (see Fig. 3), but the derived spectral index α = 0.7 ± 0.5 by Rodríguez et al. (2012) is consistent with our derived index between the two ALMA flux points (0.52 ± 0.48) and also with the index between ALMA Band 6 and VLA 10 GHz (0.53 ± 0.06). These indices clearly indicate that in B331 free-free emission dominates the measured flux, though a thermal dust contribution from a disk cannot be excluded.

Based on these observations, we assume for our further analysis that the fluxes of B243 and B275 are dominated by thermal dust emission from a circumstellar disk, and those of B268 and B331 by free-free emission. For the latter two objects (the upper limit of) the Band 7 flux was used for fitting. Their fit results and mass estimates based on ALMA data are placed in parentheses in Tables 7 and 8.

The serendipitous detections B331 NW and B331 SW have spectral indices well above 2, consistent with (partially) optically thin dust emission from a circumstellar disk. B275 SE (α = 1.4 ± 0.55) and B243 SW(α = 0.6 ± 0.8) have indices consistent with a large contribution from free-free emission, though, again, a contribution from dust in a disk cannot be excluded. We discuss the serendipitous detections in Sect. 5.3.

4.2 Fit results and disk mass estimates

The best fit values of the thin disk modeling (Sect. 3.5) are listed in Table 7 and the resulting models are shown in Fig. 3. The results for B331 and B268 are placed in parentheses because the fitted ALMA Band 7 flux point is either likely dominated by free-free emission (B331) or an upper limit (B268). Fitting simultaneously the short- and long-wavelength data points required a relatively shallow column density decrease (q > −1.5, except for B243), and high inner-disk surface densities (all best fit Σ_i are on or close to their maximum bound). This reflects that models with a nominal density structure (q = −1.5) and a small disk size (R_out ⪅ 50–70 au) that fit the short-wavelength flux underpredict the ALMA flux. We found the reverse of this effect in the PRODIMO models, which overpredict the short-wavelength fluxes while fitting the ALMA points.

The course of the dust-to-gas mass ratio in the disk as parameterized by R_turn and β (Eq. (5) and Fig. 2) differs per object. This is visualized in Fig. 4, which shows the dust temperature and best fit dust-to-gas mass ratio in the disk as a function of disk radius. Only for B331 is β high, and the dust-to-gas ratio approaches a step function with R_turn indicating where in the disk the dust abundance rises significantly. In all other cases β is small, so that the dust abundance rises approximately linearly, and R_turn indicates where it bends toward its limiting value at 10⁻². Figure 4 also shows the normalized cumulative flux in the NIR (4 µm) and in ALMA Band 6 (1.3 mm), illustrating where in the disk the emission originates at these wavelengths. From this we see that for B275, B243, and B268 the NIR flux is dominated by hot dust (∼1000 K) in the inner parts of the disk (≲5–10 au). For B331 there are two contributions: one from very tenuous hot dust in the inner disk, and a second from a small dense ring of cooler dust (∼300 K). The ALMA flux mostly originates from dust beyond R_turn, where the dust abundance is highest and the disk has more surface area. The emitting dust temperatures vary from object to object. For B275 and B331 the best fit models indicate a lack of extended cool dust (≲50–100 K), that is, the ALMA flux is dominated by higher temperatures (∼100–200 K). Since the flux of B268 is an upper limit, these results suggest that only B243 can be said to have an extended disk containing cool dust.

Table 8 summarizes all disk mass estimates for the four intermediate- to high-mass PMS stars. Again, the estimates for B331 and B268 that are based on ALMA fluxes are placed in parentheses, as the ALMA fluxes for these objects are likely dominated by free-free emission. The estimate based solely on the extrapolation of the inner disk (≲1 au ) CO bandhead fitting result leads to the lowest gas and dust masses (M_gas ∼ 10⁻⁵– 10⁻⁴ M_⊙). The best fit PRODIMO model leads to a mass estimate close to the values from the extrapolation (M_tot = 5 × 10⁻⁴ M_⊙). As this is the lowest-mass model calculated in the grid and it overestimates the NIR emission, this can be considered an upper limit in the context of estimates based on PRODIMO modeling.

The dust masses determined directly from the ALMA fluxes through Eq. (2) and from modeling the SEDs with the thin disk model lead to similar dust masses, that are at most an order of magnitude higher than the previous two estimates. The slightly lower dust mass estimate in Band 7 simply reflects how the measured spectral index differs from that assumed for the dust opacity (Sect. 3.3). The fit results illustrated in Fig. 4 suggest that the inner disks are dust depleted or, in the case of B331, that there is a dust gap⁴. In general, gas and dust gaps do not have to coincide (e.g., van der Marel 2023). By varying the dust-to-gas ratio separately from the gas column density, the models allow for the gas to decouple from the dust. If, instead, a gas-dust coupling was enforced or if the gas were even more depleted than the dust, this would result in lower gas masses than those presented in column 5 of Table 8. Thus, the gas masses derived from using the adapted P23 model to fit the SEDs can be seen as upper limits. This does not change the fact that the estimation heavily relies on the adopted results from the CO bandhead fitting. Finally, we recall that the nondetection of CO lines with ALMA suggests an upper limit on the gas mass of ∼2 × 10⁻⁵ M_⊙, in the optically thin case, or ∼6 × 10⁻³ M_⊙, when assuming an optically thick compact disk (Sect. 3.2).

5 Discussion

We present the first ALMA detections of disks around well- studied intermediate- to high-mass PMS stars at distances well beyond 1 kpc. An analysis of their properties reveals that they are mainly low-mass compact disks, possibly dust-depleted in the inner regions, that likely do not contribute significantly any more to the final stellar mass nor to the formation of any significant planetary systems (Sect. 5.1.1). Two of our target sources are dominated by free-free emission, indicating the presence of ionized material close to the star (Sect. 5.1.2). We place these findings in the context of trends observed in low-mass SFRs and among Herbig stars, and find tentative evidence for an influence of the massive star formation environment on disk formation and evolution (Sect. 5.2). We also report four serendipitous detections, which demonstrates the capabilities of ALMA to probe disks in high-mass SFRs (Sect. 5.3).

5.1 Detection of compact millimeter emission toward intermediate- to high-mass PMS stars

5.1.1 Detection of disks

Of the four intermediate- to high-mass PMS stars in M17 presented in Sect. 2.1, three were unambiguously detected with ALMA in Bands 6 and 7. B268 was marginally detected in Band 6, providing an upper limit on the Band 7 flux and spectral index. The spectral indices presented in Sect. 4.1.2 suggest that of the four detections, two (B275 and B243) are likely related to cool dust emission from a disk, while the other two (B331 and B268) have a significant contribution from free-free emission, limiting the constraints on possible outer disk masses to upper limits. Different approaches were used to estimate the upper limits to the dust and gas masses in the disks. For the gas mass, these estimates were partially based on earlier results from fitting CO overtone bandhead emission from the inner disks of the targets.

The results detailed in Sect. 4.2 and Fig. 4, obtained by fitting the SEDs with the adapted P23 model, suggest that the inner disks are dust depleted. In Sect. 3.5.1 we mention that the variation of the dust-to-gas ratio could also reflect a scenario in which the hot surface layers of the inner disk and the geometrical orientation determine the level of NIR flux rather than the amount of dust itself. The P23 1D disk model cannot distinguish between these scenarios (or a combination of them) and if the inner disks are not as dust depleted as the best fit models suggest, this could increase the dust mass. However, this increase will be minor since the issues only affect the inner part of the disk, while the bulk of the dust still resides in the outer parts. Moreover, the gas masses reported in Table 8 (and the remarks made on the results in Sect. 4.2) would remain unaffected.

The overall picture emerging from the reported results is that the studied objects are surrounded, if at all, by small circumstellar disks (R_out ≲ 50–70 au) with low dust masses (at most a few M_⊕). Estimations of the gas mass in the disks range from ∼10⁻⁵ to about ∼10⁻³ M_⊙. It is possible, even likely, that the disks are more compact than assumed in the modeling and optically thick in the CO lines as discussed in Sect. 3.2. If this is the case, the gas mass might be as high as ∼6 ×10⁻³ M_⊙. In any case, the low disk masses imply that all stars in the sample have essentially finished their main accretion phase; that is, they have stopped accreting amounts of mass that still could contribute significantly to the final mass of the central star. This is in line with the well-detected photospheres of these objects and the previously reported lack of jets, outflows, or other signs of strong accretion bursts (Derkink et al. 2024).

In this context the presence of CO overtone bandhead emission in the spectra, which is often associated with accretion and earlier stages of formation (see Poorta et al. 2023, and references therein), may be surprising. In addition, the hydrogen emission lines in the spectra of these targets are often linked to accretion and are used to estimate accretion rates in Herbig stars (e.g., Wichittanakom et al. 2020). The accretion rates reported for objects with CO bandhead emission are very diverse, ranging from 10⁻³ to 10⁻⁸ M_⊙ yr⁻¹ (Poorta et al. 2023). Accretion rates on the order ∼10⁻⁷–10⁻⁸ M_⊙ yr⁻¹ are consistent with the mass upper limit estimates.

Our modeling suggests that the inner disks are dust depleted, reminiscent of gap-like or ring-like structures in transition disks (TDs) among lower-mass PMS stars (van der Marel 2023). Herbig stars are traditionally classified into Group I and Group II disks, originally based on their SED appearance (Meeus et al. 2001). The basic characteristics of these groups are summarized in Table 9⁵. Group I disks can be understood to be the high- mass equivalent of T Tauri TDs, with a large inner cavity and a millimeter-bright, extended, massive outer disk component. Classification in terms of Group I or II for our objects is ambiguous. Based on their near- to mid-infrared colors and their inner dust gap, B275 and B331 qualify as Group I. However, in terms of their low millimeter-brightness and their limited size, all the objects in our sample qualify as Group II.

We discuss the lack of extended massive disks in millimeter emission further when comparing our results to previously studied disk populations in Sect. 5.2.

5.1.2 Free-free emission

While some contribution of cold dust to the millimeter emission of B331 (and possibly B268) cannot be excluded, the spectral indices indicate a dominant role for free-free emission. The spectral indices found for B331 (0.5–0.7) are consistent with a stellar wind or disk wind. However, (a contribution of) other mechanisms such as optically thin free-free emission from ionized gas close to the star or the base of a radio jet cannot be excluded (see Sect. 4.1). With a mass of ∼10 M_⊙, B331 is the most massive source in the sample. There are no indications of a strong stellar wind in the X-shooter spectra (Ramírez-Tannus et al. 2017). As can be seen in the SED of B331 (top right panel of Fig. 3) it lacks the NIR excess seen in the other sources, the infrared excess only starts at ∼4–5 µm. This translates into the gap-like feature modeled in the dust abundance (Fig. 4); that is, the dust is absent up to ∼45 au, after which the dust-to-gas ratio steeply rises to 10⁻². Backs et al. (2023) also model this gap and, in the absence of long wavelength data, classify the object as a Group I source. All in all, the spectral index and the derived disk morphology combined could suggest a photo-evaporative disk wind. This is consistent with photo-evaporation mechanisms where the disk is being sublimated farther away from the star and a mid-disk gap works its way both outward and inward (Alexander et al. 2006; Hollenbach et al. 1994). Similar inner gaps have been observed for massive YSOs by Frost et al. (2021), who also attribute this, at least in part, to photo-evaporation.

The fact that B268 is only marginally detected indicates that the dust disk of this object is more dissipated than that of the other objects. It is noteworthy that both its stellar characteristics and its infrared SED are very similar to those of B243, for which a millimeter disk is most clearly detected (Table 1 and Fig. 3). It is also interesting that the CO bandhead emission for this source is very pronounced (Poorta et al. 2023), but that the estimated mass of its disk based on this emission is the lowest of all the objects in the sample (Cols. 5 and 6 in Table 8). This emphasizes that the NIR view of a disk does not necessarily predict its properties probed at longer wavelengths. The most likely explanation of the observations of this source is still that its disk is more compact than that of the other detected sources.

5.2 Disk masses and sizes compared to other disk populations

5.2.1 Low-mass star-forming regions and Herbig stars

In Figs. 5 and 6 we compare our objects with samples from the low-mass SFRs Lupus (Ansdell et al. 2016) and Upper Sco (Barenfeld et al. 2016), and with a sample of Herbig AeBe stars within a distance of ∼450 pc studied by Stapper et al. (2022). The dust masses for all the objects in the quoted studies were determined from ALMA Band 6 and/or 7 continuum flux using Eq. (2), although sometimes with different assumptions for the dust temperature than in this study. For the purpose of comparison, we therefore use the dust mass derived from Eq. (2), and because of the possible free-free contamination discussed earlier we use (the upper limits on) the dust masses based on the Band 7 flux.

In the four panels of Fig. 5 we show the dust mass in the disk versus stellar mass (upper left), age (upper right), stellar luminosity (lower left), and outer dust radius (lower right). The parameters used for this comparison are summarized in Table 10. Compared to the other samples, our sample is small. The purpose of this comparison is therefore to see to what extent the higher- mass objects in M17 follow previously identified trends.

In low-mass SFRs two main trends have been observed with respect to dust masses among Class II disk populations (Lada 1987). The first is an increase with stellar mass, and the second is a decrease with stellar or region age (e.g., Pascucci et al. 2016; Ansdell et al. 2017; Barenfeld et al. 2016). The stellar mass trend can be seen in the Lupus and Upper Sco data points in the upper left panel of Fig. 5; the age trend is seen when comparing between those regions globally (upper right panel). The Herbig sample of Stapper et al. (2022) follows the positive correlation with stellar mass, but not the negative correlation with older age. The authors explain this by suggesting that the effects of mass and age can be degenerate. As also found in previous studies (e.g., van der Marel & Mulders 2021), not only do stars with higher mass start out with a more massive disk, but they also appear to retain this mass for a longer time. van der Marel & Mulders (2021) find that the existence of dust traps, caused by structures such as an inner gap or ring-shaped clearing, is necessary to retain an observable millimeter dust reservoir for a longer time. Both Stapper et al. (2022) and van der Marel (2023) propose that such structures may form more readily in disks with higher mass.

There is a large spread in dust mass in the shown samples, and as such the dust masses of the M17 objects are not outliers. However, the M17 sample studied here is both the most massive and the youngest sample in question, and yet of all the samples it has the lowest mean disk mass, contrary to the trends described above. The mean values are listed in Table 11 and are shown in Fig. 5. In addition, the models suggest a millimeter dust gap, but for the M17 objects this does not come with an extended high-mass millimeter dust disk. On the other hand, in the intermediate- to high-mass range the evolutionary timescale of the star may become small with respect to the disk evolution timescales of low-mass stars. Higher-mass stars may sooner deplete their disks due to short formation timescales, high accretion rates, strong ultraviolet (UV) illumination, and outflows from the star and/or disk. To allow for these effects in the age comparison across the stellar mass ranges, in Fig. 6 we plot the age of the objects “normalized” by the PMS lifetime (i.e., zeroage main-sequence (ZAMS) age) corresponding to their stellar mass. The dust mass in this plot is normalized by the stellar mass to better isolate the effects of age. Almost all the plotted objects appear on their PMS (age/ZAMS-age <1). As may be expected, the M17 objects are further in their PMS evolution than the low-mass star populations, where the PMS lifetime is (very) long (10⁷ to almost 10⁹ Myr) compared to disk evolution timescales. Compared to the Herbig sample from Stapper et al. (2022) the M17 objects are similar or slightly earlier in their evolution. Though there are some Herbig stars that fall in the same range, with respect to all the samples the stellar-mass-normalized dust masses of the M17 sample are very low.

The age plot in Fig. 5 includes the serendipitous discoveries which, as discussed in Sect. 5.3, are likely lower-mass YSOs with disks comparable in mass to the original targets. As argued in that section, considering the density of YSOs in the region, the lack of higher-mass disk discoveries may hint at a disk population with relatively low masses, especially for a SFR as young as M17.

One possible explanation for these findings is that the formation environment plays a crucial role in the initial disk mass distribution and in the disk evolution and lifetime. This explanation is supported by studies in other massive SFRs. Eisner et al. (2018), for example, studied the disk population (92 disks) in the Orion Nebula Cluster and found compact disks (<60 au), with low masses and with only a weak correlation between stellar mass and disk mass compared to low-density SFRs (see Table 11 for the mean disk and stellar masses of their sample). They suggest that the disks are affected by photo-ionization and, possibly, stellar encounters. Similar population studies in Orion suggest that even intermediate far-UV radiation fields from A0 and B stars have a significant impact on the evolution of protoplanetary disks (van Terwisga et al. 2020; van Terwisga & Hacar 2023). Investigating mass and accretion in the disk population of the Lagoon Nebula (M8), Venuti et al. (2024) also find evidence for faster disk evolution in the central regions with respect to the outskirts of this massive SFR.

Finally, the upper limits on the disk radii are consistent with the observed trend between disk mass and radius in Stapper et al. (2022). Given the discussion above it would appear likely that the detections comprise optically thick compact disks (⪅30 au), which are able to shield themselves to survive the high UV-radiation environment. The mass estimate even in this case remains rather low. Next to external evaporation and encounters in the dense environment after disk formation, this might also be due to processes such as competitive accretion limiting the amount of mass the disks could assemble during disk formation (e.g., Bate 2018).

Fig. 5 Disk dust mass as a function of stellar mass (upper left), age (upper right), stellar luminosity (lower left), and outer dust radius (lower right), based on ALMA continuum data for different samples. The low-mass PMS populations of the Lupus (Ansdell et al. 2016) and Upper Sco (Barenfeld et al. 2016) SFRs are shown in orange and green, respectively. For these populations the age range of the SFR is indicated by plotting each object as a horizontal line in the age plot. The Herbig Ae/Be sample from Stapper et al. (2022) is shown in blue and gray. The gray points are special cases in their sample that are both higher mass and younger; the labeled objects are discussed in the text. In black we show the PMS stars in M17 studied in this work. In the upper panels the mean dust mass, and the mean stellar mass (left) and age (right) of each sample are indicated with horizontal and vertical dashed lines. Upper limits on the dust mass are indicated by down-pointing triangles. Upper limits on the outer dust radius (lower right) are indicated by left-pointing triangles.

Fig. 6 Normalized disk dust mass vs. normalized age for the same samples as in Fig. 5. The dust mass is normalized to the stellar mass of each object. The age is divided by the PMS lifetime of a star of the respective mass, set by the ZAMS age from the MIST PMS tracks (Dotter 2016). For the objects in Lupus (orange) and Upper Sco (green), for which no individual age is known, the age range is indicated by plotting each object as a horizontal line. The dashed lines indicate the mean values of both quantities for the different populations.

5.2.2 Herbig stars of higher mass

In addition to the influence of the SFR environment on the disks, another (not mutually exclusive) explanation for the low disk masses of the M17 targets could be a faster disk evolution and dissipation in the higher stellar mass ranges. So far, in our comparison we have taken into account the bulk of the Herbig stars in Stapper et al. (2022), which have a mean mass of ∼2 M_⊙ (blue marks in Fig. 5; see also Table 11), significantly lower than the masses of the M17 targets. Stapper et al. (2022) marked a few objects in their sample as exceptions (gray marks in the plots) because their Herbig nature was ambiguous, three of these due to their high mass or very young nature. These three exceptions (MWC297, ZCMa, and HD58647) are labeled in the plots because they are similar to the M17 objects in age, mass, and luminosity (lower left panel in Fig. 5 and Table 12).

The small number of objects of this nature emphasizes the dearth of YSOs in the higher mass ranges (≳4 M_⊙) that have been subject to multiwavelength studies. Because of this, it is not possible at the moment to derive statistical properties, and the comparison remains on an object-to-object basis. The difference between the labeled sources and the M17 targets is that the former are relatively isolated objects, while the latter are located in a dense cluster environment. We now discuss each labeled source separately.

The two objects, ZCMa and HD58647, are similar in mass (∼3.8 M_⊙); they both have outflows associated with accretion. ZCMa has a FU Orionis variable companion and an asymmetric outflow (Baines et al. 2006). HD58647 is observed to have a disk wind from an accretion disk (Kurosawa et al. 2016) and its dust mass is similar to the M17 objects. This serves as an example that low observable dust mass and ongoing accretion can both be present.

MWC297 is a high-mass star (∼14.5 M_⊙) in the Herbig catalog of Vioque et al. (2018), but was earlier classified as a ZAMS star of ∼10 M_⊙ by Drew et al. (1997). Manoj et al. (2007) detected a remnant formation disk (<80 au) with dust mass similar to that reported by Stapper et al. (2022) (∼100 M_⊕). They note the possibility of a disk wind or stellar wind based on a VLA spectral index of ∼0.6. Interestingly, while a hot gas component (∼1500 K) from the inner disk (∼12 au) was detected in CO ro-vibrational fundamental and overtone emission, no cold gas was detected in CO lines at millimeter wavelengths (Sandell & Vacca 2023; Manoj et al. 2007).

MWC297 is very similar to B331 in stellar mass, age, and luminosity. Moreover, both have CO overtone bandhead emission in the NIR, both lack cold CO gas, and they have similar VLA spectral indices. The only significant difference is that the dust disk of MWC297 is almost two orders of magnitude higher in mass. Perhaps this points, again, to a role for the environment in the disk evolution of high-mass YSOs.

5.3 Serendipitous discovery of sources

In the four fields surrounding the original targets, four serendipitous detections were made in both ALMA Bands. All these radio sources have NIR counterparts. We collected NIR and MIR photometry from the VizieR Catalogues⁶. The flux points are listed in Table C.1 and plotted in Fig. 3. The source positions and projected distances from the closest original target are listed in Table 5.

We estimated the expected number of extra-galactic background sources as a function of field-size by adopting Eq. (A10) from Anglada et al. (1998) with the ALMA 12m primary beam θ_A ≈ 20.6″ × (300/ν) GHz. This yields (similar to Eq. (A11) from Anglada et al. 1998) $$\begin{array}{l}<N>=0.00196\left\{1-\exp \left[-0.00123{\left(\frac{{\theta }_{\text{F}}}{\mathrm{arcsec}}\right)}^2{\left(\frac{v}{300\text{\hspace{0.17em}GHz}}\right)}^2\right]\right\}\\ {\left(\frac{v}{300\text{\hspace{0.17em}GHz}}\right)}^{-2}{\left(\frac{ {\text{S}}_0}{\text{mJy}}\right)}^{-0.75}{\left(\frac{v}{5\text{\hspace{0.17em}GHz}}\right)}^{-0.525},\end{array}$$

with S₀ = 0.165 mJy the minimum detectable flux density at the field center, taken to be five times the sensitivity in Band 6; ν = 225 GHz the Band 6 frequency; and θ_F = 14″ the field diameter. We find the expected number of extra-galactic sources in each field to be negligible, <N > ≈ 2 × 10⁻⁴. We now discuss the known characteristics for each source.

The two sources B275 SE and B243 SW are listed in the Massive Young star-forming Complex Study in Infrared and X- rays (MYStIX; Broos et al. 2013) catalog as an X-ray (Chandra) source with NIR counterpart, and classified as a “young star in a massive SFR.” Both objects also feature in a study by Yang et al. (2022), who use machine learning techniques to assign a probability (P) to the classification of X-ray sources based on multiwavelength features (X-ray, NIR-MIR photometry and colors). They classify B275 SE as a low-mass X-ray binary (LMBX) with P = 0.4 or as a YSO with P = 0.3. Considering the evolved nature of LMBXs and the young age of M17, the first classification seems rather unlikely. B243 SW is classified as a YSO with P = 0.75. This object is also labeled as a variable source in Gaia DR3 Gaia Collaboration (2023).

Based on these findings, both B275 SE and B243 SW are likely low-mass YSOs. The spectral indices (α = 1.4 ± 0.55 and α = 0.6 ± 0.8, respectively) indicate a contribution from free-free emission, for B275 SE probably next to dust emission from a circumstellar disk. In the case of B243 SW, the presence of both variability and X-rays may indicate that the free-free emission is related to an outflow driven by ongoing accretion (see, e.g., Rota et al. 2024). However, even without correcting for the contribution of free-free flux, the derived dust mass for this detection is rather low (see Table 6).

B331 NW and B331 SW have spectral indices well above 2, consistent with (partially) optically thin dust emission from a circumstellar disk. B331 NW (α = 2.3 ± 0.6) has a relatively low projected separation from B331 (1.4″, ∼2300 au) and is therefore not resolved in most photometric surveys. In the reported data from Pan-STARRS1 (Chambers et al. 2016) and MYStIX crowded fields (King et al. 2013) it is barely resolved, and the measured fluxes are likely contaminated by B331. It is clearly observed as a distinct source in the previously mentioned acquisition images from LBT. B331 SW (α = 2.4 ± 0.5) is the faintest detection in the NIR with only one survey (King et al. 2013) reporting fluxes, yet it is the second brightest detection (after B331) with ALMA. Perhaps this object is of lower mass, partly obscured by its disk, or simply more embedded.

All in all, it is likely that all four serendipitous detections relate to (lower-mass?) YSOs and that at least two have a dust disk. None of these objects is resolved, and the measured fluxes are similar to those of the original targets, leading to similar values for (the upper limits on) the dust mass and radius. For the reported masses in Table 6 a temperature of T = 150 K was used. Assuming lower dust temperatures for these sources (e.g., T_dust = 20 K) would result in dust mass estimates about one order of magnitude higher (i.e., a few tens of M_⊕), closer to the averages found in low-mass SFRs (Fig. 5). However, for the two objects with free-free emission this would still represent an upper limit.

The overall picture from these detections, under the assumption that they are YSOs, is that their disks are small and on the lower-mass end for an age of ∼0.6 Myr (see upper right panel Fig. 5). While the discoveries testify to the source density in the region, it is remarkable that brighter disks that would have been detected if present in these fields, are lacking. The number of sources is too small to be conclusive, but it is not unlikely that we are catching the brightest (i.e., most massive) disks and that deeper observations would yield an even fainter population. As suggested before for the original targets it is possible that the disks are compact, optically thick, and more massive than derived under an optically thin assumption. That would support the notion that in a dense high-UV environment such as M17 perhaps this is the kind of disk that survives, when disks survive at all.

6 Summary and conclusions

In this paper we presented ALMA millimeter continuum detections of four intermediate- to high-mass PMS stars in the massive star-forming H II region M17 (distance 1.7 kpc), and report on the serendipitous discovery of four additional sources, likely low-mass YSOs. All the detected sources are unresolved at a resolution of ∼0.07″, constraining the outer radius of the emission to a maximum of ∼60 au. We used the spectral index between ALMA Bands 6 and 7, and (the upper limits on) the VLA centimeter flux to determine the origin of the millimeter emission, that is, free-free radiation or thermal dust emission. We derived (upper limits on) the dust masses from the ALMA continuum flux and an upper limit on the gas mass from the nondetection of the rotational ¹²CO (J = 2–1) line. We combined the ALMA data with near- to mid-infrared photometry and spectroscopy and applied different models to estimate the total disk masses of the target sources. In the following we summarize our results:

• The disks around the four target sources as well as the four serendipitous discoveries contain a dust mass of a few M_⊕ at most;

• Estimates of the upper limit on the total gas mass in the disks around the intermediate- to high-mass PMS stars vary between ∼2 × 10⁻⁵ M_⊙ and, when assuming an optically thick compact disk, ∼6 × 10⁻³ M_⊙ ;

• Spectral indices suggest that the fluxes from two target sources and two serendipitous sources have a significant contribution from free-free emission from ionized material. In the higher-mass sources (B331 and B268) this may indicate a small-scale H II-region, the base of a radio jet, or a photo-evaporative disk wind. In the likely low-mass sources (B275 SE and B243 SW) the combination with X-ray emission may indicate an accretion driven outflow;

• Our modeling suggests that the inner disks of some of the targets are dust depleted, reminiscent of gap- or ring-like structures, despite the presence of hot and dense gaseous disks, as evidenced by CO bandhead emission from these objects;

• Comparison to disk population studies in low-mass SFRs and among Herbig stars indicates that our sample is both the most massive and the youngest. However, contrary to commonly observed trends, it has the lowest mean disk dust mass.

All in all, we conclude that the studied targets are surrounded by low-mass, compact disks that likely no longer offer a significant contribution to either the final stellar mass or the formation of a planetary system. Along with the four serendipitous discoveries, our findings offer tentative evidence of the influence of the massive star formation environment on disk formation, lifetime, and evolution. The detections presented hint at a rich YSO population in M17, with disks detectable by ALMA. This, along with the thousands of known infrared excess sources in M17 (Lada et al. 1991; Jiang et al. 2002), offers the perspective of disk population studies that not only include a significant number of intermediate- to high-mass targets, but also allows their study alongside low-mass YSOs in a massive SFR.

Data availability

Supplementary material is available on https://doi.org/10.5281/zenodo.14501009.

Acknowledgements

The authors express their gratitude to the anonymous referee for their helpful comments and insights toward improving this paper. JP acknowledges support from NWO-FAPESP grant 629.004.001 (PI L. Kaper). This paper makes use of the following ALMA data: 2019.1.00910.S, 2018.1.01091.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. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the VizieR catalog access tool, CDS, Strasbourg, France (DOI : 10.26093/cds/vizier). CHR acknowledges the support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Research Unit “Transition discs” – 325594231 and the support by the Excellence Cluster ORIGINS which is funded by the DFG under Germany’s Excellence Strategy – EXC-2094 – 390783311. CHR is grateful for support from the Max Planck Society.

Appendix A Observational details

Appendix B Full fitting results

In Tables B.1 and B.2 we provide the detailed results of the 2D Gaussian fitting procedure described in Sect. 2.3. Where the values for the deconvolved major/minor axes are missing, the CASA fitting task could not deconvolve the source from the beam, due to insufficient resolution. The cases where the deconvolved axes are provided, the resulting values are significantly smaller than the corresponding beam axis and likely not representative of true source sizes.

Appendix C Photometry

Appendix D Detected source images

Figure D.1 provides zoomed-in images around all the detected sources.

Fig. D.1 Zoomed-in images of all detected sources. The target sources are labeled in red; the serendipitous detections are labeled in white, and are named after the target in whose image they are detected (see Sect. 2.3). The flux scales of each image are slightly different. Contour levels: 3, 5, 8, 10 × noise rms (between 0.022 – 0.040 mJy in Band 6 and 0.056 – 0.083 mJy in Band 7). No radio continuum flux is detected for B268 in Band 7 and B331 NE in Band 6; they are included as illustrations.

References

All Tables

All Figures

Fig. 1 Overview of the M17 region in different wavelengths. The left panel shows a NIR color composite image of M17 based on 2MASS data: J (blue), H (green), and K (red), with large-scale 1.3 mm continuum emission contours (archival ALMA Band 6 data, project code 2018.1.01091.S). Contour levels: 5, 15, 30, and 70 × noise rms (respectively 0.018, 0.053, 0.11, and 0.25 Jy). The middle panel shows only the ALMA mosaic, with the beam size in the bottom left corner. In both panels the rectangular region indicates the field in which all four targets in this study (indicated in red) are found. The rightmost panel is a zoomed-in image of this region with all targets and serendipitous discoveries marked in black. To give an impression of the scales, the radius of each object dot is 1″, which corresponds to five times the image size of each panel in Fig. D.1.

In the text

	Fig. 2 Illustration of the logistic function (Eq. (5)) used for varying the dust-to-gas mass ratio Md/Mg in the disk. The limiting values are fixed to the same values as in the models used for fitting: (Md/Mg)base = 10−8 and (Md/Mg)lim = 10−2. The value of the free parameter β is varied here to illustrate its effect.
In the text

Fig. 3 Best fit SEDs, NIR and MIR photometry points, and ALMA and VLA flux points for our four original targets and four sources serendipitously detected in the covered fields. Each panel shows the flux data for the original target and the Kurucz model for its stellar continuum in blue. The best fit model SED from the thin disk modeling (described in Sect. 3.5) is plotted in orange. Each panel also shows the flux data for the serendipitous detection(s) made in the field of the respective original target (in green and red); these data are not fitted. The downward arrows indicate the upper limit derived for nondetections (mostly VLA 10 GHz). Only B331 (upper right panel) is detected at VLA wavelengths. The dotted lines indicate the slopes between the radio flux points for this object, with a spectral index ∼0.67 for the VLA points at 4.96 GHz and 8.46 GHz (Rodríguez et al. 2012), and ∼0.52 between the ALMA and VLA 10 GHz flux points (Yanza et al. 2022). There is an offset between the data from Rodríguez et al. (2012) and (Yanza et al. 2022).

In the text

Fig. 4 Dust temperature (blue) and dust-to-gas mass ratio (black) as a function of disk radius resulting from the thin disk model fit (described in Sect. 3.5) for each object. The dotted black vertical line indicates the fitted inflection point Rturn of the logistic function of the dust-to-gas mass ratio. In red is the normalized cumulative flux distribution to illustrate where in the disk the majority of the thermal dust emission originates for different wavelengths and how this relates to the local dust abundance and temperature. The red dashed line is for the flux at ∼4 µm, the red solid line for the flux at ∼1.3 mm (ALMA Band 6).

In the text

Fig. 5 Disk dust mass as a function of stellar mass (upper left), age (upper right), stellar luminosity (lower left), and outer dust radius (lower right), based on ALMA continuum data for different samples. The low-mass PMS populations of the Lupus (Ansdell et al. 2016) and Upper Sco (Barenfeld et al. 2016) SFRs are shown in orange and green, respectively. For these populations the age range of the SFR is indicated by plotting each object as a horizontal line in the age plot. The Herbig Ae/Be sample from Stapper et al. (2022) is shown in blue and gray. The gray points are special cases in their sample that are both higher mass and younger; the labeled objects are discussed in the text. In black we show the PMS stars in M17 studied in this work. In the upper panels the mean dust mass, and the mean stellar mass (left) and age (right) of each sample are indicated with horizontal and vertical dashed lines. Upper limits on the dust mass are indicated by down-pointing triangles. Upper limits on the outer dust radius (lower right) are indicated by left-pointing triangles.

In the text

Fig. 6 Normalized disk dust mass vs. normalized age for the same samples as in Fig. 5. The dust mass is normalized to the stellar mass of each object. The age is divided by the PMS lifetime of a star of the respective mass, set by the ZAMS age from the MIST PMS tracks (Dotter 2016). For the objects in Lupus (orange) and Upper Sco (green), for which no individual age is known, the age range is indicated by plotting each object as a horizontal line. The dashed lines indicate the mean values of both quantities for the different populations.

In the text

Fig. D.1 Zoomed-in images of all detected sources. The target sources are labeled in red; the serendipitous detections are labeled in white, and are named after the target in whose image they are detected (see Sect. 2.3). The flux scales of each image are slightly different. Contour levels: 3, 5, 8, 10 × noise rms (between 0.022 – 0.040 mJy in Band 6 and 0.056 – 0.083 mJy in Band 7). No radio continuum flux is detected for B268 in Band 7 and B331 NE in Band 6; they are included as illustrations.

In the text