© M. Ackermann Pyerin et al. 2021

1 Introduction

The formation of planets in protoplanetary disks is one of the main subjects of study in modern astrophysics. A very intriguing part of this topic is how newborn planets interact and actively shape the disk they are embedded in. This idea of planet-induced substructures is currently driven by recent observations with advanced telescopes, which have revealed that circumstellar disks can exhibit various prominent substructures, such as rings, gaps, asymmetries, or spirals (e.g., ALMA Partnership 2015; Long et al. 2018; Andrews et al. 2018; Andrews 2020; Cieza et al. 2021). Some of these substructures are thought to be indeed created by planet-disk interactions, rendering the observed features indicators of planet formation (e.g., Pinilla et al. 2012, 2015a; Zhu & Stone 2014), and may be used to constrain the properties of an embedded planet candidate, such as the mass and the orbital separation (e.g., Rice et al. 2006).

The subject of this study is the circumstellar disk surrounding the Herbig star HD 100546. This target is found at a distance of 108.1 ± 0.5 pc (Brown et al. 2021), with a right ascension of 11^h 33^m 25.3^s and a declination of −70° 11′ 41.2″ (Brown et al. 2021). Observations at multiple different wavelengths detected this disk, which spans 390 ± 20 au outward from the star in CO emission and 230 ± 20 au in continuum emission (Walsh et al. 2014). The following substructures were revealed: a ring of emission with a radius of 27 au (Hubble Telescope, Ardila et al. 2007), and a second very faint ring at 200 au (Atacama Large Millimeter/sub-millimeter Array (ALMA), Walsh et al. 2014).

While several other mechanisms have been proposed to produce ring-like substructures in protoplanetary disks, such as the dead-zone outer edge mechanism (Pinilla et al. 2016), particle growth by condensation near ice lines (Zhang et al. 2016; Stammler et al. 2017), and magnetic disk winds (Suriano et al. 2019), in this work we assume the planet-disk interaction to be the main driver for the observed substructures of HD 100546. In this mechanism, an embedded planet creates a pressure trap outward of its orbit that triggers an accumulation of dust by preventing the material from moving inward. This has been motivated in previous works, where simulations were executed to prove the similarity between structures created by planets and the observed structures of HD 100546 (Walsh et al. 2014; Montesinos et al. 2015; Quanz et al. 2015; Pinilla et al. 2015b; Fedele et al. 2021). For example, Pinilla et al. (2015b) estimated that a model with two high-mass planets can qualitatively reproduce the two rings observed in this disk. The inner planet is located at 10 au and the outer planet at 70 au. Furthermore, planet candidates were previously identified in this system by direct imaging (Quanz et al. 2013, 2015), but could not be confirmed because of confusion with disk features (Currie et al. 2014, 2015).

In this study, our aim is to further constrain the properties of the potential embedded planets in the HD 100546 circumstellar disk. Particularly, we wish to determine whether one or two planets are required to explain the observed substructures and to obtain their masses and distances from the central star. To do this, we calibrated and analyzed recent ALMA observations of the system to identify the main features and focused on resolving the outer ring that was suggested in Walsh et al. (2014); Pinilla et al. (2015b). On the numerical side, we carried out hydrodynamic multifluid simulations (gas and dust) with FARGO3D (Benitez-Llambay & Masset 2016; Masset 2000), and then generated synthetic observations from the outputs using RADMC3D (Dullemond et al. 2012). In this way, we tested different planet setup scenarios until we reached the best possible agreement between the synthetic and real observations by comparing their respective images and radial profiles.

The layout of this paper is as follows. In Sect. 2, we present the self-calibrated 0.9 mm ALMA observations of HD 100546. In Sect. 3, we describe the method of our hydrodynamic models and planet setups and also describe how synthetic observations can be obtained from them. In Sect. 4, we present the synthetic observations that we compare with the ALMA observations. In Sect. 5, we discuss what can be inferred about the HD 100546 system from these comparisons. Finally, Sect. 6 summarizes our main findings.

2 ALMA observations

2.1 Calibration

The data sets studied in this work include ALMA observations of HD 100546 at 0.9 mm wavelength (ALMA Band 7) under ALMA projects 2016.1.00497.S. (PI: A. Pohl) and 2015.1.00806.S (PI: J. Pineda), to which we refer as short-baseline dataset (SB) and long-baseline dataset (LB), respectively. For the SB dataset, the project aimed to measure linear polarization in the emission, and therefore the correlator was configured to observe dust continuum in its four spectral windows centered at 336.495 GHz, 338.432 GHz, 348.494 GHz, and 350.494 GHz, each with a bandwidth of 2 GHz. The polarization setup of this project was tuned to measure the polarizations (XX, XY, YX, and YY), but we only used the total intensity information. For the LB dataset, the correlator was configured to observe two spectral windows with dust continuum at 331.988 GHz and 343.488 GHz, and the two remaining windows were aimed at the molecular lines ¹²CO (J : 3−2) and ¹³CO (J : 3−2). The bandwidth of the dust continuum spectral windows in the LB dataset was also 2 GHz and 0.469 GHz for the remaining windows. The details of all the data are summarized in Table A.1.

After ALMA standard pipeline calibration, we used CASA 5.6.2 to handle and self-calibrate the datasets. In order to extract the dust continuum emission from the spectral windows with molecular line emission, we flagged the channels located at ± 25 km s⁻¹ from each targeted spectral line. Similarly to the calibration of the DSHARP sources (see Andrews et al. 2018), the remaining channels from all spectral windows were averaged into 125 MHz channels.

