Steven D. Doty, Niranjan Reji, Angela Sommerer
Department of Physics and Astronomy, Denison University, Granville, USA.
DOI: 10.4236/ijaa.2025.154025   PDF    HTML   XML   71 Downloads   361 Views   Citations

Abstract

We have performed self-consistent, three-dimensional, continuum radiative transfer modeling in the disk of HD 100546. Following observations and previous work, the model includes gaps in the disk. The model also includes variations in dust size, composition, and distribution. A best fit model is achieved over a wide wavelength range from 2 - 1000 mm. While we generally confirm the results of previous models, we find evidence for dust settling and shattering in the outer disk. We also find evidence for dust segregation in the inner disk, and inferential evidence for an irregular and tenuous inner region ~1 AU, which may be due to early disruption of the disk. Such settling is an important process in planet formation. The location of the settling in the outer disk, and the existence of segregation and irregularity in the interior, suggest that the race between disruption, settling, and planet formation is on-going in HD 100546, and that the timescales for the various processes vary throughout the disk at a given epoch.

Keywords

Protostellar Disks, Dust Evolution

Share and Cite:

1. Introduction

The process of star-formation due to the collapse of cores within giant molecular clouds, and the intermediation formation of disks around the resulting protostars is relatively understood. As the luminosity and magnetic field of the protostar increase, the disk is disrupted and dispersed. It is in this intervening time between disk formation and disk disruption that planetary systems can be formed.

While over 6000 exoplanets have been observed to date, observing them in the process of formation is difficult. Consequently, the understanding of planet formation requires significant reliance on models. However, in the past decade, gaps and rings within protoplanetary disks have been observed. One particular example is the 10 Myr old Herbig Be star HD 100546. Previous work has observed a large disk of hundreds of AU, with potential signs of a tenuous inner disk r < 4 AU, and a gap from ~4 - 13 AU [1], with potential other rings between ~15 - 40 AU and again between ~110 - 250 AU [2].

In this work, we model continuum data publicly available for the SED of HD 100546, and apply an unbiased machine-fitting to the spectrum. The data are discussed in Section 2. The modeling approach and fitting in Section 3. The results are described in Section 4, and in Section 5 we discuss the implications of these fit.

2. Data

We obtained observational data for the 2.4 to 196.8 micron range from the ISO data archive. ISO covered the 2.4-to-45-micron range using its Short-Wave Spectrometer [3] and covered the 45-to-196.8-micron range using its Long-Wave Spectrometer [4]. This data was complemented by data from the Herschel Space Observatory—we obtained data from the Herschel Science Archive covering the 194-to-672-micron range obtained using the SPIRE (Spectral and Photometric Imaging Receiver) on the telescope [5]. These data were stitched together and cleaned of outlier data points using the astropy sigma_clip routine. Outliers beyond 3 standard deviations away from the rest of the data were iteratively removed. Subsequently, data points with a signal to noise ratio below 3 were clipped as well. As seen in Figure 1, the ISO and Herschel data are able to be joined without discontinuities once calibration differences are accounted for. More recent ALMA observations [6] [7] are also consistent with the SED. No epoch-based filtering was applied because the ALMA/ISO/Herschel overlap overlap smoothly within their uncertainties, indicating that inter-epoch variability and unresolved instrumental offsets should not meaningfully contaminate the derived disk structure.

Figure 1. The SED data utilized in this study. The errorbars are the uncertainties quoted by the source measurements. The details of the data sources are described in the text.

3. Source and Model

3.1. Source

The source HD 100546 is well observed. It has been established as an axisymmetric, ring system with mild gas irregularities [6]. Later observations [7] show that when integrating over the disk, the dust continuum is steady, even if localized gas motions occur. More recent CO observations [8] show a Doppler-flip in ¹²CO (2-1) which can be viewed as a transient, localized gas event (a “disk eruption”)—either vertical motion, a disk wind, or planet-induced perturbation—superposed on the previous steady dust structure. Indeed, an estimate of the dust mass involved in the disk eruption based on the line strength and dust/gas ratio of 1:100 yields an involved dust mass of ~10⁻⁴ x the total dust mass within 28 AU, consistent with only a small effect on the SED.

