Modeling the Accretion Disk around the High-mass Protostar GGD 27-MM1

Understanding how high-mass stars form and evolve is one of the hot topics in astrophysics, due to the strong role that these objects play in the life of a galaxy. However, the study of high-mass protostars is difficult because of their fast evolution (∼10⁵ yr) to the main sequence and their large distances and high obscuration.

It is well known that low-mass stars are formed through an accretion disk that transports gas and dust from the envelope to the protostar and through a jet that removes the excess angular momentum (Shu et al. 1987; McKee & Ostriker 2007). Disks around nearby solar-type stars have been studied to great extent and detail (e.g., Williams & Cieza 2011; Testi et al. 2014; Hartmann et al. 2016; Andrews et al. 2018), but at the moment the number of disk studies of more distant and massive stars is still comparatively very small.

Accretion disks around massive stars are a plausible mechanism that can alleviate the radiation pressure problem, hence allowing an accretion flow to continue once photoionization has started (e.g., Tan et al. 2014; Klassen et al. 2016; Kuiper & Hosokawa 2018). However, their physical properties are still uncertain (see, e.g., Beltrán & de Wit 2016 for a recent review). Flattened, disk-like structures have been observed in a few massive, young stars (Beltrán et al. 2011, 2014; Sánchez-Monge et al. 2014; Johnston et al. 2015; Sanna et al. 2019; Zapata et al. 2019). However, these structures are large (∼1000–10,000 au) and have masses considerably larger than the central protostar, and it is difficult to envisage them as real accretion disks. Therefore, higher angular resolution and more sensitive observations are required to better characterize the physical parameters and the role of these rotating structures (e.g., Ginsburg et al. 2019; Maud et al. 2019). Three of the best examples in the literature of massive protostars associated with a clearly defined jet and a compact (few hundred astronomical units), angularly resolved dusty disk candidate are Cepheus A HW2 (Patel et al. 2005), GGD 27-MM1 (Girart et al. 2018), and G11.92-061 MM1a (Ilee et al. 2018).

The HH 80–81 objects (at a distance of 1.4 kpc; see appendix) are associated with a spectacular (∼10 pc long), highly collimated radio jet (Martí et al. 1993; Heathcote et al. 1998; Masqué et al. 2015), which is powered by a massive early B-type protostar IRAS 18162−2048 (GGD 27-MM1; Fernández-López et al. 2011a; Girart et al. 2017). This protostellar radio jet is the first one where polarized emission due to relativistic electrons has been detected, showing the presence of a magnetic field aligned with the jet (Carrasco-González et al. 2010; Rodríguez-Kamenetzky et al. 2017). Very Large Array (VLA) continuum observations at 7 mm reveal a cross-shaped morphology, which was interpreted as two overlapping structures that could correspond to the radio jet and a disk of ∼200 au radius (Carrasco-González et al. 2012), in agreement with the upper limit imposed by the 1.3 mm continuum dust emission observations (Fernández-López et al. 2011a). The size of the putative disk coincides with theoretical predictions of the centrifugal radius based on the spectral energy distribution (SED) fitting of some high-mass protostar regions (De Buizer et al. 2005).

In addition, observed velocity gradients in the molecular gas perpendicular to the HH 80–81 radio jet have been interpreted as rotating motions (Fernández-López et al. 2011b; Carrasco-González et al. 2012; Girart et al. 2017). Definitive evidence of a compact disk around IRAS 18162−2048 comes from Atacama Large Millimeter/submillimeter Array (ALMA) observations at 1.14 mm with an angular resolution of ∼40 mas (∼56 au), which reveal a compact dust disk clearly perpendicular to the radio jet (Girart et al. 2018; see Figure 1).

Zoom In Zoom Out Reset image size

Because of the similarities with disk–protostar–jet systems in low- and intermediate-mass protostars, in this paper we analyze the GGD 27–MM1 disk by applying models that have successfully explained disks around low-mass stars. The main goal is to investigate if the assumptions that are usually adopted for disks around low-mass stars can be roughly extrapolated to the case of massive stars. For example, disks around low-mass stars are much less massive than the central protostar, and therefore they are usually gravitationally stable. In the high-mass case, it is not clear if disks are stable (Maud et al. 2019) or unstable (Motogi et al. 2019; Zapata et al. 2019).

The paper layout is as follows. In Section 2 we present the ALMA observations of the disk around GGD 27–MM1, describing the disk model in Section 3. In Section 4 we present the main observational properties of the region. In Section 5 we explain the fit procedure, and Sections 6 and 7 correspond to results and discussion, respectively. We list the main conclusions in Section 8 and present the new distance estimation in the appendix.