The self-calibration of the SB data was made in two stages: First, we self-calibrated the shortest baseline observations, taken in April 2017, with solution intervals of 360 s, 150 s, 60 s, and 24 s for the phase calibration and 150 s for the amplitude calibration. This calibrated dataset was joined with the remaining observations from SB taken on October 2016 and was again self-calibrated by repeating the same solution intervals. The LB data were self-calibrated independently from the SB data, with solution intervals of 360 s and 150 s for the phase calibration and 360 s for the amplitude.

Before combining the datasets, we corrected their phase centers such that the center of all observations was the center of the LB dataset. After self-calibration, we reduced the data volume by averaging the continuum emission into one channel per spectral window and 30 s of time binning. Finally, we applied the JvM correction to the images generated from the self-calibrated datasets. This correction accounts for the volume ratio ϵ_v between thepoint spread function of the images and the restored Gaussian of the CLEAN beam, as described in Jorsater & van Moorsel (1995) and also used in Andrews et al. (2021). We find ϵ_v = 0.69 for the SB and LB images and ϵ_v = 0.88 for the images with SB and LB combined.

2.2 Observational results: Outer ring

The inclusion of the different baselines allowed us to cover and resolve all the relevant scales of the disk. As shown in Fig. 1, we detect the bright compact disk around HD100546 (previously published in Pérez et al. 2020), as well as an extended faint ring that is only detectable when the SB information is included in the image reconstruction. When the disk is deprojected using the geometry from Casassus & Pérez (2019), with an inclination of 45° and a position angle (PA) of 150°, we find that the outer disk peaks at ≈ 1.85″ ± 0.1″ from the disk center, which is 200 ± 11 au at the distance of the source. Alternatively, we also imaged the outer structure after subtracting the contribution to the visibilities from the inner disk emission, in an attempt to reduce the high brightness dynamic range between the rings. This image is shown in panel c of Fig. 1.

To measure the flux from the outer ring, we calculated the flux density from the image generated using only the SB dataset with an elliptical mask with the inclination and position angle of the source, and we integrated the emission between 1.1″–2.5″, which effectively masked the inner disk emission. We detect the ring to be 4.5 ± 0.2 mJy, which does not include the 10% ALMA flux uncertainty.

We also estimated the optical depth τ of the emission to estimate the dust mass content in this region. We assumed that τ = − ln(1 − T_bright∕T_phys), where T_bright is the brightness temperature of a blackbody, and T_phys is the physical temperature of the midplane, which we assumed to be T_phys = 20 K constant. The temperature profile of the outer ring peaks at 0.1 K, and this very low temperature translates into an equally low estimated optical depth of τ_peak = 0.005. This is highly uncertain, however, because the underlying estimate of T_phys has a high level of uncertainty.

The dust mass of the model was estimated by assuming optically thin emission and assuming that the flux (F_ν) received at 0.9 mm is being emitted by dust with a constant temperature of T = 20 K (as e.g., Ansdell et al. 2016; Pinilla et al. 2018). We followed Hildebrand (1983), $$M_{\text{dust}}=\frac{d^2{F}_{\nu }}{{\kappa }_{\nu }{B}_{\nu }(T)},$$(1)

where d is the distance to the star, ν is the observed frequency, B_ν is the Planck function, and ${\kappa }_{\nu }=2.3\times {(\nu /230\text{GHz})}^{0.4}{\text{cm}}^2{\text{g}}^{-1}$ is the frequency-dependent mass absorption coefficient (as in Andrews et al. 2013).

By using this F_ν = 4.5 mJy for the faint outer ring in Eq. (1), we obtained a dust mass of 0.7 ± 0.07 M_⊕, where we assumed a conservative uncertainty of 10%. The error ofthis estimate is, in addition to the uncertainty of the ALMA flux, mainly caused by the mass absorption coefficient, because thisparameter depends on the dust grain size distribution and composition, which are unknown for protoplanetary disks.

Fig. 1 HD 100546 protoplanetary disk dust continuum emission at λ = 0.9 mm. Panels a and b are the images generated using only the SB and LB datasets, respectively, with SB being on asinh scale and saturated to show the faint outer ring. Panel c shows the outer ring structure when the SB and LB datasets are combined after subtraction of the visibilities from the inner disk. The scale bar represents 20 au at the distance of the source, the lower left ellipse represents the beam size of each image, and the numbers to the right of the ellipse indicate thebeam size in milliarcseconds. Panels a and b were generated with a robust parameter of 0.5, and panel c was generated with a robust parameter of 1.8.

3 Numerical simulations

In this section, we describe our numerical approach and the underlying assumptions we used to model the HD 100546 system. First, we derive estimates for the locations and masses of the potential embedded planets in Sect. 3.1. Second, we describe ourhydrodynamic simulations in Sect. 3.2, with the two initial conditions we explored and our modifications of the code to manually taper the mass of the planets over time (Sects. 3.2.1 and 3.2.2, respectively). Third, we describe how synthetic observations can be obtained from the multifluid FARGO3D simulations in Sect. 3.3. Finally, we outline our exploration of the parameter space in Sect. 3.4.

3.1 Planet parameters

3.1.1 Planet locations

We deduce the orbit of the planet candidates by analyzing the radial profile of the combined ALMA image. To do this, we need to first understand how a massive planet affects the disk emission.

If a planet is massive enough, it can open a gap in the gas distribution along its orbit. This creates a dust trap at the outer edge of the gap, where inward-drifting dust grains accumulate. This accumulation can then be observed as a ring (e.g., Pinilla et al. 2015b). This ring is located outward of the planetary orbit. Additionally, the gas surface density decreases toward the position of the planet, and hence it is expected that the planetary orbit should correspond to the local minimum in the radial profile of the observed flux.