3.2. Model Description

The dust radiative transfer modeling is undertaken with the RADMC3D radiative transfer code [9]. This code uses a Monte Carlo approach to propagate photon packets to self-consistently solve for the temperature distribution within the source and the emitted radiation as a function of wavelength and position. The code inputs include the dust properties, density profiles in three dimensions, source properties, etc. The code has been carefully validated and used widely in the astrophysical community. In the model implementation here, 148 radial grid cells were utilized, with 64 in the polar direction. Both grids are approximately logarithmic, with more points where densities and/or their gradients are larger. Axial symmetry is assumed in concert with the adopted disk structure. A wavelength grid of 640 wavelengths is used in approximately logarithmic spacing, with more points added in regions where the dust opacities vary significantly. The disk is assumed to be inclined at an angle of 45 degrees, consistent with observations [7] [10].

We varied the physical structure of our model as we improved our fit to the data above. Our current model splits the circumstellar disk around HD 100546 radially into 3 regions specified by an inner radius which is fit, an outer radius which is fit, and two demarcation lines one-third and two-thirds of the way between the inner and outer radius. These two demarcation lines allow for the disk to be divided into three regions for the purpose of varying grain properties, and are based on experience gained in the fitting process.

We allow for variation in the mix of dust species within each region. We adopt three dust compositions: amorphous Olivene [11] [12] and Ossenkopf & Henning (OH) [13] dust with thin ice mantles and with no ice mantles. The OH component with no mantle is used up to the radius at which disk temperatures drop to around ~100 K and is replaced with the thin mantle component outside of that radius.

We assign each dust species i a non-negative weight d_i. The fraction of dust mass in species i is the weight divided by the sum of all weights for all species, namely.

$f_i=d_i/{\displaystyle {\sum }_j{d}_j}$ (1)

The species surface density is simply the product of each f_i, with the dust surface density of the disk. Starting points for the star mass, star temperature, the circumstellar disk’s surface density power law, the gas-to-dust ratio, disk mass, and the disk’s outer radius were taken from Bruderer [14]. The star’s radius is derived from the temperature and mass using the Stefann-Boltzmann relation. The form for the structure of the disk is largely derived from Chiang & Goldreich [15]. The scale height is calculated by computing the grazing angle for the disk,

$\alpha \sim 0.4R_{*}/a+ad\left(H/a\right)/da$ (2)

using it to calculate the temperature of the disk,

$T_e\sim {\left(\alpha /2\right)}^{1/4}{\left(R_{*}/a\right)}^{1/2}{T}_{*}$ (3)

and then we calculate the scale height using the previously calculated quantities.

$h/a={\left(T_e/T_c\right)}^{1/2}{\left(a/R_{*}\right)}^{1/2}$ (4)

The scale height for each dust species is multiplied by a flaring factor, h_i, allowing us to vary scale heights for various species and observe whether our model hints at dust settling. Our surface density profile is normalized to the mass of the circumstellar disk

${\displaystyle \sum \left(R\right)} = M_{disk} \left(p+2\right)R^p/\left[2\text{π}\left(R_{out}^{p+2}- R_{in}^{p+2}\right)\right]$ (5)

which we obtain from assuming that the radial surface density is proportional to some power p of the radial distance from the star, and then we normalize the total mass of the disk to be equal to M_disk.

3.3. Fitting Process

Fitting the parameters listed above to the observational data was done using gradient descent, once we chose an error function. If s_i is the step size by which we vary parameter i, and the current value of parameter i is c_i, we simply compute an SED each for

• Our base set of parameters, with no modifications

• Our base set of parameters, with i modified by plus or minus s_i

Once we calculate the error for each resulting SED for parameter i, we identify the value of i for which we obtain the least error and set the parameter to this new base value. During each “iteration” of our code, this is done for every parameter that is being optimized. All step sizes are multiplied by a constant of our choosing (0.85) once each iteration is completed. The code terminates once the parameters being optimized stop being updated to new values that lower the error.

If the i’th model data point is m_i and the i’th observational data point is o_i, then for our error function, we calculate the expression shown below for wavelengths 1) below 45 microns 2) above 45 and below 80 microns and 3) above 80 microns.

