Issue |
A&A
Volume 671, March 2023
|
|
---|---|---|
Article Number | A13 | |
Number of page(s) | 18 | |
Section | Interstellar and circumstellar matter | |
DOI | https://doi.org/10.1051/0004-6361/202244817 | |
Published online | 28 February 2023 |
Massive pre-main-sequence stars in M17
Modelling hydrogen and dust in MYSO disks
1
Anton Pannekoek Institute for Astronomy, University of Amsterdam,
Science Park 904,
1098
XH Amsterdam,
The Netherlands
e-mail: f.p.a.backs@uva.nl
2
Institute of Astrophysics, Universiteit Leuven,
Celestijnenlaan 200 D,
3001
Leuven,
Belgium
3
Max Planck Institute for Astronomy,
Königstuhl 17,
69117
Heidelberg,
Germany
4
Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München,
Scheinerstr. 1,
81679
München,
Germany
5
Max-Planck-Institut für extraterrestrische Physik,
Giessenbachstrasse 1,
85748
Garching,
Germany
6
Kapteyn Astronomical Institute, University of Groningen,
Postbus 800,
9700
AV Groningen,
The Netherlands
Received:
26
August
2022
Accepted:
20
December
2022
Context. The young massive-star-forming region M17 contains optically visible massive pre-main-sequence stars that are surrounded by circumstellar disks. Such disks are expected to disappear when these stars enter the main sequence. The physical and dynamical structure of these remnant disks are poorly constrained, especially the inner regions where accretion, photo-evaporation, and companion formation and migration may be ongoing.
Aims. We aim to constrain the physical properties of the inner parts of the circumstellar disks of massive young stellar objects B243 (6 M⊙) and B331 (12 M⊙), two systems for which the central star has been detected and characterized previously despite strong dust extinction.
Methods. Two-dimensional radiation thermo-chemical modelling with PRODIMO of double-peaked hydrogen lines of the Paschen and Brackett series observed with X-shooter was used to probe the properties of the inner disk of the target sources. The model was modified to treat these lines. Additionally, the dust structure was studied by fitting the optical and near-infrared spectral energy distribution.
Results. B243 features a hot gaseous inner disk with dust at the sublimation radius at ~3 AU. The disk appears truncated at roughly 6.5 AU; a cool outer disk of gas and dust may be present, but it cannot be detected with our data. B331 also has a hot gaseous inner disk. A gap separates the inner disk from a colder dusty outer disk starting at up to ~100 AU. In both sources the inner disk extends to almost the stellar surface. Chemistry is essential for the ionization of hydrogen in these disks.
Conclusions. The lack of a gap between the central objects and these disks suggests that they accrete through boundary-layer accretion. This would exclude the stars having a strong magnetic field. Their structures suggest that both disks are transitional in nature, that is to say they are in the process of being cleared, either through boundary-layer accretion, photo-evaporation, or through companion activity.
Key words: circumstellar matter / stars: massive / stars: pre-main sequence / stars: formation
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
The evolution of high-mass protostars and their descendants, the massive young stellar bbjects (MYSOs), is expected to be very fast, ~105 yr (e.g. Hosokawa & Omukai 2009), and, until very late in the build-up process, hidden from view as it unfolds deep within dusty natal clouds. Though much is still unclear about the mechanism that leads to the assembly of a massive star, most theories agree on the need for a dense and massive accretion disk (e.g. Tan et al. 2014; Beltrán & de Wit 2016), possibly in combination with bipolar polar outflows at relatively early phases of formation. Likely, the massive disks that develop are complex in structure with multiple processes acting to cause instabilities, for example magneto-rotational instabilities leading to turbulence and associated angular momentum redistribution (e.g. Balbus & Hawley 1998) and gravitational instabilities giving rise to the formation of companions (e.g. Kratter et al. 2010; Oliva & Kuiper 2020).
The inner parts of the disk are particularly interesting for at least two reasons. First, from these regions the accretion of mass and redistribution of angular momentum takes place. Important open problems are the gas content of the inner disk reservoir and the accretion efficiency of this material close to the arrival of the star on the main sequence. The assembly of angular momentum is particularly interesting in view of the buildup and outcome of the surface rotation velocity, the latter property being key to reaching a new level of understanding of main sequence single and binary star evolution (Wang et al. 2022). Second, the mechanism(s) disrupting the circumstellar disk likely act first fairly close to the central star. The main candidates for this process are disk winds, outflows resulting for example from photo-evaporation (e.g. Gorti & Hollenbach 2009; Owen et al. 2010; Gangi et al. 2020), stellar winds, and stellar or planetary companion formation (e.g. Müller & Kley 2013). Our understanding of disks that are in the process of clearing themselves out (so-called transitional disks) stems mostly from low- and intermediate-mass stars, where the process takes 10–20% of the total YSO lifetime (Furlan et al. 2011). Given the short duration of the MYSO phase (104 – 105 yr; Hosokawa & Omukai 2009), it remains to be seen whether such a modest fraction of the MYSO lifetime will suffice, or whether perhaps the majority or even all MYSO disks are transitional in nature.
Atomic recombination and CO bandhead emission lines are the main (gas-phase) diagnostics of the properties, structure, and kinematics of the inner dust-free regions of the disks of MYSOs. The hydrogen line-formation mechanism near the end of the MYSO phase is complex, as the stars are likely not accreting enough to cause significant hydrogen ionization as a result of viscous heating and shocks, and the central stars' photospheres are not hot enough to produce sufficient ionizing radiation and associated hydrogen photo-ionization.
Evidence is mounting that the hydrogen emission originates from very close to the stellar surface. Koumpia et al. (2021) present an interferometric K-band study of a sample of MYSOs spatially resolving the dust and Brγ emission in objects ranging from 9 to 16 M⊙. They find Brγ emission to originate from ~0.9 to 6 AU from the central star and the 2.2 μm emission to originate from ~1.4 to 6.8 AU from the central star depending on the object. They conclude that the hydrogen emission originates from a more compact region than the 2.2 μm dust emission, with the continuum dust emission being consistent with originating from the dust sublimation radius. This is in line with the findings of Caratti o Garatti et al. (2016) who find a Brγ emitting region of 6–13 AU and a further out continuum emitting region of 17 AU for a 20 M⊙ MYSO. Additionally, they detected an outflow in Brγ. Kraus et al. (2010) also find a dust free cavity analogous to those found in lower mass YSOs compatible with the dust sublimation radius. Gravity Collaboration (2020) detect CO emission in the inner gaseous disk of the massive (15 M⊙) YSO NGC 2021 IRS 2, which is well within the dust sublimation radius. The dynamical nature of this hot gas close to the central star is still a matter of debate. MYSOs could have weak or absent magnetic fields similar to their main-sequence counterparts in which few stars are observed to have strong magnetic fields (e.g. Fossati et al. 2015; Wade et al. 2016). Thus accretion might proceed through boundary layer (BL) accretion rather than magnetospheric accretion. Therefore, disks might potentially extend all the way to the stellar surface, becoming dust free inside of the dust sublimation radius. Alternatively, the nearby gas could result from bipolar outflows (Frost et al. 2021).
On the basis of infrared spectroscopy, photometry, and high-spatial resolution imaging of optical scattered light from small dust grains, disks around YSOs have been classified into two groups (e.g. Garufi et al. 2017), reflecting large-scale morphological properties. Group I sources show large extended disks with a pronounced gap in the centre. This results in low or lacking near-infrared excess due to the absence of hot dust. Such inner dust-free cavities extending beyond the dust sublimation radius are observed in MYSOs (Frost et al. 2021). Group II sources show strong NIR-excess, but limited excess at longer wavelengths. This suggests the presence of hot dust near the star, but the absence of cooler dust further out, which can be caused by a small truncated disk, or a self-shadowed disk. Though this classification does clearly imply a different structure of the inner disk (where the hot dust is located), it does not link to a potential gas reservoir in these inner regions. Scrutinizing this hot gas, if present, may help to understand the processes shaping the inner disk regions.
In this paper we present the modelling of hydrogen emission lines and dust continuum emission from two MYSOs using the radiation thermo-chemical code PRODIMO (Woitke et al. 2009, Woitke et al. 2016; Kamp et al. 2010, Kamp et al. 2017; Thi et al. 2011). We fitted X-shooter data of the hydrogen lines (Ramírez-Tannus et al. 2017) and archival photometry of the dust emission for two sources in M17. The H II region containing the stars and the data are described in Sect. 2. In Sect. 3 we introduce the model and the modifications applied to the model to enable our analysis. Section 4 shows the results; these are discussed in Sect. 5. Conclusions are outlined in Sect. 6.
2 Description of the data
M17 is one of the most luminous H II regions in the Galaxy (3.6 × 106 L⊙; Povich et al. 2009). It is located in the Carina-Sagittarius spiral arm at a distance of 1.98 ± 0.14 kpc (Xu et al. 2011). The young cluster NGC 6618 (~1 Myr; Hanson et al. 1997; Broos et al. 2007; Povich et al. 2009), containing 16 O-type stars and over 100 B-type stars, is located at the centre of the H II region (Chini et al. 1980; Hoffmeister et al. 2008), while the surrounding molecular cloud hosts pre-main-sequence (PMS) stars (Hanson et al. 1997; Ramírez-Tannus et al. 2017). On the basis of a near-infrared (NIR) excess in the spectral energy distribution (SED) and the presence of hydrogen and CO bandhead emission lines, Hanson et al. (1997) identified a sample of massive (6-20 M⊙) young stellar object (YSO) candidates embedded in the molecular cloud. Using X-shooter (Vernet et al. 2011) mounted on the ESO Very Large Telescope, Ramírez-Tannus et al. (2017, hereafter RT17) carried out a spectroscopic follow-up of 12 of these objects, providing the spectroscopic confirmation that these objects are massive pre-main-sequence stars. The stellar parameters were determined from photospheric line fitting or SED fitting. In the nine cases where the objects have visible photospheres, they modelled the stellar spectrum using the non-LTE radiative transfer code FASTWIND (Puls et al. 2005; Rivero González et al. 2012) together with the genetic algorithm PIKAIA (Mokiem et al. 2005; Charbonneau 1995). Otherwise SED fitting was used to determine stellar properties. The obtained temperatures, luminosities, and projected rotational velocities allowed them to place these MYSOs in the Hertzprung-Russell Diagram (HRD) and compare their position with pre-main-sequence (PMS) evolutionary tracks from Hosokawa & Omukai (2009). Two of the observed objects (B111 and B164) have spectra characteristic of O-type stars and are located at the zero-age main-sequence (ZAMS). Two other objects (B215 and B289) show an IR excess longwards of 2.3 μm but do not show emission lines nor NIR excess in the X-shooter spectral range. This is a sign that they are surrounded by dusty disks. Six of the observed objects (B163, B243, B268, B275, B331, and B337) have H I, O I, and Ca II double-peaked emission lines as well as CO bandhead emission and a NIR excess in their SED; this indicates that they are surrounded by a gaseous disk with a dust component. RT17 conclude that they are MYSOs with disks that are probably a remnant of the assembly process on their way to become main-sequence stars with masses between 6 and 15 M⊙ consistent with having undergone high mass-accretion rates (Macc ~ 10−4−10−3 M⊙ yr−1), hence with formation timescales of ~105 yr.
2.1 Target objects, observations, data reduction, and contaminants
In this paper, we focus on modelling the circumstellar disks of two MYSOs in M17: B243 and B331, whose stellar properties are summarized in Table 1 (RT17). These two targets show pronounced disk features, namely strong and clearly detected double-peaked hydrogen emission lines up to high order in the Paschen series and an infrared excess. B331 also shows Brackett series lines. Their X-shooter spectra cover the optical to NIR (300–2500 nm) with a spectral resolution of 5100 in the UVB arm, 8800 in the VIS arm, and 11 300 in the NIR arm. The observed spectra are flux calibrated using spectrophotometric standards from the ESO database.
The H I circumstellar emission in the spectra of the two sources can be contaminated with various other features, including nebular emission, stellar absorption, circumstellar emission by metallic species, and interstellar absorption. Some of these can and need to be corrected for in order to isolate the hydrogen emission from the disk itself.
The spectrum of B331 shows strong nebular emission, which contaminates the hydrogen circumstellar emission by adding a strong central component. Most of this emission was corrected for using multiple Gaussians fitted to off-source nebular emission as described by van Gelder et al. (2020). The nebular emission in B243 was not as substantial, and did not require further correction. Additionally, the studied hydrogen lines are blended with other emission lines originating from the disk, such as weak Ca II λ8498, λ8542, λ8662 Å emission in B243, blending the blue wings and peaks of Pa-16, -15, and -13. It is beyond the scope of this study to investigate the emission mechanism of these blending lines. Consequently, we cannot accurately predict their strengths and we refrain from correcting for their presence.
The H I cirsumstellar emission lines are superimposed on the stellar hydrogen absorption lines. The measured flux needs to be corrected for this absorption to isolate the disk contribution. The applied procedure to do so is described in Sect. 2.3.
Finally the spectrum also shows relatively wide absorption features originating from the interstellar medium. These so-called diffuse interstellar bands (DIBs) are likely due to large carbonaceous molecules such as (e.g. Sarre 2006; Cox et al. 2017; but see Campbell et al. 2015). An example of such a feature is seen in Fig. 1 at λ8630 Å. We exclude regions with known DIBs when estimating the continuum level.
Slight variations in the continuum of the normalized spectrum are still visible after applying the correction for stellar and nebular emission. These variations, due to unknown features, irregularities in the response curve, flux calibration, and uncertainties in the correction and normalization method, are the main uncertainty on the data. We conservatively estimate these uncertainties to be approximately 5% of the continuum on average and assume this to be normally distributed.
Stellar properties and extinction parameters of B243 and B331 updated from RT17.
![]() |
Fig. 1 Illustration of the normalization method. Top panel: the flux calibrated observed spectrum of B243 in black. The photospheric model spectrum is indicated by the orange dashed line and the continuum of this model by the green dash-dotted line. The model spectrum is matched to the observed spectrum using a least squares method in the unhatched wavelength ranges. At the top the location of the Paschen lines and a DIB are indicated. Bottom panel: the photosphere corrected and normalized observed spectrum (in solid black) and photospheric model (in dashed blue). The horizontal blue line indicates unity. |
2.2 Photometric data
The photometric data used in this work are taken from the literature and catalogues. For both stars we included data from DENIS (Fouqué et al. 2000), 2MASS (Skrutskie et al. 2006), and Spitzer GLIMPSE (Reach et al. 2005). WISE photometric data were included for B243 (Wright et al. 2010), as well as points from Kassis et al. (2002) for B331. SOFIA photometry of B331 at 20 and 37 μm Lim et al. (2020) show very high flux and, due to SOFIA's low spatial resolution, likely contains emission from outer material not belonging to the circumstellar disk. Therefore, it was not included in the analysis. Additionally, the spatial distribution of the SOFIA data does not agree with the 2MASS and Spitzer data, with the SOFIA data being significantly offset with respect to the others.
2.3 Normalization and correction of stellar features
To correctly model the disk emission lines contaminant stellar absorption features had to be removed, that is the flux absorbed in the stellar atmosphere had to be added back to the observed flux. To do so, a good estimate of the continuum is essential. However, in and around the hydrogen line series it is not possible to accurately make this estimate based on the observed spectra alone, as the Paschen line wings blend to form a 'secondary' continuum flux almost everywhere in the studied part of the spectrum; see the top panel of Fig. 1. As the stellar properties of B243 and B331 were already determined by RT17, see Table 1, we chose to simultaneously normalize and correct for the stellar features by matching model spectra to the flux calibrated observed spectra, for which we used a BT-NextGen atmospheric model computed with the PHOENIX code (Hauschildt et al. 1999; Allard et al. 2012).
The model spectrum, with Teff = 13 000 K and log g = 4.0, is rotationally broadened and convolved with the spectral resolution of the observation. Then the flux is scaled to match the observed spectrum using a least squares algorithm only taking wavelength ranges into account free from circumstellar or interstellar features (non-hatched regions in Fig. 1). As the calibrated spectrum has slight flux variations this scaling is done piece wise to find the best agreement between model and observation. The thus obtained scaled model is continuum subtracted to obtain the stellar features, which are then added to the observed spectrum. The resulting spectrum now only contains stellar continuum, disk features, and interstellar features. This spectrum is then divided by the stellar continuum from the model to get a normalized spectrum. Figure 1 shows an example of this process in part of the Paschen series regime of B331. A more extended part of the spectrum of both B243 and B331 is shown in Fig. A.1. The uncertainty on the stellar parameters will affect this normalization. Overestimating the surface gravity of the star results in wider hydrogen lines, which would result in stronger wings in the circumstellar emission profile. Underestimating the gravity would result in absorption features in the wings. The uncertainty in the temperature slightly affects the strength of the hydrogen lines.
3 Model
The double peaked hydrogen emission lines observed in B243 and B331 suggests that the emission originates in a circumstellar disk. Neither of the two objects show strong indication of outflows present in the X-shooter slit or other images. Therefore, we assume that the bulk of the emission originates from the central star and a circumstellar disk. A disk wind may be present, possibly contributing to the observed hydrogen emission, however, this is not modelled.
The disks are modelled using the radiation thermo-chemical code PRODIMO1 (Woitke et al. 2009, 2016; Kamp et al. 2010, 2017; Thi et al. 2011). PRODIMO self-consistently solves the chemistry, gas and dust heating and cooling, as well as line and continuum radiative transfer. The disk is assumed to be azimuthally symmetric and in steady state, and to follow a Keplerian rotation profile.
The applied chemical network accounts for 13 elements, namely H, He, C, N, O, Ne, Na, Mg, Si, S, Ar, Fe, and Polycyclic Aromatic Hydrocarbons (PAHs). In the present work, the latter only contribute to the free electron density. A total of 100 chemical species can form and interact. For the elemental abundances we adopt the Proto-Sun based values of Lodders (2003). These abundances are not depleted to compensate for the metals that are trapped in dust grains, as the chemistry we focus on takes place in the dust free inner region of the disk. The chemical reactions and rate coefficients are taken from the UMIST 2006 database (Woodall et al. 2007).
The dust temperature and continuum emission are solved for the given stellar and interstellar radiation field and dust opacities. The adopted dust composition is a mix of olivine (Mg2SiO4), amorphous carbon, and vacuum. Dust settling is taken into account following Dubrulle et al. (1995) adopting a turbulent viscosity parameters α = 10−3. A complete description of the treatment of solids can be found in Woitke et al. (2016).
Next, the gas heating and cooling are solved simultaneously with the chemistry as these are coupled. The gas temperature is based on 100 heating and 93 cooling processes implemented in PRODIMO revision number 3053.
Because the (free) parameters describing the density structure in the disk are found to be most influential on the studied lines (see Sect. 3.2), we briefly describe their relations. The relevant parameters are the total dust and gas mass Mdisk; inner and outer radius of the disk, Rin, and Rout, respectively; radial column density exponent ϵ; reference vertical scale height H0 specified at a reference radius R0, and scale height flaring power index β. The inner and outer radii specify the extent of the disk. The radial profile of the vertically integrated column density, Σ, of the disk is given by
(1)
where Σ0 is determined such that the integration of this profile over the radial extent of the disk results in Mdisk. The vertical density structure at a given radius r is defined by
(2)
and using Eq. (1)
(3)
is solved to find the midplane density ρ0(r). The vertical scale height H is a function of the distance to the central star, that is
(4)
where, H0 is the scale height at the reference radius, R0, and β is the flaring power.
The disk is in Keplerian rotation. It is irradiated by a central star of mass M★, effective temperature Teff, and luminosity L★. Additionally, an X-ray luminosity is included referred to as LX which is plays an important role in the chemistry and thermal balance of the gas. A mass accretion rate Ṁ can be specified as well. This causes viscous heating in the disk. For B243 we include a minor mass accretion, however, the heating resulting from this is minor, <1%, compared to that of the radiation field. The accretion rate was chosen such that the disk would have a reasonable lifetime given its mass. In B331 no effect was seen as result of similar accretion rates, therefore no accretion heating was implemented.
We divide the disk in two regions, referred to as inner disk and outer disk. The inner disk is located closest to the star and is purely gaseous. This zone has a negligible dust mass as temperatures are too high for solid particles to exist for any appreciable amount of time. The hydrogen line emission originates in this part of the disk. The outer disk also contains dust, hence determines the infrared continuum emission. The inner radius of the outer disk results from a fit to the SED. For the outer disk the model parameters are given the label '2'; for instance, the total mass of the outer disk is M2,disk.
3.1 Model modifications
PRODIMO is originally designed to model T Tauri and Herbig disks, that is planet forming disks around cooler less-luminous stars, where a photospheric EUV radiation field is not expected to significantly affect disk processes. Therefore, hydrogen-ionizing photons with energies between 13.6 and 100 eV were previously ignored. For hotter stars the EUV contribution becomes stronger, and needs to be considered. PRODIMO does include X-ray radiation (Aresu et al. 2011; Rab et al. 2018), defined as photons with energies >100 eV, as the ionizing properties of such photons are essential for the chemistry in the disk (Glassgold et al. 1997). In order to properly treat hydrogen ionization, we extended the code to also include the EUV radiation field.
The standard atomic model for hydrogen in PRODIMO considers 25 (n, ℓ)2 energy levels up to n = 5. This does not allow for a prediction of the observed Paschen and Brackett series, reaching up to n = 16. Therefore, the model atom was replaced by a new model using the NIST atomic database (Kramida 2019) to include the energy levels up to n = 20 and the transitions between these states. The collisional excitation cross sections are calculated using Jefferies (1968) as
(5)
with E0 the energy of the transition between energy levels u and l, k the Boltzmann constant, T the temperature of the gas, flu the oscillator strength of the transition, and g1 and gu the statistical weights of the lower and upper levels.
For the treatment of photo-ionization, see for example Mihalas (1978), we use
(6)
with v the frequency of a given photon, the charge of the ion Z = 1 for hydrogen, gII the bound-free Gaunt factor, which we assume to be 1, and ν0(n) the ionization frequency for energy level n.
![]() |
Fig. 2 Line profile of Pa-14 from the model of B243 for various parameter values. In each of the panels one parameter is varied, while the others remain fixed at an inner radius of 0.06 AU, a reference scale height of 0.5 AU, a total disk mass of 10−3 M⊙, and an inclination of 75°. |
3.2 Model fitting and model grid
We fitted the emission lines and SED separately, by varying certain parameters of the inner and outer disk, respectively. In the inner disk the hydrogen line profiles were fitted using an interpolated 4-dimensional grid of models. The free parameters of this grid are inner radius (Rin), reference scale height (H0), inner disk mass (Minner), and inclination (i). These parameters have a strong effect on the resulting line profiles. This is illustrated in Fig. 2, where large variations in line strength and shape are visible. Increasing the inner radius of the disk results in weaker and narrower lines, as the disk moves further from the star where it is less illuminated and orbital velocities are lower. The smallest sampled scale heights result in weaker and wider lines. As the material is confined closer to the mid-plane of the disk, less of the stellar radiation reaches the disk; moreover, it is harder for this radiation to reach more (radially) distant parts of the disk due to the high density. Consequently, line emission emanates from regions with relatively high Keplerian velocities. As the volume density of the disk decreases with increasing scale height and light reaches larger radial distances, the peak separation decreases. The line strength initially increases, because more radiation is intercepted by the disk. However, at the highest scale heights, the decrease in volume density starts to outweigh the increased illumination and the line strength decreases again. Increasing the mass of the disk results in a higher density and thus a stronger emission line, until it saturates from which point on the line strength will decrease. Changing the inclination angle changes the projected velocity distribution of the emission.
Other parameters such as the flaring power, β, and column density exponent, ϵ, are degenerate with the varied parameters in how they influence the line profiles. This degeneracy results from the small area of the disk (see Sect. 4.4) from which the H I emission lines originate. Therefore, the surface density exponent and disk mass affect the density of the emission region in a similar way. The flaring power affects the scale height of this region in a similar way as the reference scale height. A fiducial flaring power (of e.g. 1.2) results in an extremely large vertical extent in the outer parts of the inner disk as a relatively large scale height is required close to the star. To prevent this, we adopt a value of β = 0. 5 for the inner disk.
The outer disk is modelled using a similar grid, varying the outer radius, scale height, and mass. An example of the effects of the parameters on the SED is shown in Fig. 3. The inclination has a modest effect on the SED. For consistency with the inner gaseous disk we set this parameter to the best fit value of the line profiles. The scale height has a strong effect on the continuum disk emission. For the relatively massive disk shown here, a larger outer radius allows for cooler dust causing stronger mid-infrared emission. However, for a relatively low disk mass (not shown here) the redistribution of mass that results from assuming a larger outer radius causes the inner regions of the disk to become optically thin, reducing the near-infrared flux. The mass of the disk affects the full spectral range. For a higher mass the emission increases until the disk becomes optically thick. This happens first at shorter wavelengths.
The parameter values of the calculated grid points are listed in Table 2. The range of adopted values allow for both much weaker and much stronger lines than observed. As the range can cover several orders of magnitude we opt for a relatively coarse grid and interpolate between the models using an N-dimensional linear interpolation algorithm (Virtanen et al. 2020). With this interpolation a significantly higher resolution grid consisting of 3 125 000 points is calculated. This finer grid is used to fit the data. The interpolated data are verified to be accurate by calculating full models at interpolated grid points. Typically, the interpolated line profiles diverge by <1% from full model calculations, however in some cases the stronger lines can differ by up to 10%.
A summary of input parameters for PRODIMO is shown in Table 3. Parameters not listed are identical to the PRODIMO standard parameter values obtained in the DIANA project (Woitke et al. 2016). Full input files can be found online3. Parameter values listed in the table that are held fixed are the result of exploratory fitting. The fixed parameters are either degenerate with other parameters (see Sect. 3.2) or do not impact the fitting strongly, therefore should not be considered well constrained. Notably, the outer radius of the inner disk is poorly determined as the H I lines form within the first one to two AU of the inner disk, see Sect. 4.4.
The hydrogen emission lines studied in this work are for B243: Pa-β through P-16, save for Pa-8 and Pa-10 which are in regions of strong telluric absorption. The same lines are studied for B331 as well as Br-γ, Br-10, Br-11 Br-12, Br-14, and Br-16. The disk parameters are determined by fitting a subset of the lines simultaneously. The subset aims to include the lines most representative of the inner region of the disk. This subset consists of Pa-16 to Pa-9 for B243 and for B331 Br-16 to Br-11 is added to that set. Additionally, each line if fitted separately allowing the inner radius, scale height, and mass to vary. The inclination is fixed to the best fit value of the combination of lines.
![]() |
Fig. 3 Models of the SED for B243 for various parameter values. In each of the panels one parameter is varied, while the others are held fixed at an outer radius of 10 AU, a reference scale height of 0.04 AU, a disk mass of 10−4 M⊙, and an inclination of 15°. |
Overview of the model grid used to analyse B243 and B331.
4 Results
We present the results of our fitting efforts of the photometric data and X-shooter spectra using PRODIMO. First, we discuss the source of the hydrogen ionization in Sect. 4.1. Subsequently, we describe the fitting results for B243 and B331 in Sects. 4.2 and 4.3. Correlations between parameters are briefly presented in Appendix C.1. Finally, we present the derived physical properties of the modelled disks in Sect. 4.4. All results are also available online3.
Only the disk contribution to the emission is fitted. For the stellar contribution we use the best fit properties as listed in Table 1. We report the scale heights at the inner radius of the disk rather than at the reference radius, which is used to set up the grid of models. The mass is translated to the column density at the inner radius of the disk, as the latter is better constrained. These quantities are calculated using the equations in Sect. 3. An overview of the best fit parameters is shown in Table 4.
Overview of the input parameters for both the central star and the disk.
Best fit model parameters and their 1∈ uncertainties.
4.1 Source of hydrogen line emission
The spectra of both B243 and B331 show significant hydrogen line emission. The mechanism through which these lines form is recombination, as the temperatures in these disks are not high enough to allow collisional excitation of the relevant levels. Photo-ionization by the central stars cannot account for the hydrogen ionization rates (see Sect. 4.4) required to produce the observed line strengths – the stars are simply too cool (Teff being 13 500 and 13 000 K). Therefore, an alternative source of ionization in the inner part of the disk is required. In the simulations we find charge exchange reactions between ionized sulfur (S+) and hydrogen, that is
(7)
to be the dominant source of hydrogen ionization. Sulfur is ionized by UV photons with energies of at least 10.3 eV and then collides with neutral hydrogen, ionizing it. Other charge exchange reactions with hydrogen also take place, however these do not contribute to the total hydrogen ionization significantly.
4.2 B243
4.2.1 Hydrogen emission lines
All of the Paschen lines in B243 show clear double peaked emission consistent with being formed in a rotating disk; see Fig. 4. Higher Pa-series lines become weaker and display more pronounced wings and a larger velocity separation of the peaks, see Fig. 12 and RT17. The range in line strength is large, Pa-β peaking about three times over continuum while for Pa-16 this is only about 10%.
The blue lines in Fig. 4 show fits to each line individually. These match the observed profiles very well, save for the blue wing and blue peak of Paschen 16, 15, and 13, which are contaminated with superimposed Ca II triplet emission (see Sect. 2.1). The corresponding model parameters are shown in orange in Fig. 5 as function of oscillator strength (bottom axis) and line identifier (top axis). The parameters derived from the Paschen lines show fair agreement with one another though we note that for this star the uncertainties are sizable, in part due to some level of degeneracy in several of the model parameters. For instance, the inner radius Rin and scale height H both impact the line broadening – increasing Rin leads to narrower lines requiring a lower scale height and higher mass (see Fig. 2). We find that the Paschen lines indicate a small inner radius, that is the disk extends to the stellar surface to within the order of a stellar radius (R★ = 0.035 AU). The surface density at the inner radius is at least ~100g cm−2. The vertical scale height of the disk at its inner radius is ≳2 R★, which is significantly larger than implied by vertical hydrostatic equilibrium. Such a puffed-up inner zone is quite common in HAeBe stars (e.g. Dullemond et al. 2001) and may originate from its direct exposure to stellar light (see also Sect. 5.2). The radial extent of the inner disk cannot be constrained from the hydrogen lines, hence, such constraints – if any – should come from the dust modelling.
The horizontal dashed line in Fig. 5 shows the best fit parameters when analysing Pa-9 and Pa-11 to Pa-16 simultaneously. The corresponding line profiles are shown in orange in Fig. 4. We find a best fit inclination of 75°, which is used for the individual line fits. The best value for the column density at the inner radius of the disk lies lower than the confidence intervals of some of the individual fits. This is due to a local minimum in the χ2. The combined fit recovers all line strengths – that overall vary by a factor of ~30 – to within a factor of two. It does reveal a systematic trend in that the model under-predicts the higher (i.e. weaker) Paschen lines and over-predicts the lower Paschen (i.e. stronger) lines. We discuss this discrepancy in Sect. 5.
![]() |
Fig. 4 Best fitting models in blue to the observed lines in black for B243. The blue shaded regions indicate models that fall within the 1σ confidence interval. Each line is fitted individually with the inclination set to 75°, the respective parameters are shown in Fig. 5. The orange profiles indicate the best fit when fitting Pa-16 to Pa-9 simultaneously. Notice the factor ~30 range in line strengths covered by the Paschen series. |
![]() |
Fig. 5 Parameter values of the inner disk with their uncertainties for each modelled line of B243 and B331 as function of their oscillator strength. The horizontal dashed lines indicate the best fit parameters for B243 when fitting Pa-9 and Pa-11 to Pa-16 simultaneously and for B331 when fitting Pa-9, Pa-11 to Pa-16 and Br-11 to Br-16. The corresponding line profiles are shown in Figs. 4 and 7 for B243 and B331, respectively. For each fit an inclination of 75° is used. |
![]() |
Fig. 6 Observed and best fit SED of B243 and B331. For each star, the scatter points indicate the dereddened observed photometry and the black dashed lines the Kurucz model used as input for the stellar spectrum in PRODIMO. The best fitting models are drawn in green, with the shaded region indicating the 1σ uncertainty. For B331 the photometry from Lim et al. (2020) is indicated as upper limits; these are not considered in the fitting procedure. |
4.2.2 SED
The NIR excess of B243 is well fitted with dust particles emitting at an almost constant temperature of ~1500 K; see Fig. 6. This temperature is about the condensation temperature of silicate-based and carbonaceous grains placing the inner radius of the dust disk at 3 AU distance from the star. We use the best fit inclination of the inner gaseous disk of 75°. We find that the outer radius of the disk is poorly constrained; the best fit lies at ~6.5 AU from the central star. The disk may appear truncated due to some self-shadowing. Longer wavelength photometry is needed to better investigate the extend of the disk. A very small truncated dust disk, containing only a small amount of hot dust. The scale height at the inner rim of the dusty outer disk is R★, therefore the dust is more confined to the disk mid-plane than is the gaseous material at the inner edge of the inner disk, consistent with the inner most parts of the disk being puffed up. We find a lower limit of the surface density at this inner rim of the dusty disk of ~9 g cm−2, the material becoming optically thick and insensitive to further increases of the surface density.
In summary, the gas and dust analysis points to B243 having a disk with a best fit outer radius of 6.5 AU of which the inner part, up to 3 AU, is dust free. The mass contained in this small disk is at least ~2 × 10−5 Μ⊙ when adopting a surface density exponent ϵ = −1 for both the inner and outer disk, or ~3 × 10−4 Μ⊙ when assuming ϵ = 0 for the inner disk. As we lack photometric and spectroscopic information at wavelengths longwards of 10 μm we cannot exclude the presence of a (tenuous) cold disk further out, but a gap or self-shadowed region is should be present based on the best fit model.
4.3 B331
4.3.1 Hydrogen emission lines
Most of the hydrogen emission line profiles are clearly double peaked and show relatively broad line wings. The latter correspond to high rotational velocities of up to ~200 km s−1. This suggests a Keplerian rotating disk extending down to the stellar surface. We find a best fitting inner radius of <0.12 AU (<1.2 R★) for nearly all lines. The scale height at the inner rim is quite similar as for B243 and points to the inner disk being puffed up. The surface density at the inner radius is well constrained at 1.1 × 103 g cm−2. A similar column density has been found for NGC 2021 IRS based on CO emission (Gravity Collaboration 2020).
For B331 too the observed Paschen and Bracket series span a large range in line strengths, of about a factor of 20 (Fig. 7). The dashed blue lines in Fig. 5 show the best fit parameters taking into account the same set of Paschen lines as for B243 and, in addition, Br-11 to Br-16. With this set, the line strengths are again reproduced to within a factor of two save for Paß and Brγ for which the difference is somewhat larger; see the orange line in Fig. 7. None of the models within the confidence interval reach the upper limits. The best fit yields an inclination of about 75°.
Similar to B243 we find a systematic trend in that the observed line strengths vary less with Paschen and Bracket line number than do the model predictions. We return to this in Sect. 5.
4.3.2 SED
The SED of B331 shows an IR excess at λ > 3 μm; see Fig. 6. We use the best fit inclination of the inner disk of 75°. The PRODIMO dust model results in a maximum temperature of ~400 Κ positioning this warm dust at ~100 AU from the star. The outer radius of the dusty disk remains unconstrained, simply because we lack photometric or spectroscopic constraints at wavelengths well beyond 10 μm. The amount of 3–10 μm flux is a function of the scale height of the inner rim of the outer dusty disk. We find for this scale height AU. The surface density at this location should be at least 19 g cm−2, with improving fits for higher values up to at least ~700 g cm−2.
In summary, the overall disk properties of B331 are quite different from those of B243. B331 lacks hot dust; warm dust seems only present starting at about 100 AU. This suggests that for this source a disk gap is created. The inner boundary of this gap, if present, is poorly constrained (the hydrogen spectral lines indicate it should start at a distance >2 AU). The properties of the inner gaseous disk up to 3.5 AU appear quite similar to B243, both having a puffed-up inner rim.
4.4 2D Disk structure of B243 and B331
Figures 8 and 9 show the hydrogen ionization fraction, density, and temperature structure of the best fitting inner disk models for B243 and B331, respectively. The distance axis is linear in the case of B243 and logarithmic in case of B331. The hydrogen ionization fraction throughout the disks is low, even at the inner rim of the gaseous disk (where it reaches ≲10−3). This is essentially due to the fairly low effective temperatures of the two stars. The ionizing photons do not penetrate far into the disk as the mostly neutral (hydrogen) medium becomes quickly optically thick for Lyman continuum radiation. B331 is almost an order of magnitude more luminous than B243, which is why for this star the thin somewhat ionized shell reaches deeper disk layers. For the same reason, the disk of B331 is heated to higher temperatures further out.
The black contours in the figures indicate the origin of half of the Pa-16 line flux. The formation region stretches out to about a radial distance somewhat less than one AU, that is hydrogen lines probe the very inner part of the disk only – dominated by the puffed up inner zone. Similar to the work by Koumpia et al. (2021), Caratti o Garatti et al. (2016), and Kraus et al. (2010), we find the hydrogen emission to originate from a region closer to the central star than the continuum emission, based on the location of the inner radius of the dust disk. We find the extent of the hydrogen emission region to reach up to ~ 1 AU, which is in line with the most compact emitting regions found in interferometric studies of Koumpia et al. (2021) and Caratti o Garatti et al. (2016), but is more compact than what they typically derive.
![]() |
Fig. 7 Best fitting models in blue to the observed lines in black for B331. The blue shaded regions indicate models that fall within the 1 σ confidence interval. Each line is Atted individually with the inclination Axed at 75°, the respective parameters are shown in Fig. 5. The orange proAles indicate the best At when Atting Pa-16 to Pa-9 and Br-16 to Br-11 simultaneously. Notice the factor ~20 range in line strengths covered by the Paschen series. |
![]() |
Fig. 8 Hydrogen ionization, density, and temperature structure of the inner region of B243. The parameters used match the values indicated by the full red line in Fig. 5. The black contour indicates the approximate origin of 50% of the Pa-16 line emission based on vertical escape probabilities. The orange circle on the left indicates the central object. |
5 Discussion
Before summarizing and discussing the disk properties of the M17 members B243 and B331, we show their positions in the Hertzsprung-Russell diagram (HRD) in Fig. 10. RT17 estimate an age for the M17 star-forming region of less than 1 Myr. This, combined with both sources showing both stellar absorption and disk emission features, strongly suggests that the stars are premain-sequence objects on Thomson or Henyey (a.k.a. radiative) tracks; in the final phase of their formation contracting towards the zero-age main sequence. The luminosity (and associated mass of ~12 M⊙) classifies B331 as an MYSO source. The estimated mass (of ~6 M⊙) of the lower luminosity source B243 is more representative of the high-mass end of Herbig Be sources. The total MYSO lifetime, defined as the time from passing the birthline until arrival on the main sequence, is less than 105 yr for B331 and a few times 105 yr for B243 (Hosokawa & Omukai 2009); both stars have progressed considerably in this evolution. Given these short evolutionary time scales and the observation that higher mass main-sequence stars in M17 lack disk signatures (RT17), one may anticipate that the disks of B331 and B243 are in the process of being cleared.
A schematic diagram illustrating the main properties derived for the two disks is shown in Fig. 11. Both stars feature a gaseous disk that almost (within a stellar radius or ~0.1 AU) reaches the stellar surface. The hydrogen emission lines originate from within the first AU of these disks, signifying that the full extend and properties of the gaseous disk are not probed by these diagnostics. The disk of B243 contains hot dust of ~1500 K, whereas the hottest grains surrounding B331 have a much lower temperature of 400 Κ implying a significant dust free inner zone spanning ~100 AU. The extend of the dust free inner zone depends on the continuum optical depth of the gaseous disk. A very opaque gas disk would move the dust disk closer to the star. The best fit model for the dust disk of B243 suggests an outer radius of ~6.5 AU, however, this is poorly constrained. The SED modelling suggests B243 is consistent with having a Group II disk and B331 a Group I disk. However, as our diagnostics are not sensitive to cool dust (our longest wavelength point is at 8 μm), we cannot rule out the presence of a (large) gap and cool dusty outer disk or a self-shadowed region in B243. Without the 8 μm point the outer radius of the hot dust disk would not be constrained. We conclude that both disks appear disrupted and identify them as transitional disks. Though not subject to analysis, both disks show pronounced CO-bandhead emission (RT17, Poorta et al., in prep.).
The gaseous inner disks of both stars have significantly larger scale heights close to the star than would be expected from standard hydrostatic calculations, consistent with puffed-up inner rims (Dullemond & Dominik 2004). This is further discussed in Sect. 5.2. We adopt a Keplerian rotation profile and did not find any indication that parts of the disk regions probed by the Paschen lines experience significant sub- or super-Keplerian motion. Still, we cannot fully exclude deviations from Keplerian motion as the rotation velocity is somewhat degenerate with the inner radius and inclination of the disk.
The small gaseous disk inner radii (of less than 2 R★) differs from the situation in lower mass YSOs. In these lower mass counterparts the disk material does not directly reach the stellar surface, but gas loops towards the surface via magnetoshperic accretion (e.g. Ingleby et al. 2013). This leaves a gap between the disk and star. Fairlamb et al. (2015) study a large sample of Herbig AeBe stars and find the UV excesses of early type Be stars to be inconsistent with magnetospheric accretion (see also Wade et al. 2007). This suggests that higher mass stars lack strong magnetic fields. Another indication for the absence of a strong magnetic field is consistent with Keplerian rotation. Strong magnetic fields can cause the disk to co-rotate with the stellar surface. This results in lower than Keplerian velocities close to the star and would likely result in too low velocities to reproduce the observed line profiles. In conclusion, the disks reaching almost the stellar surface and being compatible with Keplerian rotation close to the inner rim is in line with a lack of magnetospheric accretion, consistent with the general consensus that most higher mass stars lack a strongly magnetized envelope (Mottram et al. 2007; Oudmaijer et al. 2017).
![]() |
Fig. 9 Same as Fig. 8, but now for B331 on a logarithmic spacial scale. The parameters used match the values indicated by the dashed red line in Fig. 5. |
![]() |
Fig. 10 Adapted from Fig. 9 of RT17. Positions of B243 and B331 in the Hertzsprung-Russel Diagram. The solid blue line on the left indicates the zero age main sequence. The dashed red line on the right indicates the birth line. The solid line connecting the birth line and the zero age main sequence indicate the MIST pre-main sequence evolutionary tracks of stars with masses indicated on the right (Dotter 2016). The colour of this line indicates the relative time between birth and reaching the main sequence. The dash dotted grey lines indicate the isochrones corresponding to 3 × 103, 104, 3 × 104, 105 and 3 × 105 yr from top to bottom respectively. |
![]() |
Fig. 11 Schematic overview of the geometric properties of the disks of B243 and B331. |
5.1 Uncertainties in stellar properties
Table 1 shows the stellar parameters of B243 and B331. For B243 they result from quantitative spectroscopy; for B331 from SED fitting. The main uncertainties on the properties of these stars are the temperature, extinction (AV and RV), and luminosity. The extinction parameters do not affect this work significantly as we fitted the normalized line profiles, and the extinction at the NIR is modest. The luminosity affects the dust sublimation radius of B243 and location of warm dust for B331. A higher luminosity also results in higher hydrogen line flux (but not necessarily a higher normalized line flux, as the continuum flux is also increased).
The surface gravity of B331 is unconstrained; we assume a value log g = 4.0. The gravity of the star mostly affects the wings of the circumstellar H I emission profiles through the normalization process (see Sect. 2.3), a too high gravity giving a broader emission profile. This would have to be compensated by a lower scale height or smaller inner radius.
5.2 Puffed up inner rim
At the inner rim, both B243 and B331 have a scale height significantly larger than expected based on the hydrostatic equilibrium approximation H ~ cs/ΩK, with cs the sound speed and ΩK the Keplerian rotation frequency. For T = 6000 K, r = 0.2 AU and M = 6 M⊙ the adopted scale height for B243 is 15 times the hydrostatic value. This 'puffed' up region is unlikely to extend far in to the disk. Despite such puffed up regions have been found before in efforts to model the inner region of disks (e.g. Woitke et al. 2009; Dullemond et al. 2001). A possible explanation for the large scale height could be that a disk wind contributes to the observed emission. However, a disk wind would possibly have a velocity profile distinct from a puffed up Keplerian disk. However, the velocity may still be dominated by orbital motion if the emission originates from close to the disk, with only a slight broadening of the profile.
We also investigated the effect of mass accretion on the inner rim structure. To this end, we performed test calculations of the best fit models to both B243 and B331 in which we use an accretion rate of 10−4 M⊙ yr−1. In these cases the line emission gets significantly stronger with 5 and 1.5 times more line flux for B243 and B331, respectively. The increase in line strength could be compensated by lowering the scale height. This would require a significant gas reservoir to be present to feed these accretion flows.
5.3 Origins of the hydrogen line emission and the role of chemistry
The hottest parts of the disks of B243 and B331 reach temperatures of ~7000 K and hydrogen ionization fractions of ~10−3. We find charge exchange reactions to be the main source of hydrogen ionization. These reactions are facilitated by various heavier elements, in particular sulfur. Therefore, the abundance of these metallic species are of importance to the line formation. A lower metallicity will result in weaker lines. We note that the reaction rate coefficient of the charge exchange reaction with sulfur is highly uncertain. In PRODIMO an estimate for the reaction rate coefficient of 5 × 10−12 cm3 s−1 is supplied. However, Butler & Dalgarno (1980) find an upper limit of 3 × 10−15 cm3 s−1 at a temperature of 104 K for the reaction rate coefficient. The hydrogen emission line strength scales nearly linearly with the reaction rate coefficient, with a reaction rate reduced by a factor of 50 resulting in 40 times weaker lines. This would significantly affect the derived properties of the disk, requiring substantially higher disk masses.
Presently, the central stars' temperatures are too low to produce sufficient ionizing photons. Therefore, we conclude that the stars appear to be in a phase of formation in which chemistry dominates the production of ionized hydrogen. We note that as the stars move closer to the zero-age main sequence (ZAMS) photo-ionization may take over as the main mechanism of hydrogen ionization, particularly so for B331 (which may reach a ZAMS temperature of about 30 000 K), if a close-in gaseous disk should prevail to such a late stage.
The key role of chemistry also highlights the importance of the chemical network and the corresponding reaction rates. In this work we use the UMIST 2006 reaction network. We note that employing the UMIST 2012 reaction network results in ~20% weaker hydrogen emission lines. The charge exchange reaction rate coefficients of hydrogen and sulfur are not included in these networks, but instead are supplied by a separate PRODIMO input. We also do not find a significant difference in hydrogen line emission between the fiducial large and small chemical networks available in PRODIMO (Kamp et al. 2017).
5.4 Origin of the hydrogen line emission and accretion
Hydrogen emission lines are generally considered to be tracers of accretion when observed in (massive) young stellar objects (e.g. Ilee et al. 2014; Fairlamb et al. 2017; Mendigutía et al. 2011). This raises questions as to the possibility and nature of accretion in the disks. Fairlamb et al. (2017) speak of a correlation of line luminosities with accretion rate, while also noting that accretion need not be the physical origin of the lines. Their empirical relation between hydrogen line strength and accretion rate implies an accretion rate of ~10−5 and ~10−3 M⊙ yr−1 for B243 and B331, respectively. If a gap of at least tens of AU indeed separates the inner disks from more distant disk material (if present at all), these rates would deplete the inner disks, up to a few AU, of the order of years. The very small chance of detecting the systems at the exact moment of inner disk accretion suggests that either the inner disk is being fed from a more extensive outer reservoir, or, that the hydrogen emission is not an accurate tracer of accretion in these systems.
To expand on the latter possibility, the empirical relation mentioned above is derived using calibrations based on magnetospheric accretion modelling. Following from our previous discussion of disk geometry and magnetic fields, as well as drawing from recent literature on accretion luminosities in Herbig AeBe stars (e.g. Wichittanakom et al. 2020; Grant et al. 2022; Mendigutía 2020) it is likely that higher mass YSOs (≳4 M⊙) accrete by a different mechanism. An obvious alternative is boundary-layer accretion, where disk material slows down in a transition layer between the disk and the star, releasing its kinetic energy in that boundary layer. Lynden-Bell & Pringle (1974) relate the mass accretion rate to a black body emission from an annulus with temperature TBL. They assume all orbital kinetic energy is dissipated through radiation from the boundary layer. However, this assumes the layer to be optically thick and the radiation from the central star heating up the material is ignored. Despite those assumptions, using the accretion rates above we find TBL to be larger than the temperature in the inner rim of the disk for both B243 and B331. A more detailed modelling of boundary-layer accretion including the stellar radiation field and more accurate heating and cooling processes would be required to accurately link the temperature of the inner disk and the accretion rate of the star. Currently, there is no model that relates (hydrogen) line emission to accretion rate in the BL regime.
5.5 Disk effect on the central star
The disks, given their small inner radii and the possibility of ongoing accretion, might still affect the properties of their central star on its way to the main-sequence. Following the arguments made in the previous section and based on the derived disk masses, it is likely that the central stars have accumulated the mass at which they will start central hydrogen fusion. However, the inner disk mass may contain a sizable amount of angular momentum. Here we estimate the effect of the disk on the surface rotational velocity of the star once it reaches the ZAMS.
Assuming conservation of angular momentum for the stars and no further accretion we can calculate the surface rotational velocity upon arrival on the ZAMS. On the basis of their assumed stellar masses, we adopt radii of 2.5 R⊙ and 4.2 R⊙ at the ZAMS for B243 and B331, respectively (Brott et al. 2011). The radial density structure of the stars is approximated using the solution of the Lane-Emden equation for a polytropic index of 1.5 in both the current state and at the ZAMS. We further assume solid body rotation. Taking an inclination of 75° converts the v sin i of 110 km s−1 into vrot = 114 km s−1 for the present-day state of B243. After contracting to the main sequence this results in a rotational velocity of 342 km s−1, which is ~0.50 vcrit, where vcrit is the critical rotation rate. If we adopt a column density at the inner radius of the disk of 1000 g cm−2 and a relatively flat column density exponent of ϵ = −0.5, the angular momentum of a 5 AU inner disk would be ~17% of that of the star. The current vrot has not been constrained for B331. Doing the same calculation for this star under the same assumptions and taking the current v sin i = 110 km s−1, the same value as for B243, yields a ZAMS spin velocity of ~592 km s−1, which is ~0.80 vcrit. The angular momentum of the disk is then ~7% of that of the star.
These estimates suggest that the inner disk is unlikely to still contribute strongly to the ZAMS rotational velocity of the central star; its value is essentially controlled by stellar contraction. Spin velocities of 0.5–0.9 vcrit upon arrival on the main sequence are in line with previous findings for stars in the same mass range (Huang et al. 2010). For stars with a weak or lacking magnetic connection to the disk (see above) gravitational torques are expected to limit the ZAMS spin velocity to about half critical (Lin et al. 2011). This may be an indication that the current projected spin rate of B331 is less than the adopted 110 km s−1.
5.6 Disk disruption scenarios
The leading mechanisms responsible for the dispersal of disks around young stars are photo-evaporation (e.g. Gorti & Hollenbach 2009; Owen et al. 2010), stellar or planetary companion formation (e.g. Müller & Kley 2013), and possibly stellar winds for the most massive YSOs (Bik et al. 2006). Given the relatively high temperatures of the MYSO sources studied here relative to pre-main sequence stars of lower mass, photo-evaporation may be a contender for disk dispersal. Owen et al. (2010) provide a simple scaling relation for the characteristic radius where thermal evaporation by EUV light starts, after which more inner parts of the disks are cleared on a timescale of a few times 104 yr assuming the disk dispersal time is only a weak function of stellar mass as suggested by Owen et al. (2010). For our sources this radius is at about 100 AU; it corresponds well to the size of the dust free gap or zone we observe in both sources. In this mechanism, the innermost parts of the disks (i.e. those probed by H lines) survive the longest (Alexander et al. 2006). So, photo-evaporation as a disk dispersal mechanisms seems congruent with the derived disk properties for both B243 and B331. We remark however that photo-evaporation models do not take into account a puffed-up inner disk, as found for both stars, that may extinct a sizeable amount of the ionizing radiation, limiting the EUV-illumination of further out regions and hence efficient disk dispersal. Gorti & Hollenbach (2009) study the effect of X-ray, EUV, and FUV radiation on the photo-evaporation of disks. They find the photospheric FUV radiation to be the dominant energy source driving mass-loss from the outer regions of the disk, and EUV radiation to only affect the inner regions of the disk. We expect the inner gaseous disk to be optically thin to FUV radiation and the central stars to be bright in the FUV, therefore it would be a potential energy source for photo-evaporation.
Ongoing photo-evaporation would result in observable spectral features such as forbidden oxygen emission and H2 emission (e.g. Gangi et al. 2020). Though we have not detected any H2 emission, B243 shows forbidden [O I] λ630 nm emission. Its behaviour is different on-source than off-source pointing at a possible circumstellar origin (Derkink et al., in prep.). The [O I] emission of B331 does not allow us to make this distinction.
Gaps in proto-planetary disks are commonly observed in the low mass counterparts of MYSOs and were first identified in Strom et al. (1989). Kim et al. (2013) collate 105 of these disks. They propose the gaps are likely due to companions orbiting the central object, and that an origin linked to grain growth, turbulence or photo-evaporation processes is less likely. Companions formation through gravitational instabilities in efficiently cooling disks appears feasible on the timescales of MYSO formation, at least theoretically (e.g. Rice et al. 2003; Oliva & Kuiper 2020). In this scenario the companion will create a gap in the disk, preventing material to flow through its orbit efficiently. A small amount of gas may travel past the gap to sustain the observed gaseous inner disk while blocking dust, creating a dust free inner disk (e.g. Lubow & D'Angelo 2006). This scenario too seems in line with the derived disk properties for both stars. We note that such companion forming disks are expected to develop a spiral structure, which, if indeed the correct scenario, may perhaps play a role in explaining the (modest) discrepancies between predicted and observed H line profiles discussed below.
We conclude that both disk dispersal mechanisms appear reasonable candidates to explain the observed disk structures in B243 and B331.
![]() |
Fig. 12 Orbital velocity of the hydrogen gas based on the peak separation of the observed (circles) and simulated (stars) emission lines as function of their oscillator strength. Green markers indicate Paschen lines and yellow markers Brackett lines. The velocity is determined by fitting double gaussians to the profile. The simulated profiles are from Figs. 4 and 7. |
5.7 Deviating trends in hydrogen series properties
Figures 4 and 7 show that a single model cannot fully recover the line strengths of all modelled Paschen and Brackett series lines: lower (higher) series lines are predicted stronger (weaker) than observed. Differences remain within a factor of two though, while the range of line strengths spans a factor ~20–30. In addition to this, we observe a modest discrepancy in the trend of peak separation versus line strength.
Figure 12 shows that the disk velocity decreases for increasing oscillator strength for both B243 and B331. As in our Keplerian disk model all line emission originates from the optically thin inner disk region (i.e. from the first AU) no such trend is predicted. For classical Be stars a similar trends in found as in our systems (e.g. Kraus et al. 2012). These authors could spatially resolve the emitting region of different groups of hydrogen lines for one such Be disk and find spatially different formation regions. This indicates that in this Be disk the hydrogen lines are optically thick (unlike in our MYSO modelling).
Though we cannot pinpoint the cause of the mismatch in peak separation trend, part of it may possibly be explained by deviations from azimuthal symmetry in the inner disk. It also remains to be investigated how the inclusion of (boundary-layer) accretion and disk winds would affect the modelled lines. Massive YSO sources in M17, including the objects studied here, show variability in their spectra originating in their circumstellar disks (Derkink et al., in prep.). This variability includes excess line emission moving from the red to the blue part of the line on timescales consistent with the rotation period of the inner disk, underlining that likely there are processes active in the inner disks that break axial symmetry. Possible features include spirals, warps, large scale clumps or other, more chaotic, forms of asymmetries.
6 Summary
We have investigated the inner disk regions of two MYSO sources, B243 (with an estimated mass of 6 M) and B331 (12 M) in the star-forming region M17, using the thermochemical code PRODIMO. Likely, these inner disks are remnant structures from the star assembly process.
Our main diagnostics are (double-peaked) hydrogen lines of the Paschen and Brackett series and near-IR photometry. The first allow us to probe the kinematics and structure of the gaseous inner disk and the second the thermal emission of hot dust. The effective temperatures of the central stars are such that photo-ionization of hydrogen is much less important than ionization through charge exchange reactions, so it is chemistry that dominates the H line-formation.
Our main findings are:
A small puffed-up gaseous disk extends to very close to the stellar surface in both Group II source B243 and Group I source B331. In B243 a dust free cavity of 3 AU is present with hot dust of 1500 K at the inner rim of the dusty disk. We find a best fit outer radius of ~6.5 AU, but this is poorly constrained. The inner disk of B331 is dust free and probably of similar dimension, cool dust indicating an outer disk starting at about 100 AU with a dust free (and possibly gas free) gap in between.
The inner disk extending to almost the stellar surface suggests that some accretion might still be ongoing, likely through boundary-layer (BL) accretion. Magnetospheric accretion is less likely as the geometry of this accretion mechanisms suggests an inner gaseous disk gap. A lack of magnetospheric accretion is in line with the general consensus that a different accretion mechanism is at work in higher mass stars. A first order approach indicates that it would be interesting to investigate the BL mechanism for higher mass Herbig Be stars and MYSOs with detailed line modelling. The presence of disk winds remains an open question.
The inner disk contains too little mass and too modest angular momentum to significantly change the final (i.e. ZAMS) mass and final spin velocity of the stars. Concerning the latter, contraction towards the main sequence is the main effect impacting the spin velocity in the remainder of pre-main sequence evolution (see also Ramírez-Tannus et al. 2017). The angular momentum of the inner disk is of the order of 10% of that of the central star.
The disk structures of both sources strongly suggest that the disks are in the process of being cleared, that is they are transitional disks. Possible disk dispersal mechanisms are photo-evaporation and stellar or planetary companions formation. The derived properties of the disks are compatible with both scenarios, specifically the presence of a small inner disk and a disk gap towards a more distant (~100 AU) outer disk (in B331 and possibly in B243).
The study presented here does not fully characterize the disks orbiting B243 and B331; but it does present the first detailed 2D thermo-chemical radiative transfer modelling of hydrogen lines in such sources. High resolution, longer wavelength imaging with for example ALMA would greatly add to our insight in the properties of the outer regions of the disk. Whether such regions (if they exist at all in B243) would still play a role in the formation process of the central object remains to be seen, given the proximity of the stars to the ZAMS. The fate of such outer disks may simply be that they are dispersed by the concerted action of the star's H II region and stellar wind upon arrival on the ZAMS (e.g. Geen et al. 2021). Studies of line variability may contribute to our understanding of dynamical process in the inner disks, possibly due to the presence of companions. The combination of these different approaches may greatly help in unravelling the architecture of companion systems around massive stars, their possible migration (Ramírez-Tannus et al. 2021), and pre-main sequence or early main sequence merging with the primary star (Wang et al. 2022); topics that constitute new and exciting problems in massive star formation.
Acknowledgements
We thank the anonymous referee for insightful comments that helped to improve this manuscript. This publication is part of the project 'Massive stars in low-metallicity environments: the progenitors of massive black holes' with project number OND1362707 of the research TOP-programme, which is (partly) financed by the Dutch Research Council (NWO). This work is based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 089.C-0874(A). We thank SURF (www.surf.nl) for the support in using the Lisa Compute Cluster. This research has made use of NASA's Astrophysics Data System. This work makes use of the Python programming language (Python Software Foundation; https://www.python.org/), in particular packages including NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), and Matplotlib (Hunter 2007).
Appendix A Normalized and corrected spectra
The normalized and corrected spectra of B243 and B331 are shown in Fig. A.1. The spectra have only been normalized around the modelled lines. The continuum of B243 at wavelengths longer than 1500 nm deviates due to IR excess from dust emission.
![]() |
Fig. A.1 Normalized and corrected spectra of B243 (top panel) and B331 (bottom panel), focusing on the Paschen and Brackett series. The spectrum of B243 is not normalized at wavelengths longer than 1500 nm due to the normalization method's incompatibility with the IR excess. |
Appendix B Model test – HII region
The new hydrogen atom description and extended radiation field were tested by simulating an H II region and comparing the numerically calculated ionization edge with the analytical Strömgren radius. We simulate an O-type star with Τ = 38 000 Κ and L★ = 1.65 × 105 L⊙. We surround this star with large gaseous disk with a constant radial density in the mid-plane (ρ ~ 10−20 g cm−3). The disk does not contain dust or poly-cyclic aromatic hydrocarbons (PAHs). The disk material has (protocolar abundances. Fig. B.1 shows the ionization fraction in the simulated H II region. A clear ionization edge is visible. The radius at which n(H) = n(H+) matches the analytically calculated Strömgren radius, which is indicated by the dashed white line. The surface layer of the disk has a lower density, which moves the ionization edge to a larger distance from the star. The theoretical prediction of the Strömgren radius is based on the number of ionizing photons emitted and the recombination rate of electrons and protons which depends on the temperature and density of the circumstellar medium. We only consider case Β recombination with the coefficient
(B.1)
by Draine (2011), with . The radius at which the medium becomes neutral, the Strömgren radius, can then be described as
(B.2)
with nH the hydrogen particle density and Qion the rate at which hydrogen ionizing photons are produced. In this system Qion ~ 5 × 1048 s−1.
As the ionization edge of the model and the Strömgren theory agree, we conclude that the hydrogen photo-ionization and recombination balance is calculated correctly in the model.
![]() |
Fig. B.1 H II region simulated with PRODIMO. The colour indicates the ratio between neutral and ionized hydrogen. The white dashed line indicates the analytically calculated Strömgren radius based on the radiation field, temperature and density of the medium. |
Appendix C Correlations and grid behaviour
The parameters chosen to be varied have some degeneracy, which can result in large uncertainties in their best fit value. Additionally, the coarseness of the grid can result in what we will refer to as a 'noisy' χ2 statistic. This noise is the result of the correlation between parameters and misalignment in the step size of the grid and correlation strength. This can be seen in the top panels of each column of Fig. C.1 and Fig. C.2. When exploring a parameter, the optimal values of the remaining parameters can only be approximated to a limited degree. The optimal value is likely to lie between grid points. The size of this deviation varies for each point and results in a varying accuracy to which the minimum χ2 can be determined. This affects the accuracy of the optimal fit values and their uncertainties. In order to minimize this effect the coarse model grid is interpolated.
![]() |
Fig. C.1 Corner plot indicating the χ2 value for the model fit to the Pa-12 emission line of B243 as function of each of the parameters and each parameter combination. The dashed line in each top panel indicates the 1σ threshold. |
Appendix C.1 Correlations
We find correlations between the mass and scale height, inclination and scale height, and the inclination and mass as indicated in Fig. C.1 and Fig. C.2. An increase in mass can result in a similar line profile by increasing the scale height and inclination of the disk. This results in degeneracy allowing a relatively large range of parameter values to fit to the data. This degeneracy combined with the coarseness of the grid results in noisy χ2 values as function of a given parameter. When the parameter value moves to a different grid point, the optimal value of another parameter might lie between two grid points and will be missed. This results in a higher χ2 value. Therefore, the accuracy of the minimum χ2 for a given parameter is limited.
![]() |
Fig. C.2 Corner plot indicating the χ2 value for the model fit to the Pa-12 emission line of B331 as function of each of the parameters and each parameter combination. The dashed line in each top panel indicates the 1σ threshold. |
References
- Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006, MNRAS, 369, 229 [CrossRef] [Google Scholar]
- Allard, F., Homeier, D., & Freytag, B. 2012, Phil. Trans. R. Soc. London Ser. A, 370, 2765 [Google Scholar]
- Aresu, G., Kamp, I., Meijerink, R., et al. 2011, A&A, 526, A163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 [Google Scholar]
- Beltrán, M. T., & de Wit, W. J. 2016, A&ARv, 24, 6 [Google Scholar]
- Bik, A., Kaper, L., & Waters, L. B. F. M. 2006, A&A, 455, 561 [CrossRef] [EDP Sciences] [Google Scholar]
- Broos, P. S., Feigelson, E. D., Townsley, L. K., et al. 2007, ApJS, 169, 353 [Google Scholar]
- Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Butler, S. E., & Dalgarno, A. 1980, A&A, 85, 144 [NASA ADS] [Google Scholar]
- Campbell, E. K., Holz, M., Gerlich, D., & Maier, J. P. 2015, Nature, 523, 322 [NASA ADS] [CrossRef] [Google Scholar]
- Caratti o Garatti, A., Stecklum, B., Weigelt, G., et al. 2016, A&A, 589, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Charbonneau, P. 1995, ApJS, 101, 309 [Google Scholar]
- Chini, R., Elsaesser, H., & Neckel, T. 1980, A&A, 91, 186 [NASA ADS] [Google Scholar]
- Cox, N. L. J., Cami, J., Farhang, A., et al. 2017, A&A, 606, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Dotter, A. 2016, ApJS, 222, 8 [Google Scholar]
- Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton: Princeton University Press) [Google Scholar]
- Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [Google Scholar]
- Dullemond, C. P., & Dominik, C. 2004, A&A, 417, 159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957 [NASA ADS] [CrossRef] [Google Scholar]
- Fairlamb, J. R., Oudmaijer, R. D., Mendigutía, I., Ilee, J. D., & van den Ancker, M. E. 2015, MNRAS, 453, 976 [Google Scholar]
- Fairlamb, J. R., Oudmaijer, R. D., Mendigutia, I., Ilee, J. D., & van den Ancker, M. E. 2017, MNRAS, 464, 4721 [NASA ADS] [CrossRef] [Google Scholar]
- Fossati, L., Castro, N., Schöller, M., et al. 2015, A&A, 582, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Fouqué, P., Chevallier, L., Cohen, M., et al. 2000, A&AS, 141, 313 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Frost, A. J., Oudmaijer, R. D., de Wit, W. J., & Lumsden, S. L. 2021, A&A, 648, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Furlan, E., Luhman, K. L., Espaillat, C., et al. 2011, ApJS, 195, 3 [Google Scholar]
- Gangi, M., Nisini, B., Antoniucci, S., et al. 2020, A&A, 643, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Garufi, A., Meeus, G., Benisty, M., et al. 2017, A&A, 603, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Geen, S., Bieri, R., Rosdahl, J., & de Koter, A. 2021, MNRAS, 501, 1352 [Google Scholar]
- Glassgold, A. E., Najita, J., & Igea, J. 1997, ApJ, 480, 344 [Google Scholar]
- Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 1539 [Google Scholar]
- Grant, S. L., Espaillat, C. C., Brittain, S., Scott-Joseph, C., & Calvet, N. 2022, ApJ, 926, 229 [NASA ADS] [CrossRef] [Google Scholar]
- Gravity Collaboration (Caratti o Garatti, A., et al.) 2020, A&A, 635, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Hanson, M. M., Howarth, I. D., & Conti, P. S. 1997, ApJ, 489, 698 [NASA ADS] [CrossRef] [Google Scholar]
- Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
- Hauschildt, P. H., Allard, F., & Baron, E. 1999, ApJ, 512, 377 [Google Scholar]
- Hoffmeister, V. H., Chini, R., Scheyda, C. M., et al. 2008, ApJ, 686, 310 [NASA ADS] [CrossRef] [Google Scholar]
- Hosokawa, T., & Omukai, K. 2009, ApJ, 691, 823 [Google Scholar]
- Huang, W., Gies, D. R., & McSwain, M. V. 2010, ApJ, 722, 605 [Google Scholar]
- Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
- Ilee, J. D., Fairlamb, J., Oudmaijer, R. D., et al. 2014, MNRAS, 445, 3723 [Google Scholar]
- Ingleby, L., Calvet, N., Herczeg, G., et al. 2013, ApJ, 767, 112 [Google Scholar]
- Jefferies, J. T. 1968, Spectral Line Formation (Waltham: Blaisdell Publishing) [Google Scholar]
- Kamp, I., Tilling, I., Woitke, P., Thi, W. F., & Hogerheijde, M. 2010, A&A, 510, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Kamp, I., Thi, W.-F., Woitke, P., et al. 2017, A&A, 607, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Kassis, M., Deutsch, L. K., Campbell, M. F., et al. 2002, AJ, 124, 1636 [NASA ADS] [CrossRef] [Google Scholar]
- Kim, K. H., Watson, D. M., Manoj, P., et al. 2013, ApJ, 769, 149 [Google Scholar]
- Koumpia, E., de Wit, W. J., Oudmaijer, R. D., et al. 2021, A&A, 654, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Kramida, A. 2019, APS Meeting Abstracts, N09.004 [Google Scholar]
- Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 708, 1585 [NASA ADS] [CrossRef] [Google Scholar]
- Kraus, S., Hofmann, K.-H., Menten, K. M., et al. 2010, Nature, 466, 339 [Google Scholar]
- Kraus, S., Monnier, J. D., Che, X., et al. 2012, ApJ, 744, 19 [NASA ADS] [CrossRef] [Google Scholar]
- Lim, W., De Buizer, J. M., & Radomski, J. T. 2020, ApJ, 888, 98 [NASA ADS] [CrossRef] [Google Scholar]
- Lin, M.-K., Krumholz, M. R., & Kratter, K. M. 2011, MNRAS, 416, 580 [NASA ADS] [Google Scholar]
- Lodders, K. 2003, ApJ, 591, 1220 [Google Scholar]
- Lubow, S. H., & D’Angelo, G. 2006, ApJ, 641, 526 [Google Scholar]
- Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [Google Scholar]
- Mendigutía, I. 2020, Galaxies, 8, 39 [Google Scholar]
- Mendigutía, I., Calvet, N., Montesinos, B., et al. 2011, A&A, 535, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Mihalas, D. 1978, Stellar Atmospheres (New York: W.H.Freeman & Co Ltd) [Google Scholar]
- Mokiem, M. R., de Koter, A., Puls, J., et al. 2005, A&A, 441, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Mottram, J. C., Vink, J. S., Oudmaijer, R. D., & Patel, M. 2007, MNRAS, 377, 1363 [NASA ADS] [CrossRef] [Google Scholar]
- Müller, T. W. A., & Kley, W. 2013, A&A, 560, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Oliva, G. A., & Kuiper, R. 2020, A&A, 644, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Oudmaijer, R. D., Ababakr, K. M., & Fairlamb, J. R. 2017, Mem. Soc. Astron. Ital., 88, 605 [Google Scholar]
- Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415 [Google Scholar]
- Povich, M. S., Churchwell, E., Bieging, J. H., et al. 2009, ApJ, 696, 1278 [NASA ADS] [CrossRef] [Google Scholar]
- Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Rab, C., Güdel, M., Woitke, P., et al. 2018, A&A, 609, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ramírez-Tannus, M. C., Kaper, L., de Koter, A., et al. 2017, A&A, 604, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ramírez-Tannus, M. C., Backs, F., de Koter, A., et al. 2021, A&A, 645, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Reach, W. T., Megeath, S. T., Cohen, M., et al. 2005, PASP, 117, 978 [NASA ADS] [CrossRef] [Google Scholar]
- Rice, W. K. M., Armitage, P. J., Bonnell, I. A., et al. 2003, MNRAS, 346, L36 [NASA ADS] [CrossRef] [Google Scholar]
- Rivero González, J. G., Puls, J., Najarro, F., & Brott, I. 2012, A&A, 537, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Sarre, P. J. 2006, J. Mol. Spectr., 238, 1 [Google Scholar]
- Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [Google Scholar]
- Strom, K. M., Strom, S. E., Edwards, S., Cabrit, S., & Skrutskie, M. F. 1989, AJ, 97, 1451 [Google Scholar]
- Tan, J. C., Beltrán, M. T., Caselli, P., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson: University of Arizona Press), 149 [Google Scholar]
- Thi, W. F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711 [Google Scholar]
- van Gelder, M. L., Kaper, L., Japelj, J., et al. 2020, A&A, 636, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Vernet, J., Dekker, H., D’Odorico, S., et al. 2011, A&A, 536, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
- Wade, G. A., Bagnulo, S., Drouin, D., Landstreet, J. D., & Monin, D. 2007, MNRAS, 376, 1145 [NASA ADS] [CrossRef] [Google Scholar]
- Wade, G. A., Neiner, C., Alecian, E., et al. 2016, MNRAS, 456, 2 [Google Scholar]
- Wang, C., Langer, N., Schootemeijer, A., et al. 2022, Nat. Astron., 6, 480 [NASA ADS] [CrossRef] [Google Scholar]
- Wichittanakom, C., Oudmaijer, R. D., Fairlamb, J. R., et al. 2020, MNRAS, 493, 234 [Google Scholar]
- Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197 [CrossRef] [EDP Sciences] [Google Scholar]
- Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [Google Scholar]
- Xu, Y., Moscadelli, L., Reid, M. J., et al. 2011, ApJ, 733, 25 [NASA ADS] [CrossRef] [Google Scholar]
A reproduction package can be found at https://doi.org/10.5281/zenodo.7024846
All Tables
All Figures
![]() |
Fig. 1 Illustration of the normalization method. Top panel: the flux calibrated observed spectrum of B243 in black. The photospheric model spectrum is indicated by the orange dashed line and the continuum of this model by the green dash-dotted line. The model spectrum is matched to the observed spectrum using a least squares method in the unhatched wavelength ranges. At the top the location of the Paschen lines and a DIB are indicated. Bottom panel: the photosphere corrected and normalized observed spectrum (in solid black) and photospheric model (in dashed blue). The horizontal blue line indicates unity. |
In the text |
![]() |
Fig. 2 Line profile of Pa-14 from the model of B243 for various parameter values. In each of the panels one parameter is varied, while the others remain fixed at an inner radius of 0.06 AU, a reference scale height of 0.5 AU, a total disk mass of 10−3 M⊙, and an inclination of 75°. |
In the text |
![]() |
Fig. 3 Models of the SED for B243 for various parameter values. In each of the panels one parameter is varied, while the others are held fixed at an outer radius of 10 AU, a reference scale height of 0.04 AU, a disk mass of 10−4 M⊙, and an inclination of 15°. |
In the text |
![]() |
Fig. 4 Best fitting models in blue to the observed lines in black for B243. The blue shaded regions indicate models that fall within the 1σ confidence interval. Each line is fitted individually with the inclination set to 75°, the respective parameters are shown in Fig. 5. The orange profiles indicate the best fit when fitting Pa-16 to Pa-9 simultaneously. Notice the factor ~30 range in line strengths covered by the Paschen series. |
In the text |
![]() |
Fig. 5 Parameter values of the inner disk with their uncertainties for each modelled line of B243 and B331 as function of their oscillator strength. The horizontal dashed lines indicate the best fit parameters for B243 when fitting Pa-9 and Pa-11 to Pa-16 simultaneously and for B331 when fitting Pa-9, Pa-11 to Pa-16 and Br-11 to Br-16. The corresponding line profiles are shown in Figs. 4 and 7 for B243 and B331, respectively. For each fit an inclination of 75° is used. |
In the text |
![]() |
Fig. 6 Observed and best fit SED of B243 and B331. For each star, the scatter points indicate the dereddened observed photometry and the black dashed lines the Kurucz model used as input for the stellar spectrum in PRODIMO. The best fitting models are drawn in green, with the shaded region indicating the 1σ uncertainty. For B331 the photometry from Lim et al. (2020) is indicated as upper limits; these are not considered in the fitting procedure. |
In the text |
![]() |
Fig. 7 Best fitting models in blue to the observed lines in black for B331. The blue shaded regions indicate models that fall within the 1 σ confidence interval. Each line is Atted individually with the inclination Axed at 75°, the respective parameters are shown in Fig. 5. The orange proAles indicate the best At when Atting Pa-16 to Pa-9 and Br-16 to Br-11 simultaneously. Notice the factor ~20 range in line strengths covered by the Paschen series. |
In the text |
![]() |
Fig. 8 Hydrogen ionization, density, and temperature structure of the inner region of B243. The parameters used match the values indicated by the full red line in Fig. 5. The black contour indicates the approximate origin of 50% of the Pa-16 line emission based on vertical escape probabilities. The orange circle on the left indicates the central object. |
In the text |
![]() |
Fig. 9 Same as Fig. 8, but now for B331 on a logarithmic spacial scale. The parameters used match the values indicated by the dashed red line in Fig. 5. |
In the text |
![]() |
Fig. 10 Adapted from Fig. 9 of RT17. Positions of B243 and B331 in the Hertzsprung-Russel Diagram. The solid blue line on the left indicates the zero age main sequence. The dashed red line on the right indicates the birth line. The solid line connecting the birth line and the zero age main sequence indicate the MIST pre-main sequence evolutionary tracks of stars with masses indicated on the right (Dotter 2016). The colour of this line indicates the relative time between birth and reaching the main sequence. The dash dotted grey lines indicate the isochrones corresponding to 3 × 103, 104, 3 × 104, 105 and 3 × 105 yr from top to bottom respectively. |
In the text |
![]() |
Fig. 11 Schematic overview of the geometric properties of the disks of B243 and B331. |
In the text |
![]() |
Fig. 12 Orbital velocity of the hydrogen gas based on the peak separation of the observed (circles) and simulated (stars) emission lines as function of their oscillator strength. Green markers indicate Paschen lines and yellow markers Brackett lines. The velocity is determined by fitting double gaussians to the profile. The simulated profiles are from Figs. 4 and 7. |
In the text |
![]() |
Fig. A.1 Normalized and corrected spectra of B243 (top panel) and B331 (bottom panel), focusing on the Paschen and Brackett series. The spectrum of B243 is not normalized at wavelengths longer than 1500 nm due to the normalization method's incompatibility with the IR excess. |
In the text |
![]() |
Fig. B.1 H II region simulated with PRODIMO. The colour indicates the ratio between neutral and ionized hydrogen. The white dashed line indicates the analytically calculated Strömgren radius based on the radiation field, temperature and density of the medium. |
In the text |
![]() |
Fig. C.1 Corner plot indicating the χ2 value for the model fit to the Pa-12 emission line of B243 as function of each of the parameters and each parameter combination. The dashed line in each top panel indicates the 1σ threshold. |
In the text |
![]() |
Fig. C.2 Corner plot indicating the χ2 value for the model fit to the Pa-12 emission line of B331 as function of each of the parameters and each parameter combination. The dashed line in each top panel indicates the 1σ threshold. |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.