From this assumption and the radial profile of our combined ALMA observations (Fig. 2), we find two promising locations forembedded planets: (1) The first location is such that it can explain the prominent inner substructure peaking at 27 au. Because we were unable to resolve the corresponding local minimum, we reverted to the estimated planet position from previous studies, specifically Pinilla et al. (2015b); Walsh et al. (2014), suggesting that the inner planet is located at 13 au. (2) Our combined ALMA image also reveals a ring with a peak at radius 200 au, which suggests that another planet is needed to explain the observed emission. We estimate the orbital separation of the outer planet to be at 143 au, corresponding to the local minimum in the observed intensity profile.

It is important to note that the assumption that exactly one planet is responsible for each gap is a simplification. Multiple planets withorbits close to each other could open a shared gap (e.g., as examined in the PDS70 disk, Bae et al. 2019). Furthermore, the behavior of opening a gap along the orbit and one ring being located outwards per planet is motivated by simulations of this specific disk. However, under different circumstances, a more complex behavior such as multiple rings induced by one planet and the planet positions coinciding with a dust ring has been found (Dong et al. 2017, 2018).

Fig. 2 Radial profile corresponding to the ALMA observation shown in Fig. 1. The gray interval marks the 1σ error obtained from statistical analysis of deviations in the image.

3.1.2 Planet mass estimates

The planet mass has a strong effect on the width of the induced gap. The precise relation of gap width and planet mass was studied by Dodson-Robinson & Salyk (2011); Pinilla et al. (2012); Rosotti et al. (2016); Fung & Chiang (2016); Facchini et al. (2018), yielding the dependence $$M_K=3M_{\star }{\left(\frac{D}{Kr}\right)}^3,$$(2)

where M_K is the planet mass estimate, which depends on the dimensionless parameter K, the stellar mass M_⋆, the orbit of the planet r, and the gap width in the millimeter dust continuum D (defined following Lodato et al. 2019 as the minimum to peak radial distance), obtained from observations. K describes how many Hill radii $R_{\text{Hill}}=r\sqrt[3]{M/3M_{\star }}$ of the planet correspond to the continuum gap width. K has been found to lie within the interval 7 ≤ K ≤ 10. This allows us to find mass estimates directly from the radial profile of the dust continuum emission (Fig. 2). We inferred the continuum gap width as the radial distance of the local flux peak of the substructure to the corresponding estimate of the planet position. This yields D₁ = 14 au for the inner and D₂ = 79 au for the outer gap, which corresponds to the mass estimates 8 M_Jup ≤ M_K,1 ≤ 24 M_Jup for the inner and 1.1 M_Jup ≤ M_K,2 ≤ 3.3 M_Jup for the outer planet.

We also took the minimum planet mass M_min into account that is required to open a gap at all, following the studies of Crida et al. (2006) and Lodato et al. (2019). First, there is the model of Crida et al. (2006). This model is a balance between the pressure and viscosity torque from the disk with the gravitational torque from the planet. The authors estimated the minimum mass to be $$M_{\text{min}}^{\text{Crida}}=\left\{q{M}_{\star }|\frac{3\sqrt[3]{3}H}{4r}{q}^{-1/3}+\frac{50\alpha c_{\text{s}}H}{\sqrt{\frac{G{M}_{\star }}{r^3}}{r}^2}{q}^{-1}=1\right\},$$(3)