$ϵ={\displaystyle \sum {\left(\frac{m_i-o_i}{\sigma }\right)}^2}$ (6)

Our error function takes the mean error per data point for each of these 3 regions and returns their sum. This is the error value our gradient descent code minimized.

4. Results

A best fit model is achieved in approximately 100 iterations. The fitting is terminated when variations do not improve the fit independently in the near-, mid-, and far-infrared/submillimeter. At termination, the fractional change in the error estimate is Δε/ε ~ 8 × 10⁻³. The model reaches the same best fit independent of the initial conditions, confirming that the appropriate minimum has been found. The resulting model parameters for the best fit are shown in Table 1.

Table 1. Final model parameters leading to a best fit of the SED of HD 100546. The dust weight parameters have been modified so that they are the percentage of dust (by mass) of that type.

Parameter Name	Varied/Fit?	Best Fit Value
Star Mass (Msun)	No	2.2
Star Radius (Rsun)	No	1.6
Star Temperature (Tsun)	No	1.8
Dust-To-Gas Mass Ratio	No	0.01
Surface Density Power Law	Yes	−1.70
Mean Molecular Gas Mass (amu)	No	2.3
Disk Mass (Gas + Dust) (Msun)	Yes	0.013
Disk Inner Radius (AU)	Yes	22
Disk Outer Radius (AU)	Yes	360
Flaring Factor (0.1 μm dust)	Yes	0.23
Flaring Factor (1 μm dust)	Yes	3.1
Flaring Factor (10 μm dust)	Yes	1.7
Flaring Factor (100 μm dust)	Yes	1.2
Flaring Factor (OH5 dust)	Yes	2.2
Dust Weight (Region 1, 0.1 μm)	Yes	17
Dust Weight (Region 1, 1 μm)	Yes	28
Dust Weight (Region 1, 10 μm)	Yes	7.9
Dust Weight (Region 1, 100 μm)	Yes	25
Dust Weight (Region 1, OH5)	Yes	22
Dust Weight (Region 2, 0.1 μm)	Yes	34
Dust Weight (Region 2, 1 μm)	Yes	7.4
Dust Weight (Region 2, 10 μm)	Yes	14
Dust Weight (Region 2, 100 μm)	Yes	40
Dust Weight (Region 2, OH5)	Yes	5.2
Dust Weight (Region 3, 0.1 μm)	Yes	33
Dust Weight (Region 3, 1 μm)	Yes	7.2
Dust Weight (Region 3, 10 μm)	Yes	2.3
Dust Weight (Region 3, 100 μm)	Yes	43
Dust Weight (Region 3, OH5)	Yes	14

The gross source structure properties are generally consistent with previous work. We find a total mass of 1.3 × 10⁻² M_sun, leading to a dust mass of 1.3 × 10⁻⁴ M_sun. In comparison, previous models find a value of M_dust ~ 10⁻⁴ M_sun [15]. We find an outer disk radius of 357 AU, consistent with the 200 - 400 AU [16] [17] and 350 AU [1]. The inner radius of our disk of 21 AU is about 50% larger than that adopted by others [15]. The surface density power law determined here of −1.7 is larger than the value of −1 adopted by previous authors [1] [15]. Clearly the gross structure is similar. The difference in the dust density power law is due partially to our choice to fit this parameter (rather than adopt a fixed value), and different choices about the grain opacities including a greater emissivity at longer wavelengths, thus requiring less dust at large distances where the low temperatures give rise to the long wavelength emission. It is heartening that the total mass of the disk is so similar, due to the fact that at long enough wavelengths, the emission should be roughly proportional to the mass of the dust. As a result, our model is consistent with the same physical constraints as other models (in particular the total mass of the dust), but with a different distribution of that dust. The difference in inner radius is discussed later in Section 5.