In this work, we use ALMA continuum observations at 1.14 mm (263.0 GHz; project ALMA#2015.1.00480.S). The Band 6 receiver with the correlator set in continuum mode (time division mode) covering the 253.0–257.0 and 269.0–273.0 GHZ frequency ranges was used. The observations were carried out with 37 antennas in the c36-8/7 configuration, which provided baselines between 82 m and 11.05 km (13–5400 kλ). The Stokes I image toward GGD 27-MM1 was generated using the resulting visibilities after the subtraction of the compact source (see Section 2.1). This was done using the CASA task tclean with a value of 0.5 for the robust Briggs weighting parameter. Because of a lack of visibilities between 150 and 300 kλ and in order to filter extended emission coming from the envelope, we used visibilities from baselines larger than 300 kλ (calibration of the data is described in Girart et al. 2018). The resulting synthesized beam has an FWHM of 45.0 mas × 38.3 mas (PA = −624). The Stokes I rms noise is 60  μJy beam⁻¹. The ALMA#2015.1.00480.S project also had a science goal at Band 7 to observe several molecular lines in the 298–302 and 310–313 GHz frequency ranges (c40-6 configuration). Here we analyze the position–velocity diagrams of two of the brightest lines detected that better trace the disk kinematics: the SO₂ 9_2,8–8_1,7 and 19_3,17–19_2,18 transitions. For these two molecular lines, the calibrated data were self-calibrated from the continuum by using all the available channels in the four observed spectral windows, except for the channels with the brightest line emission (e.g., H₂CO 4_2,3–3_1,2). Velocity channel maps were obtained (after continuum subtraction) using tclean with natural weighting, which yielded an FWHM synthesized beam of 021 × 016 (PA = −871). The channel width was ≃0.98 km s⁻¹. The rms noise achieved was 1.3 mJy beam⁻¹ per channel. The other molecular lines detected will be presented in a forthcoming paper (M. Fernández-López et al. 2019, in preparation).

2.1. The Compact (Few Astronomical Units) Source at the Disk Center

Figure 2 shows the real part of the observed visibilities as a function of the uv distance from the phase center (disk center). The flux density decreases steeply with increasing visibility radius, $r(u,v)$, for the shortest baselines ($r(u,v)$ ≲ 2000 kλ). At larger radii, the flux density decreases more smoothly up to $r(u,v)$ ≃ 4000 kλ. At this point, the flux density remains roughly constant with the visibility radius. This suggests that the emission from the longest baselines may be dominated by a very compact object. To check this possibility, we generated several images using only visibilities with a minimum visibility radius, r(u, v)_min, between 3500 and 4750 kλ and with a robust weighting of 1. The values of the minimum visibility radius used and the resulting synthesized beam are given in Table 1. Maps including only baselines longer than 4000 kλ are devoid of artifacts, due to the severe missing flux density from the disk. The compact source appears in maps with long visibility radii (see Figure 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Gaussian Fit to the Compact Source

	Image Plane			Visibility Plane
r(u, v)min	Synthesized Beam	Flux Density	Deconvolved Size	Flux Density	Size
(kλ)	(mas × mas, °)	(mJy)	(mas × mas, °)	(mJy)	(mas)
3500	26.7 × 22.4, −63	28.5 ± 0.2	16 ± 1 × 7 ± 1, 91 ± 5	29.5 ± 0.1	13.6 ± 0.1
4000	23.8 × 20.9, −68	21.0 ± 0.2	unresolved	20.8 ± 0.1	4.7 ± 0.5
4250	21.9 × 20.9, −59	20.1 ± 0.2	unresolved	19.9 ± 0.3	3.5 ± 1.0
4500	21.6 × 19.8, 77	19.7 ± 0.3	unresolved	19.3 ± 0.4	faileda
4750	21.8 × 18.3, 67	19.1 ± 0.4	unresolved	19.5 ± 0.6	faileda

Note.

^aAlgorithm did not find a solution, probably due to the small range of visibilities.

Download table as:  ASCIITypeset image

A two-dimensional Gaussian fit was performed to the different images obtained with different baseline ranges. The flux density and the deconvolved size obtained from the Gaussian fit are listed in Table 1. The images with visibility radii ≥4000 kλ show an unresolved source with a flux density of ∼19 mJy at the center of the disk. We also performed a Gaussian fit to the visibilities using the same range of visibilities as before. The flux density and size of these fits are also shown in Table 1. The fits to the visibilities with r(u, v)_min > 4000 and >4250 kλ reveal that the emission arises from a very compact region with a radius of ∼4 mas (∼5.6 au). Given the relatively short range of visibilities used in this fit, further very high angular resolution observations (∼10 mas) are needed to better constrain its size. In any case, these values imply that the brightness temperature of this compact source is probably ∼10⁴ K.

The origin of the compact source cannot be due to the thermal emission of dust grains, but rather to ionized gas (see Section 7.1). Therefore, this compact source was removed from the visibilities to obtain the map presented here, tracing only the dust emission from the accretion disk (Figure 4).

Zoom In Zoom Out Reset image size

The disk was modeled using the irradiated α-accretion disk models created by D'Alessio et al. (1998, 1999, 2001, 2006), which have been successfully used and further developed to model disks around low- and intermediate-mass stars (e.g., McClure et al. 2013; Osorio et al. 2014, 2016; Macías et al. 2018).

The models describe disks around stars with parameters typical of classical T Tauris, that is, an irradiated flared disk with two populations of grains. These two populations aim to emulate the dust growth and vertical settling predicted by dust evolution models (Dullemond & Dominik 2004). The code computes the vertical and the radial structure of the disk using the α-viscosity prescription and enforcing vertical hydrostatic equilibrium. In the past, these models have already been used to reproduce the SED of a possible disk around the high-mass protostar AFGL 2591-VLA 3 (Trinidad et al. 2003). Furthermore, the α-viscosity prescription has been also used to model quasi-steady self-gravitating disks around massive protostars under certain conditions (H/R ≤ 0.1 and  ${M}_{\mathrm{disk}}$/M_* < 0.5, with H the disk scale height, and R the radius of the disk; Forgan et al. 2016).

The model allows us to set two populations of grains with a power-law size distribution n(a) ∝ a^−3.5, where a is the grain radius. Regarding the degree of settling, we used the epsilon parameter  = ζ_small/ζ_std, where ζ_small and ζ_std are the dust-to-gas mass ratio of the small grains and the initial standard value, respectively (see D'Alessio et al. 2006). The relative abundances of the different dust components were adopted with a dust-to-gas ratio of 0.0065 corresponding to the measured abundances of silicates and graphites in the interstellar medium (Draine & Lee 1984; D'Alessio et al. 2006; Osorio et al. 2014). The remaining ratio up to the commonly used value of 0.01 would be water ice and other ices. These ices should be sublimated at the high temperatures expected in the disk, so they should have a negligible contribution in the model (see Section 7.6).

The other parameter related to the settling is Z_big, which locates, as a function of the scale height (H), the border between the populations of grains, which was fixed as Z_big = 0.1H.

The main heating sources are the stellar irradiation and the viscous dissipation, which is parameterized through α (Shakura & Sunyaev 1973) and is assumed to be constant over the disk. The viscosity effective coefficient is defined as ν_t = αc_sH, where c_s is the local sound speed and H is the hydrostatic scale height of the gas:

$$\begin{eqnarray}&&\displaystyle \frac{H}{R}=\displaystyle \frac{{c}_{{\rm{s}}}\left({T}_{{\rm{c}}}\right)}{R{\rm{\Omega }}(R)}={\left[\displaystyle \frac{{{kT}}_{{\rm{c}}}R}{{{GM}}_{\mathrm{tot}}\mu {m}_{{\rm{H}}}}\right]}^{1/2},\end{eqnarray} \tag{ 1 }$$

where k is the Boltzmann constant, T_c is the disk midplane temperature, G is the gravitational constant, M_tot is the total mass (M_* + M_disk) at every radius, μ = 2.33 is the mean molecular weight, and m_H is the hydrogen mass. Besides, the model considers accretion luminosity as part of the irradiation of the disk.

The temperature and density structures are calculated self-consistently once the stellar parameters (radius R_*, mass M_*, and temperature T_*), the dust content (abundances, distribution of grain sizes), viscosity (α), and disk mass accretion rate (${\dot{M}}_{\mathrm{acc}}$) are set. The dust opacity includes absorption and isotropic self-scattering. In an α-accretion disk model, the mass surface density is ${\rm{\Sigma }}={\dot{M}}_{\mathrm{acc}}{\rm{\Omega }}/3\pi \alpha {c}_{{\rm{s}}}{({T}_{{\rm{c}}})}^{2}$. The remaining parameters to describe the disk model are the disk radius, R_disk, and the inclination angle of the disk i.

The disk is considered to be steady, axisymmetric, and geometrically thin. Its self-gravity is neglected compared to the stellar gravity, and it is assumed to be in Keplerian rotation and in hydrostatic equilibrium in the vertical direction. The model assumes that dust and gas are well mixed and thermally coupled; thus, a unique temperature as a function of position in the disk is calculated.

In this section, we present the main observational properties that can be used to constrain the parameters of our disk model. Figure 4 shows the ALMA continuum image of the GGD 27–MM1 disk at 1.14 mm. We resolve the disk at this wavelength, obtaining a flux density of 351.30 ± 0.33 mJy and a peak intensity of 46 mJy beam⁻¹ (see also Busquet et al. 2019). The morphology observed at this wavelength is consistent with an inclined disk (49°; 0° for a face-on disk) with a radius of ∼240 au and PA = 113°. The brightness temperature of the disk in the central region reaches a value of ∼470 K (Girart et al. 2018). This is an indication of an important source of heating.

Although there are not enough data points at different frequencies with high angular resolution to build the SED of the disk, we have extensive knowledge of the region, which provides us with a series of observational constraints regarding the physical parameters of our model, as discussed below.

Bolometric luminosity. The observed value is ∼2.0 ×10⁴ L_⊙ for a distance of 1.7 kpc (Gómez et al. 2003). As shown in the appendix, the source distance is between 1.2 and 1.4 kpc, and therefore the luminosity can be recalculated to be between ∼1.0 × 10⁴ and ∼1.4 × 10⁴ L_⊙. This luminosity must be considered an upper limit for the massive protostar GGD 27–MM1 because it comes mainly from the IRAS 18162–2048 fluxes that probably encompass other sources (Fernández-López et al. 2011a).

Dynamical mass. The total dynamical mass of the star–disk system can be obtained from the molecular line emission tracing the gas motions from the disk and assuming that they behave as a rotationally supported disk (i.e., Keplerian velocity). We used the data cubes for the SO₂ 9_2,8–8_1,7 and 19_3,17–19_2,18 lines (see Section 2) to construct position–velocity (PV) maps, centered at the dust peak intensity with a position angle of 113°, that is, along the major axis of the disk. Figure 5 shows the resulting plots. The brightest blueshifted emission appears in the southeast side of the disk, while the redshifted emission arises from the northwest side of the disk. This is in agreement with previous observations at lower angular resolution and less sensitivity (Fernández-López et al. 2011b; Carrasco-González et al. 2012; Girart et al. 2017). We note that there is also significant red/blueshifted emission in the southwest/northeast side of the disk. This could be an indication of infall motions, although MHD simulation of disk formation and evolution shows that this can be a projection effect for significant inclinations (e.g., Seifried et al. 2016). In order to constrain the dynamical mass from these PV maps, we followed the procedure given by Seifried et al. (2016). This procedure fits the Keplerian profile to the 5σ contour emission and was tested in synthetic ALMA molecular line PV maps generated from MHD simulations for disks around both low- and high-mass stars. It should best work for lines that fully resolve, spatially and kinematically, the Keplerian profile in the PV maps and for cases not too close to a face-on projection. The best fit obtained yielded an inclination-corrected dynamical mass of 31 ± 1 and 21 ± 1 M_⊙  for the SO₂ 9_2,8–8_1,7 and 19_3,17–19_2,18 lines, respectively. Since the line emission mostly arises from the outskirts of the dusty disk modeled in this paper, we can consider this dynamical mass as the combined mass of the star and the accretion disk. Thus we explored stellar mass values between 15 and 25 M_⊙ and keep the total mass (star + disk) close to the dynamical mass.

Zoom In Zoom Out Reset image size

Mass accretion rate. Carrasco-González et al. (2012) estimate the mass-loss rate in the jet using the formulation given by Reynolds (1986; see Equation (2)), who model the free–free emission from an ionized jet by adopting a power-law dependence with radius (see also Anglada et al. 2018):

$$\begin{eqnarray}\begin{array}{rcl}\left[\displaystyle \frac{{\dot{M}}_{\mathrm{out}}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right] & = & 1.9\times {10}^{-6}{x}_{0}^{-1}\left[\displaystyle \frac{{v}_{\mathrm{jet}}}{1000\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right]{\left[\displaystyle \frac{{S}_{\nu }}{\mathrm{mJy}}\right]}^{0.75}\\ & & \times {\left[\displaystyle \frac{\nu }{\mathrm{GHz}}\right]}^{-0.45}{\left[\displaystyle \frac{D}{\mathrm{kpc}}\right]}^{1.5}{\left[\displaystyle \frac{{\theta }_{0}}{\mathrm{rad}}\right]}^{0.75}.\end{array}\end{eqnarray} \tag{ 2 }$$

They assume a pure hydrogen jet with a constant opening angle θ₀ ∼ 19° (conical jet), terminal velocity v_jet (∼1000 km s⁻¹), ionization fraction x₀ = 0.1, and electron temperature T_e = 10⁴ K, and that the jet axis is in the plane of the sky. They obtain a mass-loss rate of ${\dot{M}}_{\mathrm{out}}\sim {10}^{-5}$  M_⊙ yr⁻¹ for a distance of 1.7 kpc, ≃8 × 10⁻⁶  M_⊙ yr⁻¹ for the corrected distance of 1.4 kpc (see the appendix). This value is in agreement with the value obtained from CO observations of the molecular outflow associated with the jet (Qiu et al. 2019). Assuming the accretion rate ${\dot{M}}_{\mathrm{acc}}$ of the disk onto the star is ∼10 times larger than the mass-loss rate (Bontemps et al. 1996), the mass accretion rate would be ${\dot{M}}_{\mathrm{acc}}$ ∼ 8 × 10⁻⁵  M_⊙ yr⁻¹. However, Beltrán & de Wit (2016, and references therein) obtain a ratio between the mass-loss rate and mass accretion rate of approximately ∼0.3 for disks in young, high-mass stars. In this case, the mass accretion rate would be lower, ∼3 × 10⁻⁵  M_⊙ yr⁻¹.

Even though the mass accretion rate is not well determined using these methods, it allows us to define a range of values to explore. Therefore, in our modeling, we explored values of ${\dot{M}}_{\mathrm{acc}}$  in the ∼1 × 10⁻⁵ to ∼2 × 10⁻⁴ M_⊙ yr⁻¹ range.

Stellar parameters. It is expected that most of the dynamical mass (21–31 M_⊙) is stellar, otherwise the disk would be unstable and would show significant asymmetries (e.g., spiral structures) that are not observed with the present data. Such a massive star, if it were in the main sequence, should develop an H II region, which is not detected with the present observations (see Section 7.1). A possible solution to mitigate this problem is to assume that the star is inflated, with low enough temperature to not produce stellar UV radiation and create an H II region.

In that sense, Hosokawa & Omukai (2009) found a dependence of the protostellar radius with the accretion rate. They obtain that, in general, the higher the accretion rate, the larger is the stellar radius. Then the protostar has a lower maximum temperature for a certain stellar mass. This fact causes a delay in the onset of the main-sequence phase and therefore a delay in the formation of an H II region. Thus, because of the high accretion rate estimated, we decided to explore large stellar radii, between 10 and 30 R_⊙. Because the luminosity is known, this implies temperatures between 12,000 and 18,000 K. This is consistent with Johnston et al. (2013), who modeled the envelope and disk around the luminous star AFGL 2591-VLA3 (2.3 × 10⁵ L_⊙) and find a stellar radius of 90 R_⊙ with a temperature of 13,000 K.

Inner disk radius. Considering a sublimation temperature of 1400 K for the most refractory grains (D'Alessio et al. 2006) and the observed luminosity, we can locate the sublimation wall from L_tot = 4πσT⁴R² at ∼12 au of radius. Because of the high spatial resolution of the observations (∼56 au), an inner radius larger than ∼20 au was discarded. Otherwise we should marginally resolve the inner wall of the disk. The sublimation temperature usually is assumed to cover a range from 1200 to 1500 K, so we explored an inner radius range between 10 and 20 au instead of just using the 12 au that we calculated (see Table 2).

Settling degree and grain size. We can obtain some constraints on the settling of the larger dust grains at the midplane by looking at the polarization emission due to scattering from large grains. Based on polarization observations, Girart et al. (2018) do not find signs of dust settling, and they determine a dust maximum grain size (a_max) from 50 to 500 μm (see Section 7.2). We explored in our modeling cases with different settling degrees, including the no-settling case.

Taking into account the observational restrictions set out above, in the following we proceed to look for the set of parameters that best fits the observations.

We computed several grids of models by varying the parameters of the disk and the star. Through these grids, the parameters were refined until the best-fit model image was obtained.

The fitting procedure and analysis was done using Common Astronomy Software Applications (CASA) and Multichannel Image Reconstruction, Image Analysis, and Display (MIRIAD; Sault et al. 1995) data reduction software packages.

In order to properly compare the model with the ALMA data, synthetic visibilities were computed from the model images. This was done with the CASA simutil package. After that, tclean in CASA was used with the simulated visibilities to create a final image from the model, adopting the same cleaning parameters used to create the observed ALMA image. In Figure 6 we show the ALMA image (top panel), the modeled image (middle panel), and the residual map (bottom panel). The residual map was obtained by subtracting the model image from the observed one. Therefore, positive values in the residuals show regions where the model underestimates the emission.

Zoom In Zoom Out Reset image size

The radial intensity profile was computed by averaging in concentric ellipses every 001 (∼1/4 of the beam), with the inclination adopted in the disk model (see Figure 7). The best-fit parameters of the disk model were obtained initially by visual inspection of the radial intensity profile, and then, from among this first selection, we chose the best-fit model based on the minimum χ² of the residual map. In the case of the residual map, only pixels inside of an ellipse with R = 230 au were considered to calculate the χ² (see Figure 4) to avoid contamination from the outskirts of the map. Appendix B shows the variation of the disk parameters with χ² (leaving the other parameters fixed).

Zoom In Zoom Out Reset image size

In this section, we present the best-fit model obtained after exploring a wide range of values in the space of parameters of the model presented in Section 3. The best-fit parameters are shown in Table 3.

We found a massive disk of ∼5 M_⊙, that is, with 20% of the stellar mass. The radius of the disk is ∼170 au with an inclination angle of ∼49° (angle between the rotation axis and the line of sight). The mass accretion rate, which in our model is a very sensitive parameter, resulted to be ∼7 × 10⁻⁵  M_⊙ yr⁻¹. The total luminosity is ∼1 × 10⁴ L_⊙, which is in agreement with the previous estimation (see Section 4). The stellar mass is 20 M_⊙, which together with disk mass (5 M_⊙) is consistent with the estimated dynamical mass for the star–disk system (21–30 M_⊙; see Section 4).

Regarding the composition of the disk in terms of grain size, different maximum sizes of grains according with the results found in Girart et al. (2018) were tested (see Section 4). We did not find substantial differences between models with a_max of 100 and 500 μm and with 1 mm (see Figure 8). Only for small radii (<005) is the difference between the models noticeable. The grain sizes of the best-fit model go from a minimum value of 0.005 μm to a maximum of 3 μm in the disk upper layer. For those grains settled in the disk midplane, the grain sizes are between 5 and 500 μm.

Zoom In Zoom Out Reset image size

The density and temperature profiles for the best-fit model are shown in Figure 9. We found a flared disk with a maximum scale height of ∼13 au. The disk shows a temperature profile that goes from ∼1400 K at the inner edge of the disk to ∼150 K at the outer part. The small irregularities that can be seen in the internal part of the disk are caused by numerical effects and by the sublimation of the different dust components, which result in a step in the dust opacity, but they do not affect the results.

Zoom In Zoom Out Reset image size

From the mass surface density (Σ) and using the total opacity of the model (absorption + self-scattering), χ = 0.13 cm² g⁻¹ (opacity for a dust grain mixture with the physical properties described in Section 3 and with a grain size distribution that assumes grains with a maximum radius of 500 μm and considers the scattering effects), we obtain an optically thick disk at 1.14 mm for all radii, with τ = Σχ ranging from 50 to 170. Busquet et al. (2019) computed the mass of gas and dust of the disk, assuming that the 1.14 mm dust continuum emission is optically thin and the temperature distribution is uniform (T_d = 109 K). They estimated a disk mass for GGD 27 MM1 of ∼0.5 M_⊙. This mass could be considered as a lower limit, due to the optical thickness of the disk and to the fact that the opacity due to the self-scattering is not considered. We fitted the density and the temperature profiles with power-law functions (Σ(R) ∝ R^p, T(R) ∝ R^q) using the method of the minimum mean squared error (MSE). Here we present the coefficient and the power index. Equation (3) shows the behavior of the surface density, Σ(R), and the temperature at the midplane, T_c(R), approximated as power laws. The MSE of the fits are 0.023 and 0.034, respectively:

$$\begin{eqnarray}&&\begin{array}{c}\displaystyle \frac{{\rm{\Sigma }}(R)}{\left[{\rm{g}}\,{\mathrm{cm}}^{-2}\right]}\sim \ 500{\left(\displaystyle \frac{R}{\left[100\mathrm{au}\right]}\right)}^{-0.5}\\ \displaystyle \frac{{T}_{{\rm{c}}}(R)}{[{\rm{K}}]}\sim \ 300{\left(\displaystyle \frac{R}{\left[100\mathrm{au}\right]}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 3 }$$

In this work, we have studied the ALMA image at 1.14 mm of the circumstellar disk around GGD 27-MM1. We found a compact source coming from the inner radius (∼4 mas/5.6 au) that is probably due to ionized gas. By modeling the dust continuum emission from the disk (with the compact source previously subtracted), we found a massive (∼5 M_⊙) and compact (∼170 au) disk. In the following, we discuss the implications of our results, while also analyzing the gravitational stability of the disk.

7.1. Ionized Component

7.2. Restrictions from Polarization

An important source of information about the composition of the disks in terms of grain size, grain shape, and grain distribution comes from polarimetric observation.

Polarization data allow us to constrain the dust distribution in the disk (settling) and the maximum dust grain size. Girart et al. (2018), based on polarization models of Yang et al. (2017), conclude that dust settling has not yet occurred. Yang et al. (2017) compare two models with different thicknesses of the layer of large grains that are responsible for the scattering, and they propose two models in which large grains can be found up to 0.1 H' and 1 H', where H' is the (dust) hydrostatic scale height, $H^{\prime} {(R)={H}_{0}^{{\prime} }(R/{R}_{{\rm{c}}})}^{1.5-q/2}$, where R_c is a characteristic radius of the disk (dust) density distribution, and q is the temperature power-law index.

Unlike the models of these authors, our models consider the dust settling changing the dust-to-gas mass ratio between two populations through the parameter, as explained in Section 3. Small grains are in the upper layer, and big grains are in the midplane. The position of the border between these two populations can be set through the parameter Z_big.

In this work, our best-fit model has  = 1 (no dust settling), but it still has two dust populations with a border located at Z_big = 0.1H, where H is the pressure scale height of the gas at the midplane temperature (see Section 3).

We also reproduced the scenario suggested by Yang et al. (2017) to reproduce the polarization pattern observed by Girart et al. (2018). To do so, we tested a model with Z_big = 1 H and  = 1. In this case, we obtained a more massive disk (7 M_⊙) with larger mass accretion rate (1 × 10⁻⁴  M_⊙ yr⁻¹) than the model with Z_big = 0.1H (M_disk ∼ 5 M_⊙ and ${\dot{M}}_{\mathrm{acc}}$ ∼ 7 ×10⁻⁵  M_⊙ yr⁻¹; see Table 3), our fiducial best-fit model, which shows a lower value of χ².

Table 3.  Best-fit Model Parameters

Parameter	Value
M*(M⊙)	20	fitted
Teff (K)	12,000	fitted
R* (R⊙)	25	fitted
Distance (kpc)	1.4	adopted
$\dot{M}$ (M⊙ yr−1)	7 × 10−5	adopted/refineda
Mdisk (M⊙)	5	calculated
Rdisk (au)	168	adopted/refineda
H100 (au)	7	calculated
Rin (au)	14	fitted
i (°)	49	adopted/refineda
amax (μm)	500	fitted
α	0.1	fitted

Note.

^aParameter with observational constraints and uncertainty.

Download table as:  ASCIITypeset image

Furthermore, we have compared two- and one-population models. We tested models with a single dust grain population with a maximum grain size of 500 μm. Due to the lower opacity to the stellar radiation of this dust population, the one-population models would require significantly higher disk masses, even comparable to the stellar mass. These masses are inconsistent with our estimated dynamical mass (see Section 4), indicating that at least two dust populations are needed to fit the observations.

In summary, we have considered the results presented in Girart et al. (2018) concerning the size and distribution of grains. Moreover, taking into account the conceptual differences presented above (H, H'), we have explored values adjacent to them. We did not find substantial differences between both scenarios, and therefore we cannot favor either.

7.3. Stability of the Disk

We have quantified the Toomre parameter of our model to check the stability of the disk against self-gravity perturbations. Figure 9 (bottom right) shows the Toomre parameter Q (Equation (4)) for a Keplerian disk evaluated at the disk midplane temperature,

$$\begin{eqnarray}&&\begin{array}{l}Q=\,\displaystyle \frac{{c}_{{\rm{s}}}{\rm{\Omega }}}{\pi G{\rm{\Sigma }}}\\ {\rm{\Omega }}=\,{\left(\displaystyle \frac{{{GM}}_{\mathrm{total}}}{{R}^{3}}\right)}^{1/2},\end{array}\end{eqnarray} \tag{ 4 }$$

where c_s is the sound speed at the midplane temperature, Ω the Keplerian angular velocity, Σ the surface density of the disk, and G the gravitational constant.

We adopted the surface density and the temperature at the midplane as a function of the radius of our best-fit model. The stability condition (Q > 1) is satisfied up to a radius R_disk ∼ 100 au. A similar result has been found by Maud et al. (2019) in the massive O-type protostar G17.64+0.16, who reported a massive and stable disk for R_disk ≤ 150 au. Furthermore, Meyer et al. (2018) and Takahashi et al. (2016) propose an update of the Toomre criterion in which the only necessary condition for disk instability is Q < 0.6.

In addition, our best-fit model satisfies the conditions proposed by Forgan et al. (2016) in order to apply the α-viscosity prescription model in self-gravitating disks; that is, H/R < 0.1 (see Figure 9) and M_disk/M_* < 0.5.

We note that, due to the relatively high mass of the disk of our best model (5 M_⊙), self-gravity might play a nonnegligible role. In order to explore the potential effects of including self-gravity, we compared our best-fit model to the hydrodynamic simulations (including self-gravity) of massive star formation performed by Kuiper & Hosokawa (2018). Our model shows very similar results in terms of disk midplane temperature and the disk's aspect ratio to this model (Figure 6; Kuiper & Hosokawa 2018), which indicates that including self-gravity would not change our results significantly.

7.4. Residual Map

The values of the residual map are low when compared with the observed (residuals are below 5% of the peak intensity), but significant with respect to the rms noise level (σ ∼ 0.06 mJy beam⁻¹), with an intensity range between ∼2 and −2  mJy beam⁻¹ (see Figure 6). The main differences between the observed and modeled images arise from the outer parts of the disk, beyond 150 au, with an excess and a deficit of emission. This is illustrated in Figure 10, where we show flux density-position cuts along the major and minor axes of the disk (observed and modeled). We also verified that the asymmetries were not caused by a mismatch in the inclination of the system or a shift in the disk centers. In Figure 11 we present the model image and the residual map of the best-fit model (i = 49°; central panels) together with two models in which the inclination around the best fit has been modified (i = 44° and 54°). In addition, in Appendix B we show in the χ² as a function of the inclination for the best-fit model, varying the inclination angle in a range of 10° centered at 49°.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although we cannot discard intrinsic asymmetries within the disk, a plausible explanation for this excess or defect of emission at the outer parts of the disk could be a small mismatch with the flaring angle or settling of the disk. To discriminate between possible intrinsic asymmetries and modeling, we would need ALMA high angular resolution observations at other frequencies as have been carried out in low-mass YSOs (e.g., Carrasco-González et al. 2019; Macías et al. 2019).

7.5. Mass Accretion Rate and Evolutionary Stage

Trinidad et al. (2003) used the D'Alessio model to fit the SED of AFGL 2591 VLA3, a massive disk–star system located in the Cygnus X region. They find that for all models, the main heating source was the stellar irradiation for radii larger than 20 au and the viscous dissipation for smaller radii. Furthermore, ${\rm{\Sigma }}\sim {\dot{M}}_{\mathrm{acc}}/\alpha$ for radius larger than 20 au. In this scenario, they find a family of models that could explain the observed SED. This family of models corresponds to a constant value of ${\dot{M}}_{\mathrm{acc}}$/α, showing that the SED does not change as long as this ratio is maintained.

Our best-fit model yields a rate ${\dot{M}}_{\mathrm{acc}}$/$\alpha \sim 1\,\times {10}^{-4}$  M_⊙ yr⁻¹. We tested if the radial intensity profile would be significantly affected by varying the accretion rate and alpha, while keeping ${\dot{M}}_{\mathrm{acc}}$/α constant. As we show in Figure 12, the radial intensity profile at 1.14 mm is dramatically affected when the disk mass accretion rate is changed, although the ratio ${\dot{M}}_{\mathrm{acc}}$/α stays constant. The main reason for this difference is the effects of the irradiation from the accretion luminosity. Therefore, by not having this degeneration, we are able to constrain the ${\dot{M}}_{\mathrm{acc}}$.

Zoom In Zoom Out Reset image size

We would like to point out that the accretion rate could be variable with time. In fact, Martí et al. (1998), based on multiepoch VLA continuum observations, report a flux density decay of the two inner condensations in the HH 80–81 thermal radio jet. Such a flux density decay could be attributed to changes in the mass accretion rate being higher in the past. Furthermore, GGD 27–MM1 is a young source that has a faint envelope and probably an incipient (hypercompact) H II region. Some studies estimate that the timescale for the development of H II regions, with an accretion rate in the range ∼10⁻⁴–10⁻³  M_⊙ yr⁻¹, is ∼10⁵ yr (e.g., Cesaroni 2005). Considering the constant accretion rate of our best-fit model (7 × 10⁻⁵ M_⊙ yr⁻¹), for it to reach a star mass of 25 M_⊙, its age would be ∼4 × 10⁵ yr, which is larger than the development time of an H II region. Thus, higher accretion episodes in the past are necessary to explain the present situation.

The GGD 27 complex includes two compact cores, MM1 and MM2, separated by ∼7'' (∼10,000 au). Fernández-López et al. (2011a) estimated masses of MM1 and MM2 at different scales (see their Table 6). These authors show that while in MM2 most of the mass (∼75%) is found at envelope scale, in MM1 ∼70% of the mass is already at disk scale. This fact would place MM1 in a more evolved stadium than MM2, with MM1 being equivalent to a Class I low-mass star.

7.6. Temperature Distribution in the GGD 27-MM1 Disk: Comparison with Disks around Low- and Intermediate-mass Stars and Implications for the Water Snow Line

We used ALMA continuum observations at 1.14 mm, obtained with an angular resolution of ∼40 mas, to model the accretion disk around the central massive early B-type protostar exciting the powerful HH 80–81 radio jet using the α-viscosity prescription. We found an enclosed mass of 21–30 M_⊙, of which 5–7 M_⊙ can be attributed to the disk. This mass is consistent with the derived dynamical mass of 31 ± 1 M_⊙ and 21 ± 1 M_⊙ for the SO₂ 9_2,8–8_1,7 and 19_3,17–19_2,18 lines, respectively. The radius of the disk is ∼170 au, with an inclination angle of 49°. We compared the physical structure, temperature, and density profiles obtained with our model with power-law functions, showing that the GGD 27–MM1 system is a potential template for future similar studies in other high-mass protostars. In particular, we obtained a flared disk with a maximum scale height of ∼13 au and a temperature profile that ranges from ∼150 K at the outskirts of the disk up to ∼1400 K at the inner edge of the disk. The analysis of the Toomre Q parameter, evaluated at the disk midplane temperature, indicates that the disk is stable up to a radius R_disk ≃ 100 au. This work shows that the D'Alessio models can be used as a first approximation in the modeling of accretion disks around massive protostars, providing in addition several observational predictions.

We also reported the presence of an unresolved compact source at the center of the accretion disk, with a radius of 4 mas (∼5.6 au at the source distance of 1.4 kpc) and a brightness temperature of ∼10⁴ K, most likely tracing ionized gas. The origin of this compact source is uncertain; it could arise from an incipient, extremely compact H II region or from the base of the HH 80–81 radio jet. Observations at higher angular resolution would help to determine the nature of this compact source.

Finally, we have estimated a distance of 1.2–1.4 kpc to the GGD 27 star-forming region based on the Gaia DR2 catalog combined with near-infrared polarimetric data of the YSOs in the region and the extinction maps.

We are grateful to Nuria Calvet for making D'Alessio's code available to us and for very useful discussions. We would also like to deeply thank Guillem Anglada for his critical comments on this work that have undoubtedly improved it very substantially. We thank our anonymous referee for her or his very helpful comments and suggestions to improve our manuscript. We also thank Rosine Lallement for the support on STILISM. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.00480.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. N.A.-L., M.O., J.M.G, G.B., R.E., and J.M.T. are supported by the Spanish grant AYA2017-84390-C2-R (AEI/FEDER, UE). R.E. is also supported by MDM-2014-0369 of ICCUB (Unidad de Excelencia "María de Maeztu"). M.O. acknowledges financial support from the State Agency for Research of the Spanish MCIU through the "Center of Excellence Severo Ochoa" award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). S.C. acknowledges support from DGAPA, UNAM, and CONACyT, México. M.F.-L. acknowledges support from the LACEGAL project, which has received funding from the European Union's Horizon 2020 Research and Innovation Program under the Marie Sklodowska-Curie grant agreement No. 734374. C.C.-G. and R.G.M. acknowledges support from UNAM-DGAPA PAPIIT project IN 108218 and IN 104319. J.K. is supported by MEXT KAKENHI grant No. 19K14755.

Software: CASA (McMullin et al. 2007), MIRIAD (Sault et al. 1995), PYTHON, STILISM.

The GGD 27 nebulosity and the objects associated with the region (e.g., HH 80–81 and the HH 80N star-forming core; Girart et al. 1994, 2001; Heathcote et al. 1998; Masqué et al. 2011, 2013) are located in the LDN 291 large molecular cloud complex (including the dark clouds LDN 306, 314, 315, and 322) that extends ∼4° in Sagittarius (Lynds 1962; Saito et al. 1999; Reipurth et al. 2008). The distances reported in the literature range between 1.5 and 2.4 kpc (e.g., Racine 1968; Humphreys 1978; Rodriguez et al. 1980). However, the most used value in the literature is 1.7 kpc. For this distance, the spatial scale of the LDN 291 cloud complex is ∼75 pc × 19 pc, and the mass is ∼1.2 × 10⁵  M_⊙ (Saito et al. 1999).

A.1. Analysis

Here we present two different approaches to better estimating the distance to the LDN 291/GGD 27 region. The two approaches rely on the Gaia DR2 catalog (Gaia Collaboration et al. 2016, 2018). First, we used data from Gaia in combination with near-infrared polarimetric data of the YSOs located in the GGD 27 region (Kwon et al. 2016). The second approach is to use the online Structuring by Inversion the Local Interstellar Medium (STILISM¹² ) application (Capitanio et al. 2017; Lallement et al. 2018), which combines the Gaia data with extinction maps to obtain 3D dust maps of the Galaxy.

A.1.1. Method 1: YSOs in GGD 27 and Gaia

Figure 13 shows the polarization fraction as a function of the Gaia DR2 parallax for the YSOs in GGD 27 (from Kwon et al. 2016), and this figure shows that the YSO populations with counterparts in Gaia have distances clearly smaller than the distance adopted in the literature, 1.7 kpc. Most of the YSOs have parallaxes between ∼0.6 mas (1.67 kpc) and 1.2 mas (830 pc) but with significant uncertainties. However, there are two objects that have large parallaxes. One of them is [HL85] GGD 27-28 31 (Hartigan & Lada 1985). It is located at a distance of 362 ± 58 pc but has very high polarization levels (26%, 42%, and 57% in the JHK bands, respectively; Kwon et al. 2016). This star is located in front of the bright GGD 27 nebula, where high levels of circular polarization are detected. Therefore, the near-infrared linear polarization is likely coming from the nebula, and it is not related to the optical star. The other star, 2MASS J18185959-2045537, is located at 395 ± 31 pc, and it exhibits no polarization in the near-infrared, which indicates that this is a (cold) foreground star not related to the cloud. In order to estimate the distance, we calculated for the YSOs the average value of the Gaia parallaxes weighted with the uncertainty,

$$\begin{eqnarray}&&\langle {\pi }_{\mathrm{GGD}27}\rangle =\displaystyle \frac{{\sum }_{i}^{N}{\pi }_{i}/{\sigma }_{i}^{2}}{{\sum }_{i}^{N}1/{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 5 }$$

and the uncertainty is

$$\begin{eqnarray}&&\sigma (\langle \pi \rangle )=\sqrt{\displaystyle \frac{N}{{\sum }_{i}^{N}1/{\sigma }_{i}^{2}}}.\end{eqnarray} \tag{ 6 }$$

Here, π_i and σ_i are the parallax for each star and its uncertainty, respectively, from the Gaia DR2 catalog. We excluded four stars with parallaxes, within their uncertainties, larger than 0.9 mas (with distances less than ≃1100 pc). Using the rest of the sample, we obtain an average, weighted by the uncertainty, parallax of 0.801 ± 0.106 mas. Therefore, the distance to the YSO cluster is 1248 ± 166 pc.

Zoom In Zoom Out Reset image size

A.1.2. Method 2: Extinction-Gaia Data—STILISM

A recent work (Danielski et al. 2018) has correlated the Gaia distances and the corrected version of G extinction with archival ground-based data, especially with 2MASS and SDSS/APOGEA-DR14 (Capitanio et al. 2017; Lallement et al. 2018). This allows us to derive 3D maps of the extinction as a function of the distance. We used the online tool that provides the cumulative reddening curve as a function of the distance for a given line of sight. Figure 14 shows the E(B − V) extinction as a function of the distance toward the LDN 291/GGD 27 region. It clearly shows two abrupt increases of the extinction, a small one (E(B − V) ∼ 0.2 mag or A_V ∼0.6 mag) around 100 pc and a larger one (E(B − V) ∼ 0.7 mag or A_V ∼ 2.0 mag) around 1200 pc. In order to fit the two abrupt extinction jumps, we used the following approach:

$$\begin{eqnarray}&&E(B-V)=\displaystyle \frac{{A}_{0}}{1+{e}^{-{A}_{1}({D}_{}-{D}_{\mathrm{jump}})}}+{A}_{2}D\,+{A}_{3},\end{eqnarray} \tag{ 7 }$$

where D is the distance, D_jump is the distance where the jump occurs, and A_i (i = 0, ..., 3) are free parameters. We used a reduced χ² fit. The distance for the large jump is 1270 ± 65 pc. We also used this expression to estimate the first small extinction jump, obtaining a distance of 119 ± 15 pc.

Zoom In Zoom Out Reset image size

A.1.3. Distance to the GGD 27 Nebula/Molecular Cloud

The YSOs detected with Gaia are likely near the surface of the cloud; otherwise the cloud extinction would make them not visible at optical wavelengths. Indeed, most of the Spitzer YSOs from Qiu & Zhang (2009) do not have optical counterparts. The previous analysis indicates that the average distance within 1σ uncertainty is between ∼1100 and 1400 pc. The second method is even more very sensitive to the cloud's surface, giving a distance between 1200 and 1340 pc (also at 1σ level). Therefore, combining both methods, we can constrain the distance to the GGD 27 region in the range of 1200–1400 pc.

As we already advanced in Section 5, Figure 15 shows how the χ² value changes when we vary the main parameters of the disk while fixing the others to the best-fit value. The same scale for the Y-axis (χ²) is used in the four panels. The parameters under study were varied between ∼10% and ∼30% according to the observational restrictions (see Section 4). The four parameters show a minimum at the value of our best-fit model. The panel shows that the model is much more sensitive to changes in inclination and disk radius, where the parameter variations are around 10%, compared to changes in the inner radius and accretion rate, where the variations are over 30%.

Zoom In Zoom Out Reset image size

Regarding the inner radius, we situated the sublimation wall at 12 au (±2 au), based on the luminosity and the dust sublimation temperature. Inward from this border the dust cannot survive, so we cannot obtain a physically consistent model with an inner radius smaller than 14 au. In addition, an inner radius larger than ∼20 au (∼14 mas) is discarded because it should have been observed with the angular resolution of our ALMA observations.