where $c_{\text{s}}=\sqrt{\frac{k_{\text{B}}T}{\mu m_H}}$ is the isothermal sound speed, $H=\frac{c_{\text{s}}}{{\Omega }_K}$ is the disk pressure scale height, q is the planet mass in units of the stellar mass, and α is the disk viscosity from the α-disk model (Shakura & Sunyaev 1973). Then there is the model of Lodato et al. (2019), which is based on hydrodynamical simulations that include dust dynamics. It finds a minimum mass of $$M_{min}^{\text{Lodato}}=\{\begin{array}{ll}0.3M_{\star }{\left(\frac{H}{r}\right)}^3\hfill & \text{for}\text{St}<1\hfill \\ 3M_{\star }{\left(\frac{H}{5.5r\sqrt{\text{St}}}\right)}^3\hfill & \text{for}\text{St}\ge 1,\hfill \end{array}$$(4)

where St corresponds to the Stokes number, defined by $$\text{St}=\frac{\pi a_{\text{dust}}{\rho }_{\text{bulk}}}{2{\Sigma }_{\text{g}}},$$(5)

with a_dust being the dust grain size, ρ_bulk the intrinsic dust density, and Σ_g the gas surface density.

From these two criteria, we estimated the minimum planet mass by choosing the maximum of the two models, $$M_{\text{min}}=\text{max}\left(M_{\text{min}}^{\text{Crida}},M_{min}^{\text{Lodato}}\right).$$(6)

For α = 10⁻⁴, we obtain M_min,1 = 0.3 M_Jup for the inner and M_min,2 = 1.7 M_Jup for the outer planet. For α = 10⁻³, the minimum mass estimates are M_min,1 = 0.7 M_Jup for the innerand M_min,2 = 3 M_Jup for the outer planet.

We took the lowest possible mass for each of the planets for our simulations because planets have not directly been observed in this disk, and thus a lower-mass companion is more likely to exist. Therefore we determine the planet masses M for the FARGO3D setups as $$M=\text{max}\left({\text{min}}_{\text{K}}\left(M_{\text{K}}\right),M_{\text{min}}\right),$$(7)

which returns an estimate of M₁ = 8 M_Jup for the inner planet mass. For the outer companion, we obtain different estimates depending on α: α = 10⁻⁴ leads to M₂ = 1.7 M_Jup, and α = 10⁻³ implies M₂ = 3 M_Jup.

3.2 Hydrodynamic simulations with FARGO3D

We modeled this disk with the numerical hydrodynamics solver FARGO3D (Benitez-Llambay & Masset 2016; Masset 2000; Benítez-Llambay et al. 2019; Weber et al. 2019), which is designed to simulate planet-disk interactions. This code evolves a numerical fluid over time that is described with a grid of densities, velocities, and energies.

We used a multifluid approach and model the gaseous dusty disk around HD 100546 as five fluids: one fluid represents the gas, and the four remaining fluids describe the dust, where each of them corresponds to a specific grain size $a_{\text{dust}}\in \left\{0.1\text{upmum},4.6\text{upmum},0.22\text{mm},\text{and}1\text{cm}\right\}$.

The FARGO3D simulations were performed in two dimensions on a polar grid (r, φ), extending azimuthally from 0 to 2π and radially from r_min to r_max. The radial extent was determined in order to enclose the orbits of all planet candidates as well as the induced substructures without any potential boundary effects. They were calculated based on the planet Hill radius R_H as r_min = (r₁ − 3R_H,1)∕3 = 2.6 au and r_max = 3(r₂ + 3R_H,2) = 500 au. For the radial simulation boundaries, we used a power-law extrapolation as the boundary condition for the densities and standard Keplerian speed/antisymmetric boundaries for the velocities in the r/φ direction (KEPLERIAN2DDENS and KEPLERIAN2DVAZIM/ANTISYMMETRIC in the FARGO3D setup, respectively). The grid resolution was chosen such that one grid cell corresponded to one-tenth of the pressure scale height H at the two planet locations r₁ = 13 au and r₂ = 143 au, which equates to an azimuthal resolution of n_ϕ = 1024 cells and a radial resolution of n_r = 540 cells.

To determine the duration of the simulations, we need to consider the age of the HD 100546 system, which has been estimated to be 4.8 Myr (Wichittanakom et al. 2020). However, a simulation of such a duration is computationally unfeasible because we need to choose a high grid resolution to resolve the disk at the locations of the two planet candidates. We instead ran our simulations for 0.74 Myr (corresponding to 2.3 × 10⁴ orbits of the inner planet or 6.3 × 10² orbits of theouter planet), assuming that, qualitatively, the disk reaches a quasi-steady state within this time. We discuss this assumption in Appendix B.

3.2.1 Initial conditions

To model the initial distribution of the gas surface density Σ_g of the disk for the FARGO3D simulations, we used two approaches. The first approach was a simple power-law profile, $${\Sigma }_{\text{g}}(r,\phi )={\Sigma }_0{\left(\frac{r}{r_1}\right)}^{-1},$$(8)

where r₁ = 13 au is the orbit of the innermost planet, which was used as a length scale by FARGO3D, and M_disk is the total gas disk mass enclosed between r_min and r_max. The normalization constant Σ₀ is obtained such that $M_{\text{disk}}={\displaystyle {\int }_{r_{\text{min}}}^{r_{\text{max}}}}2\pi r{\Sigma }_{\text{g}}(r)\text{d}r$, where we assumed the total gas disk mass to be M_disk = M_dust∕ϵ. Here, ϵ = 10⁻² corresponds to the dust-to-gas mass ratio, taken to be equal to the interstellar value, and M_dust = 66 M_⊕ is the total dust disk mass obtained from the total observed flux of our combined ALMA dust continuum emission of 435 mJy (see Eq. (1)).

The other initial condition used in this work is the self-similar solution from Lynden-Bell & Pringle (1974), $${\Sigma }_{\text{g}}(r,\phi )={\Sigma }_0{\left(\frac{r}{R_{\text{c}}}\right)}^{-1}\exp \left(-\frac{r}{R_{\text{c}}}\right).$$(9)

This profile is a power law combined with an exponential decay to reduce the amount of material in the outer disk. Here, R_c = 80 au is the critical radius, and Σ₀ is chosen analogously to the above.

The initial dust surface densities were set following Σ_d(a_dust) = f(a_dust)ϵΣ_g, such that the dust-to-gas ratio radial profile was radially constant and equal to ϵ. The factor f(a_dust) determines the mass distribution across the different dust sizes, and it was set such that the number density of the dust behaves as $n(a_{\text{dust}})\text{d}{a}_{\text{dust}}\propto a_{\text{dust}}^{-3.5}\text{d}{a}_{\text{dust}}$, following anMRN distribution (Mathis et al. 1977).

3.2.2 Mass-tapering model

To mimic the fact that planets grow over time and avoid any numerical artifacts right at the beginning of the simulations by introducing the whole planet mass, we used the following mass-tapering approach. This formula is a variation of the mass-tapering algorithm that is included in FARGO3D. We changed the code because we needed to be able to set a different taper for each planet individually in order to study the effect of planet formation timescales. We also added a new offset parameter. Our model reads $$\frac{m_i(t)}{M_i}=\{\begin{array}{ll}0\hfill & \text{for}t<{\Delta }_i\hfill \\ \frac{1}{2}\left(1-\cos \left(\pi \frac{t-{\Delta }_i}{{\tau }_i}\right)\right)\hfill & \text{for}{\Delta }_i\le t<{\Delta }_i+{\tau }_i\hfill \\ 1\hfill & \text{for}t\le {\Delta }_i+{\tau }_i,\hfill \end{array}$$(10)

where m_i(t) is the effective mass of the planet i at a given time t, and M_i is the final planetary mass. This model has two parameters, which we refer to as the “timing parameters”: (1) Δ_i, the time at which the planet is introduced into the simulation; and (2) τ_i, the formation timescale after the planet is introduced. We fiducially set the delays Δ₁ and Δ₂ to zero and the tapers τ₁ and τ₂ to 16 Kyr. We chose 16 Kyr, corresponding to 1000 orbits of theinner planet, as an arbitrary timescale that is small compared to the simulated duration, but still long enough to prevent numerical artifacts by introducing the full planet mass from the beginning of the simulation.

We note that while our approach does increase the planet mass over time, the actual physical accretion process is ignored in this present work, that is, there is no mass transfer from the dust or gas fluids onto the planet. We also fixed the planets orbit and neglected its potential migration.

3.3 Synthetic observations with RADMC-3D

To directly compare the density distributions from FARGO3D with our self-calibrated ALMA image, we performed a radiative transfer simulation and rendered a synthetic observation, using the code RADMC-3D (Dullemond et al. 2012). This requires expanding the two-dimensional dust surface density distributions into three dimensions, which we did by using a Gaussian distribution along the vertical (z) axis, assuming hydrostatic equilibrium in this direction, $${\rho }_{\text{d}}(r,\phi ,z)={\Sigma }_{\text{d}}(r,\phi )\frac{1}{\sqrt{2\pi }{H}_{\text{dust}}(r)}\exp \left(-\frac{z^2}{2H_{\text{dust}}{(r)}^2}\right),$$(11)

with the disk height model of Dullemond & Dominik (2004), $$H_{\text{dust}}(r)=\{\begin{array}{ll}\chi H(r)\hfill & \text{for}{a}_{\text{dust}}>1\text{μm}\hfill \\ H(r)\hfill & \text{otherwise},\hfill \end{array}$$(12)

where $H=\frac{c_{\text{s}}}{{\Omega }_K}$ is the disks pressure scale height, and χ = 3% is the scale height reduction factor encoding the effect of dust settling, which is more prominent for the larger grain sizes. It therefore results in a shorter dust scale height compared to the one for smaller grain sizes that roughly follow the gas.

We assumed the disk to be composed of silicates (Dorschner et al. 1995) and used corresponding optical constants from the Jena database. To match our combined ALMA observations, the wavelength was fixed at λ = 0.9 mm, and all synthetic observations are convolved with a 0.3″ Gaussian beam similar to the actual observations.

3.4 Parameter space exploration

In order to find a model that reproduces the ALMA observations, we ran multiple simulations with different parameters (see Table 1). We started by testing how models with a single planet compare to the observation, and ran a simulation for both the initial gas surface density profiles we introduced in Sect. 3.2.1 (model I with the power law, and model II with the self-similar solution). Then, we executed simulations with two planets for both initial conditions (model III with the power law, and model IV with the self-similar solution). In all of these four simulations, we introduced the planet(s) right from the beginning of the run and used the fiducial taper of τ = 16 Kyr as well as a fiducial viscosity of α = 10⁻⁴.

Then, we explored the effect of different planet formation timescales. Because the ALMA observation reveals that the inner ring is much brighter than the outer ring, the inner pressure trap must exist very early in order to trap as much material as possible. For this reason, we always introduced the inner planet immediately at the beginning of the simulations, that is, we always set the offset Δ₁ = 0 and chose a fiducial mass taper of τ₁ = 16 Kyr, and we only varied the timing parameters for the second planet. In model V we set Δ₂ = 6 × 10⁵ yr while keeping τ₂ = 16 Kyr, while in model VI, we kept Δ₂ = 0 and set τ₂ = 6 × 10⁵ yr. This means that in models V and VI, the outer planet reaches its final mass on a timescale comparable to the duration of the entire simulation (7.4 × 10⁵ yr). In these simulations,the viscosity remained the same as for the previous ones for comparability (α = 10⁻⁴), and the initial profile is the self-similar solution in both of them because we show below that such a profile is required to reach the best agreement between synthetic and real observations.

Finally, we simulated a disk with a different viscosity α. We chose a higher viscosity of α = 10⁻³ in model VII. This simulation used the fiducial timing setup, that is, it introduced both planets from the start of the simulation with the fiducial taper (Δ₁ = Δ₂ = 0, τ₁ = τ₂ = 16 Kyr) and used the self-similar solution initial profile. It therefore directly compares to model IV. In this model, the mass estimate for the outer planet increases to 3 M_Jup.

The parameters that are common in all simulations are listed in Table 2. An overview of the exact values of the parameters that differed between the seven different simulations is given in Table 1.

4 Results

The main results we obtained for the different planet setups and disk parameters are shown in Figs. 3–6. For each setup, we present the corresponding synthetic observation and the combined ALMA image, both normalized to their respective maximum, and using the same color map. In addition, we display the azimuthally averaged normalized radial profiles as a metric to compare the two images quantitatively.

Figure 3 shows models I and II, which solely include the inner planet, with different initial density profiles (Eqs. (8) and (9), respectively). Both of these models induce a ring in the synthetic observation that matches the inner ring of the real observation in the image and the radial profile. However, these models do not create an outer ring because there is no second planet that might induce an outer trap. Comparison of the radial profiles reveals how strongly the synthetic observations deviate from the observation. Model I reaches its largest deviation at 110 au with 133σ, while it is 4σ at 200 au for model II.

In Fig. 4, we show simulations including two planets (models III and IV) with different initial density models. In contrast to the one-planet models (I and II), both of these models create two rings. The location of these rings matches the ALMA observation, and the radial profiles show a similar shape. We therefore continued our parameter search with two planets.

In model III, the outer structures are relatively bright compared to the inner ring, which is strongly overestimated: in the model, the relative peak intensity of the outer ring is 86 times higher than in the ALMA observation. This issue does not arise with model IV, which reproduces the relative peak intensity of the two substructures properly. Model IV is within a 3σ range around the observation in the normalized radial profile. We therefore continued our parameter search using two planets and the self-similar solution as initial condition.

However, model IV presents an issue regarding the symmetry of the outer features: the ring is mainly concentrated in one point and is much fainter and thinner along the remaining structure in the azimuthal direction. This differs from the observed outer ring, which is evenly bright and symmetric. To resolve this issue, we conducted further tests to quantify thepotential effect of different model parameters, specifically, the timing parameters and the disk viscosity α.

Figure 5 shows the results with different planet-mass tapering. The models are identical to model IV, except for one of the parameters of the mass-tapering model. In model V, the outer planet is introduced Δ₂ = 6 × 10⁵ yr after the inner one. This reduces the relative strength of the outer substructure by ~ 80% compared to model IV. In model VI, the outer planet mass is tapered over τ₂ = 6 × 10⁵ yr. This reduces the outer rings flux by ~40% relative to the base model IV. Both of these models are within 3σ around the ALMA detection. However, the outer ring is not more symmetrical in either one of them.

Finally, we present the result for a higher disk viscosity of α = 10⁻³ in Fig. 6 (model VII). Motivated by the previous results, model VII used the self-similar solution initial profile and the fiducial mass tapering setup (Δ_1,2 = 0, τ_1,2 = 16 Kyr). This model produces an azimuthally symmetric ring. This can be explained as follows: While vortices can live long at the edge of a gap opened by a massive planet due to the Rossby wave instability in a disk with a low viscosity of α = 10⁻⁴, they dissipate earlier with the higher viscosity α = 10⁻³ (Lin & Papaloizou 2011). However, the minimum planetary mass needs to be increased in order to open a comparable gap in this higher-viscosity scenario. We find the mass estimate for the outer planet to be 3 M_Jup (see Eq. (7)). The mass of the inner planet remains unchanged (8 M_Jup).

The outer substructure in the synthetic observation of model VII agrees well with the ALMA observation because it does not display the asymmetry seen in previous models. Its normalized radial profile resides within a frame of 3σ around the detection. Consequently, model VII leads to the best agreement between our synthetic observation and the real one.

To mimic how ALMA would see the disk that model VII generates, we took corresponding synthetic observation given by RADMC-3D and used the package galario to calculate the visibilities with the same uv-coverage as the observations. Figure 6b shows the synthetic image obtained through radiative transfer and convolved with a Gaussian beam in the image space, as we did for all the other synthetic observation, and Fig. 6c shows the simulated observation using the visibilities, imaged with the same conditions as the observations. Both images are consistent with each another andreproduce the observed radial profile, proving the observability of the structures generated in our simulations.

Fig. 3 Results for the simulations with a single planet, positioned at 13 au with a mass of 8 MJup (position marked in panel d as a dashed vertical line). The two simulation results correspond to different initial density models. Simulation I started with a power law (Eq. (8)) for the radial fall-off of Σg∕d and simulation II with a self-similar solution profile (Eq. (9)). We present them as synthetic observations convolved with a 0.3″ beam (panels b and c) in direct comparison to the modeled ALMA image (panel a).

Fig. 4 Outcome of the simulations including two planets. We positioned the inner planet at 13 au with a mass of 8MJup and the outer planet at 143 au with 1.7 MJup (positions marked in panel d as dashed vertical lines). The planets form at the same time and equally fast, but the initial density profile is assumed differently in the two compared setups: the initial density profile for simulation III was a power-law, while simulation IV started as a self-similar solution profile. The inconsistent bump, shown as the red curve centered on the position of the outer planets, is an artifact of the lacking accretion process in the simulation, as can be concluded based on the results of Bergez-Casalou et al. (2020), who have shown the effect of accretion on the gap formation process. They found that when the accretion process is included in the simulation, the gap is cleared out completely for a disk setup like ours. The accumulation at the position of the planet would therefore not be observed in reality.

5 Discussion

From the results above, it is clear that two planets are needed to explain the newly resolved outer ring in addition to the inner feature. None of the one-planet models of the disk (I and II) reproduce any strong visible substructures that would resemble the appearance of the outer ring, while the two-planet models feature two rings (e.g., VII). This conclusion is confirmed quantitatively by analyzing the errors of the normalized radial profiles: the models with a single planet have an above-significant deviation (> 3σ) from the observation, while the models with two planets IV-VII have a deviation of less than 3σ in the radial profile. This finding is consistent with the results of Pinilla et al. (2015b); Fedele et al. (2021), who both concluded that the outer substructure is likely an indicator of a second companion.

Additionally, we find that the initial density profile needs to have a stronger radial fall-off than the commonly assumed power law (Eq. (8)) in order to model the flux ratio of the inner and outer ring correctly. We have achieved this by using the self-similar solution from Lynden-Bell & Pringle (1974) as initial condition.

Our study also shows that changing the formation timing of the outer planet does affect the resulting radial profile and synthetic image (models IV, V, and VI), but all three models still reside within a 3σ range around the observation in the normalized radial profile. This means that the outer planets mass of 1.7 M_Jup in model V, where the planet started to form very recently (140 Kyr before the end of the simulation), is sufficient to successfully induce the outer substructure. This is much lower than previously estimated: Pinilla et al. (2015b) estimated that if the outer planet is much younger than the inner planet (1 Myr younger), it should be more massive (≳15 M_Jup). This difference probably arises from the two different methods. Pinilla et al. (2015b) assumed an analytical solution to introduce the gaps in the gas surface density profile and ran dust evolution models considering that the gas density remains constant over time. In our current approach, however, we modeled gas and dust evolution simultaneously, but we neglected dust growth.

Our study favors a Shakura & Sunyaev (1973) disk viscosity of α ≳ 10⁻³ because we searched for an α that avoids asymmetric accumulations in the outer ring, which is successfully achieved in model VII (Fig. 6). Models of a disk with a lower viscosity (α = 10⁻⁴) induce an asymmetric outer feature with a strong density accumulation that is absent in the ALMA image, and they also feature an accumulation of material around the position of the outer planet (see Fig. 4c). However, from the study of Bergez-Casalou et al. (2020), we suspect that this accumulation is an artifact of not simulating the accretion of the surrounding material onto the planet. This poses the question whether the shape of the outer ring itself would be as strongly affected by including accretion in the simulation as the gap. If this were the case, a simulation with a low viscosity α = 10⁻⁴ might be able to produce a symmetrical feature with the outer planet mass estimate as low as 1.7 M_Jup. Further simulations including the planet accretion process are thus needed to address this question. This lies beyond the scope of this paper, however.

In a recent study, Fedele et al. (2021) examined the HD 100546 disk planet setup with smooth particle hydrodynamics (SPH) simulations. They qualitatively reproduced the disk appearance with an inner planet of final mass of 3.1 M_Jup at 15 au and an outer planet that reached8.5 M_Jup at 110 au. While the number of planets and the positions are consistent with what we derived from the radial profile in Sect. 3.1.1 (inner planet: 13 au; outer planet: 143 au), the mass estimates differ from our results (inner planet: 8 M_Jup; outer planet: 3 M_Jup for α = 10⁻³, or 1.7 M_Jup for α = 10⁻⁴). The origin of these differences is unclear. One possible explanation is that for the outer planet mass, Fedele et al. (2021) reported an estimate based on the gap width with the same model as we used (Eq. (2)), but fixed the proportionality constant to an average value while we chose the minimum possible mass; and they have a higher estimate of 100 au for the dust gap width than we do (79 au). As a consequence, their estimated planet mass (5.5 M_Jup) is a higher than our estimate (3 M_Jup).

The results we deduced from submillimeter ALMA observations compare well to findings from observations at other wavelengths. Brittain et al. (2014) hypothesized the presence of a planet in the HD 100546 system with an orbit of 15 AU from molecular emission line observations. This companion appears consistent with the inner planet of our best model (13 AU). Furthermore, Follette et al. (2017) analyzed scattered-light observations of HD 100546, revealing multiple spiral arms. They proposed that a planet of a few Jupiter masses at 100 AU could explain these observed spiral structures. This planet roughly matches the outer planet from our model (at 143 AU).

Fig. 5 Comparison of the effect of different planet formation timescales. The planet setup is the same for both: two planets, positioned at 13 au and 143 au, with a mass of 8 MJup and 1.7 MJup. In simulation V, the outer planet started to form much later than the inner planet (Δ1 = 0, Δ2 = 600 Kyr, τ1 = τ2 = 16 Kyr). For simulation VI, the planets started to form simultaneously, but the outer planet gained mass significantly slower (Δ1 = Δ2 = 0, τ1 = 16 Kyr, τ2 = 600 Kyr).

Fig. 6 Result for a disk viscosity of α = 10−3 (model setup VII). The mass estimates in this simulation are 8 MJup for the inner and 3 MJup for the outer planet. Because this setup shows the best agreement with the modeled ALMA image of all our results, we further simulated an observation with ALMA based on the unconvolved RADMC-3D output to prove the observability of the simulated substructures. This simulated observation is shown in panel c of this figure.

6 Conclusions

We aimed to constrain some of the properties of potential planets embedded in the HD100546 circumstellar disk from the observed substructures. In particular, we were interested in constraining their masses and orbits. To do this, we analyzed recent ALMA observations of the system, and then matched them as closely as possible with our numerical model, iteratively improving the fit by systematically varying disk and planet parameters. In our numerical model, we first simulated the gas and dust in a protoplanetary disk with the multifluid code FARGO3D, and then rendered a synthetic observation from the results with RADMC-3D.

Thereby, we showed that the recent ALMA observations of HD 100546 can be reproduced by a protoplanetary disk model hosting two companions orbiting at separations of 13 au and 143 au. The best agreement between synthetic and real observations is reached with model VII, which uses the self-similar solution as initial condition and value of α = 10⁻³. This model matches the observation within 3σ of certainty in the normalized radial profile, and the resulting outer ring does not show azimuthal asymmetries. From the observed intensity profile, we estimate a mass of 8 M_Jup for the inner planet. This agrees with the results of other studies (Pinilla et al. 2015b): ~ 10 M_Jup), but it is higher than suggested by Fedele et al. (2021) (3.1 M_Jup).

In contrast to previous results, we conclude that the outer planets mass can be as low as 3 M_Jup for a disk viscosity of about α = 10⁻³, and 1.7 M_Jup for a lower viscosity (α = 10⁻⁴). Our results show that this companion can successfully reproduce the outer substructure, both in a model were it formed early in the disks lifetime and in a model where it formed very recently (as young as 140 Kyr), equally well.This estimate for the outer planet mass is low in comparison to the previous studies of Pinilla et al. (2015b) (≳ 15 M_Jup for a young planet, and ≲ 5 M_Jup for an older planet) and Fedele et al. (2021) (8.5 M_Jup). We further find that to properly reproduce the brightness ratio of the inner and outer ring, the disks initial gas surface density profile needs to be modeled as a self-similar solution profile.

Finally, the planet mass estimates need to be evaluated in terms of whether companions this massive might be directly visually detected by current instrumentation capabilities. Pérez et al. (2020) have deduced a lower limit for this in ALMA observation, depending on the star age. They found that the planets of this system would have to be more massive than 33 M_Jup for an age of 4 Myr to be detected by ALMA. This means that the detection of the proposed planets with 8 M_Jup and 3 M_Jup is still beyond our current ALMA capabilities. In addition, they are very difficult to directly image with instruments like SPHERE at the VLT (Boccaletti et al. 2020). Future high-resolution observations are therefore needed to prove or rule out the existence of embedded planets in the HD 100546 protoplanetary disk.

Acknowledgements

We would like to thank Camille Bergez-Casalou, Cornelis Dullemond and Marcelo Barraza for helpful discussions and technical insights. The Computational resources of the Max Planck Computing and Data Facility made this project possible and are gratefully acknowledged. We are also thankful for the support provided by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. We further acknowledge the efforts of Adriana Pohl, Akimasa Kataoka, and Cornelis Dullemond to motivate and propose the observation of the HD 100546 system with ALMA. We thank Ruobing Dong for revising this manuscript and providing constructivefeedback and suggestions. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00497.S., and ADS/JAO.ALMA#2015.1.00806.S. A.L.M.A. is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

Appendix A Observation details

In table A.1 we present the details of the Atacama Large Millimeter/sub-millimeter Array (ALMA) observations that were used in this work.

Appendix B Convergence and resolution tests

Our FARGO3D simulations were run for a shorter time than the age of the HD 100546 system (simulation time: 0.74 Myr, estimated system age: 4.8 Myr, Wichittanakom et al. (2020)). This was done under the assumption that the disk model reaches a quasi-steady state within the simulated time. To confirm this assumption, we illustrate the temporal convergence behavior in the gap formation process of our model in Fig. B.1 and show that our gas surface density does not vary significantly from 200 Kyr.

Fig. B.1 Evolution of the azimuthally averaged surface density profile, for the dust fluid with a grain size of adust = 5μm, from simulation IV.

Furthermore, we carried out all of our FARGO3D simulations on a two-dimensional grid that resolved one-tenth of the pressure scale height. The consequent dimensions of the grid were n_r = 540 radial and n_ϕ = 1024 azimuthal cells. To show that this resolution is sufficient and does not affect the results, we present a comparison of model VII and a simulation with the identical setup, but resolving one-twentieth scale height (n_r = 1080, n_ϕ = 2048) in Fig. B.2. The result does not change significantly by the higher resolution. The strongest deviation arises with the nonphysical material accumulation at the location of the inner planet (13 au), where the peak is higher.

Fig. B.2 Radial profiles of the surface density results after 1000 inner planet orbits for the resolution used in all the simulations presented in the paper (nr = 540, nϕ = 1024) compared to a higher resolution (nr = 1080, nϕ = 2048). The results are for the dust fluid with a grain size of adust = 5μm with simulation setup VII in both cases.

References

All Tables

All Figures

Fig. 1 HD 100546 protoplanetary disk dust continuum emission at λ = 0.9 mm. Panels a and b are the images generated using only the SB and LB datasets, respectively, with SB being on asinh scale and saturated to show the faint outer ring. Panel c shows the outer ring structure when the SB and LB datasets are combined after subtraction of the visibilities from the inner disk. The scale bar represents 20 au at the distance of the source, the lower left ellipse represents the beam size of each image, and the numbers to the right of the ellipse indicate thebeam size in milliarcseconds. Panels a and b were generated with a robust parameter of 0.5, and panel c was generated with a robust parameter of 1.8.

In the text

	Fig. 2 Radial profile corresponding to the ALMA observation shown in Fig. 1. The gray interval marks the 1σ error obtained from statistical analysis of deviations in the image.
In the text

Fig. 3 Results for the simulations with a single planet, positioned at 13 au with a mass of 8 MJup (position marked in panel d as a dashed vertical line). The two simulation results correspond to different initial density models. Simulation I started with a power law (Eq. (8)) for the radial fall-off of Σg∕d and simulation II with a self-similar solution profile (Eq. (9)). We present them as synthetic observations convolved with a 0.3″ beam (panels b and c) in direct comparison to the modeled ALMA image (panel a).

In the text

Fig. 4 Outcome of the simulations including two planets. We positioned the inner planet at 13 au with a mass of 8MJup and the outer planet at 143 au with 1.7 MJup (positions marked in panel d as dashed vertical lines). The planets form at the same time and equally fast, but the initial density profile is assumed differently in the two compared setups: the initial density profile for simulation III was a power-law, while simulation IV started as a self-similar solution profile. The inconsistent bump, shown as the red curve centered on the position of the outer planets, is an artifact of the lacking accretion process in the simulation, as can be concluded based on the results of Bergez-Casalou et al. (2020), who have shown the effect of accretion on the gap formation process. They found that when the accretion process is included in the simulation, the gap is cleared out completely for a disk setup like ours. The accumulation at the position of the planet would therefore not be observed in reality.

In the text

Fig. 5 Comparison of the effect of different planet formation timescales. The planet setup is the same for both: two planets, positioned at 13 au and 143 au, with a mass of 8 MJup and 1.7 MJup. In simulation V, the outer planet started to form much later than the inner planet (Δ1 = 0, Δ2 = 600 Kyr, τ1 = τ2 = 16 Kyr). For simulation VI, the planets started to form simultaneously, but the outer planet gained mass significantly slower (Δ1 = Δ2 = 0, τ1 = 16 Kyr, τ2 = 600 Kyr).

In the text

Fig. 6 Result for a disk viscosity of α = 10−3 (model setup VII). The mass estimates in this simulation are 8 MJup for the inner and 3 MJup for the outer planet. Because this setup shows the best agreement with the modeled ALMA image of all our results, we further simulated an observation with ALMA based on the unconvolved RADMC-3D output to prove the observability of the simulated substructures. This simulated observation is shown in panel c of this figure.

In the text

	Fig. B.1 Evolution of the azimuthally averaged surface density profile, for the dust fluid with a grain size of adust = 5μm, from simulation IV.
In the text

	Fig. B.2 Radial profiles of the surface density results after 1000 inner planet orbits for the resolution used in all the simulations presented in the paper (nr = 540, nϕ = 1024) compared to a higher resolution (nr = 1080, nϕ = 2048). The results are for the dust fluid with a grain size of adust = 5μm with simulation setup VII in both cases.
In the text