The fit of the SED is shown in Figure 2. Note that the model (line) fits the observed SED quite well over the vast majority of the observations. For most datapoints, the agreement is to much better than the uncertainty in the data. This includes both the slope in the mid-infrared, the knee around 90 - 100 μm, and the slope in the far-infrared and submillimeter. This suggests that the temperature and density distribution in the mid and outer disk are generally good, and that the adopted dust properties and distributions are generally good. There is some deviation in the 3 - 5 μm range. Emission at these wavelengths corresponds to warm-to-hot dust in the 600 - 1000 K region. This is dust that would sample the inner region of the disk, where we have a cavity. The potential implications of this discrepancy are discussed in Section 5.

5. Discussion

5.1. Overall

As noted earlier, many of the model properties are consistent with previous models of this source. In particular, the total dust mass and outer radius match previous models quite well, suggesting that the total amount of dust is consistent with previous work. However, the distribution of that dust is somewhat different.

Figure 2. The SED (datapoints) and best fit (red line). The agreement is generally quite good. In the NIR region, the fit is often within the uncertainty, except for a small region near 3 - 5 microns. Differences may be due to irregularity in the interior.

5.2. Dust Properties and Dust Settling

One significant difference between this work and previous work is the allowance for the variation in dust properties and sizes with radial position in the disk, as well as the allowance for different distributions as one moves away from the midplane. The different distributions as a function of z in the disk is represented by different flaring factors for the assumed different grain sizes. The flaring factors and thus scale heights of the 1 μm, 10 μm, and 100 μm grains all decrease as the grain size increases. The flaring factor for 1 μm grains is 3.1, for 10 μm grains is 1.7, and for 100 μm grains is 1.2. The flaring factor variation means that the largest grains are most concentrated toward the disk midplane, which is a clear sign of settling of the dust as would be expected as the large grains are less coupled to the gas. This is a key component of the core accretion model for planet formation.

On the other hand, the relative mass abundance of these grain sizes also varies with radial position. The dust weight for the 100 μm grains increases with radial position. Commensurately, the dust weight for the 1 μm and 10 μm grains decreases with increasing radial position. This would be consistent with grain growth. As the grains settle toward the disk midplane, they reach higher density environments. At these higher densities, collisions are more likely, which will deplete the smaller grain sizes as those grains grow into larger grains. In the outer regions this would also be enhanced because at low temperatures, icy grain mantles will help promote grain coagulation upon collision.

The 0.1 μm grains are, however, a different story. They are found preferentially near the midplane. There could be two reasons for this. First, in the outer part of the disk, where settling and collisions occur, there will be some amount of grain shattering. For the midplane densities expected, relative speeds of 0.1 - 10 m/s would be expected. At these speeds, there will be some amount of shattering of grains into smaller pieces [17] [18]. As a result, as grains actually enter the midplane, we would expect an increase in the number of small grains. This is exactly consistent with the surprisingly large dust weights for the smallest (0.1 μm) grains, especially in the outer part of the disk. Even moreso, it explains why the flaring factor for these smallest grains is so small, because they are preferentially formed and found at the disk midplane.

Finaly, the increase in smaller grain/larger grain mass abundances as one moves toward the inner region of the disk is also explicable. As gas falls in from the outer disk to the inner disk, it drags dust grains with it. The smaller dust grains are more strongly coupled to the gas, while the larger dust grains are more uncoupled. As a result, as the gas infalls, it is more likely to bring smaller dust grains with it, and leave larger dust grains behind. This will produce a larger ratio of small grains/larger grains in the inner disk as our model shows.

The assumption of discrete grain sizes in the models could potentially

5.3. Inner Edge Irregularity—A Potential Signature of Inner Disk Disruption

The less than great fit of the model to the observations in the 3 - 5 μm range of the SED is curious. Assuming a grain emissivity of something like Q_n = Q₀ (n/n₀)^p where 0 < p < 2 in the infrared yields a greybody peak of l_max T ~ 2400 μm K for p = 1 in the near-infrared. As a result, emission in the 3 - 5 μm range would correspond to dust in the 500 - 800 K range for smaller grains, and 600 - 1000 K for large grains. The model underprediction of flux in this wavelength range then is a result of a lack of dust in the 500 - 1000 K temperature range, as opposed to a modeling error. For unattenuated radiation, this corresponds to regions at r ~ 1 AU or less.

As a result, it appears that our model is lacking sufficient dust mass in the very interior. However, by the best fit modeling process, there is no smooth distribution of dust mass and grain sizes that will reproduce the 3 - 5 μm emission without simultaneously reducing the quality of the fit in other parts of the spectrum. The conclusion is that there is most probably irregularly distributed material in the inner 1 AU or so, with a relatively optically thin region out to the “inner edge” of the disk. It is worthwhile to note that the total mass of this missing dust is much smaller than the mass of the disk, and thus does not impact the other results of this study.

Such a structure is naturally consistent with a disrupting/dispersing disk inner edge. In this scenario, a process like a stellar wind or FUV scouring and evaporation or gas depletion by planet formation would deplete this region. It would also lead to an irregular distribution of material. Both of these are consistent with our model results.

5.4. Timing of Disk Processes

As seen above, we find evidence for ongoing settling of large grains in the upper portions of the outer disk, and production of smaller grains via grain shattering in the midplane of the outer disk. At the same time, we find potential evidence for an irregular distribution of smaller dust grains at low density very near the central source. If both of these are correct, it points out the fact that disk disruption in the center of the disk and grain evolution in the outer part of the disk can occur at the same time. As a consequence, this implies that the termination of disk evolution via disk disruption is not an instantaneous process that occurs simultaneously throughout the disk at a single given instant, but instead is a more complicated process that evolves in both space and time.

5.5. Role of Assumptions and Uncertainties on Results

The derived parameters of the best-fit model are well constrained within the limits set by the data and by the model assumptions. The principal uncertainties are +/−7% in the surface-density power-law index, +/−7% in total disk mass, and +/−10% and +/−20% - 50% in the inner and outer radii, respectively. The flaring factors for individual grain sizes are uncertain by roughly 30% - 50%, and the dust-weight fractions vary from ±10% to within a factor of two depending on grain size and radial zone. These variations modestly affect the detailed shape of the SED but do not change the overall interpretation.

Some model assumptions also contribute to systematic uncertainty. The adoption of a finite number of discrete grain sizes simplifies the underlying size distribution and may under-represent the continuum of grain properties that exists in reality; however, the chosen set (0.1 - 100 µm plus OH5) spans the range that dominates the opacity between 3 µm and 1 mm. The assumption of axisymmetry neglects possible azimuthal asymmetries seen in high-resolution ALMA images, but these affect only small portions of the disk and have minimal impact on the integrated SED. Vertical settling is parameterized through flaring factors rather than computed from full hydrodynamics, so the absolute scale heights of the different grain populations should be regarded as indicative rather than unique. Other sources of error include the choice of dust optical constants, the fixed gas-to-dust mass ratio, and the assumption that the disk is in radiative equilibrium with no transient heating. Within these limits, the overall structure inferred from the model—evidence for grain settling, shattering, and inner-disk depletion—remains robust, and the uncertainties quoted above adequately encompass the plausible range of physical conditions consistent with the observed SED.

6. Conclusion

We have modeled the emission from the disk around HD 100546. Using a self-consistent radiative transfer model and an unbiased automated fitting process, we have fit the observed spectrum and constrained disk parameters including the dust composition and distribution. We find evidence for grain settling and shattering in the outer disk, consistent with grain evolution in this region. We find evidence for grain segregation in the inner disk, consistent with infall of gas toward the central source. And we find inferential evidence of a tenuous, irregular grain population in the inner ~1 AU, consistent with disk disruption in the inner part of the disk. Taken together, these support a scenario whereby disk evolution and disruption is not an instantaneous process throughout the disk, but that it is a process which evolves at different times and over different timescales throughout the disk.

Acknowledgements

The authors would like to thank The Anderson Endowment for financial support to complete this work.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Benisty, M., Tatulli, E., Ménard, F. and Swain, M.R. (2010) The Complex Structure of the Disk around HD 100546: The Inner Few Astronomical Units. Astronomy and Astrophysics, 511, A75.[CrossRef]
[2]	Sissa, E. (2018) High-Contrast Study of the Candidate Planets and Protoplanetary Disk around HD 100546. Astronomy and Astrophysics, 619, A160.
[3]	Sloan, G.C., Kraemer, K.E., Price, S.D. and Shipman, R.F. (2003) A Uniform Database of 2.4–45.4 Micron Spectra from the infrared Space Observatory Short Wavelength Spectrometer. The Astrophysical Journal Supplement Series, 147, 379-401.[CrossRef]
[4]	Gry, C., Swinyard, B., Harwood, A., Trams, N., Leeks, S., Lim, T., et al. (2003) The ISO Handbook. European Space Agency.
[5]	Herschel Science Center (2018) The SPIRE Handbook: Herschel Explanatory Supplement, Volume IV. European Space Agency. https://www.cosmos.esa.int/web/herschel/legacy-documentation-spire
[6]	Pineda, J.E., Szulágyi, J., Quanz, S.P., van Dishoeck, E.F., Garufi, A., Meru, F., et al. (2019) High-Resolution ALMA Observations of HD 100546: Asymmetric Circumstellar Ring and Circumplanetary Disk Upper Limits. The Astrophysical Journal, 871, 48.[CrossRef]
[7]	Fedele, D., Toci, C., Maud, L. and Lodato, G. (2021) ALMA 870 μm Continuum Observations of HD 100546: Evidence of a Giant Planet on a Wide Orbit. Astronomy & Astrophysics, 651, A90.[CrossRef]
[8]	Casassus, S., Cárcamo, M., Hales, A., Weber, P. and Dent, B. (2022) The Doppler Flip in HD 100546 as a Disk Eruption: The Elephant in the Room of Kinematic Protoplanet Searches. The Astrophysical Journal Letters, 933, L4.[CrossRef]
[9]	Dullemond, C.P., Jhuasz, A., Pohl, A., Sershti, F., Shetty, R., Peters, T., Commercon, B. and Flock, M. (2012) Astrophysics Source Code Library. Astrophysics Source Code Library.
[10]	Mendigutía, I., Oudmaijer, R.D., Garufi, A., Lumsden, S.L., Huélamo, N., Cheetham, A., et al. (2017) The Protoplanetary System HD 100546 in Hαpolarized Light from Sphere/Zimpol: A Bar-Like Structure across the Di Astronomy & Astrophysics, 608, A104.[CrossRef]
[11]	Dorschner, J., Begemann, B., Henning, T., Jaeger, C. and Mutschke, H. (1995) Steps toward Interstellar Silicate Minerology. II. Study of Mg-Fe-Silicate Glasses of Variable Composition. Astronomy and Astrophysics, 300, Article 503.
[12]	Ossenkopf, V. and Henning, T.H. (1994) Dust Opacities for Protostellar Cores. Astronomy and Astrophysics, 291, Article 943.
[13]	Bruderer, S., van Dishoeck, E.F., Doty, S.D. and Herczeg, G.J. (2012) The Warm Gas Atmosphere of the HD 100546 Disk Seen by Herschel: Evidence of a Gas-Rich, Carbon-Poor Atmosphere? Astronomy & Astrophysics, 541, A91.[CrossRef]
[14]	Chiang, E.I. and Goldreich, P. (1997) Spectral Energy Distributions of T Tauri Stars with Passive Circumstellar Disks. The Astrophysical Journal, 490, 368-376.[CrossRef]
[15]	Pantin, E., Waelkens, C. and Lagage, P.O. (2000) ADONIS Observations of the HD 100546 Circumstellar Dust Disk. Astronomy and Astrophysics, 361, L9.
[16]	Panić, O., van Dishoeck, E.F., Hogerheijde, M.R., Belloche, A., Güsten, R., Boland, W., et al. (2010) Observations of Warm Molecular Gas and Kinematics in the Disc around HD 100546. Astronomy and Astrophysics, 519, A110.[CrossRef]
[17]	Blum, J. (2010) Dust Growth in Protoplanetary Disks—A Comprehensive Experimental/Theoretical Approach. Research in Astronomy and Astrophysics, 10, 1199-1214.[CrossRef]
[18]	Savvidou, S., Bitsch, B. and Lambrechts, M. (2020) Influence of Grain Growth on the Thermal Structure of Protoplanetary Discs. Astronomy & Astrophysics, 640, A63.[CrossRef]