© The Authors 2024

1 Introduction

A fundamental but still poorly constrained property of protoplanetary discs is their mass. This parameter is of paramount importance because it affects the planet formation process in many different ways. First of all, it determines the amount of material available for the formation of planets (Manara et al. 2018; Mulders et al. 2021), and, once formed, it can also influence their migration efficiency (e.g. Nelson 2018). Furthermore, it controls the interactions between gas and dust particles through drag, strongly influencing the dynamics of the dust, and therefore also its radial drift and growth (e.g. Birnstiel et al. 2010a,b; Laibe 2014; Gonzalez et al. 2017; Franceschi et al. 2022). The mass of the protopanetary disc can also affect the ionisation rate and thus the magnetic activity (e.g. Lesur et al. 2022). Finally, if the disc-to-star mass ratio is sufficiently high (i.e. q = M_disc/M_* ≳ 0.07–0.1, Kratter & Lodato 2016), the disc self-gravity (SG) affects the morphology and the stability of the system. This may potentially lead to the development of gravitational instability (GI) in the form of a spiral structure (Jeans 1902; Lin & Shu 1964; Paczynski 1978; Binney & Tremaine 1987; Lin & Pringle 1987; Bertin & Lin 1996; Kratter & Lodato 2016), which several authors claim to have observed (e.g. see the case of Elias 2-27, Pérez et al. 2016; Huang et al. 2018; Paneque-Carreño et al. 2021). However, weighing protoplanetary discs is challenging. Indeed, the main disc gaseous component, molecular hydrogen (H₂), emits faintly in cold environments, thus making its direct detection extremely difficult. Consequently, one must rely on other gas-mass tracers. A first approach consists in measuring fluxes of alternative gaseous tracers such as CO isotopologues (e.g. Miotello et al. 2017), with ¹²CO being the most abundant. Then, by assuming a model-dependent conversion factor from CO isotopologues to H₂, one can estimate the total mass. Another viable tracer is dust: the gas mass can be estimated by measuring the dust mass from millimetre flux assuming the disc is optically thin and setting a value for the gas-to-dust ratio (e.g. 100 as in the ISM; Draine 2003). However, all these methods rely on conversion factors – which lack accurate constraints – and on the assumption that the observed dust is optically thin (for a more comprehensive review, see Miotello et al. 2023). Another recently investigated tracer is hydrogen deuteride (HD; e.g. Bergin et al. 2013; Trapman et al. 2017). When the disc vertical structure is known, a more reliable estimate of the disc mass can be obtained using HD, which tends to yield higher values compared to other tracers (Trapman et al. 2017). A sufficiently large sample of measurements to test the strength of this method is still missing due to the absence of plans for a far-infrared observatory with the capability to observe HD. Recently, Trapman et al. (2022) investigated the possibility of using a combination of N₂H⁺ and C¹⁸O to determine the CO-to-H₂ ratio (Öberg et al. 2023), which nevertheless still depends on the assumed disc structure and N₂ abundance. Lastly, a handful of additional approaches have been proposed in an attempt to provide a better characterisation of this disc property. For example, by looking at multi-wavelength observations and comparing substructures to get an estimation of the Stokes number, and thus of the total disc mass or local surface density (Powell et al. 2019; Veronesi et al. 2019). Also, by investigating the morphology of a kinematic feature linked to GI, known as GI-wiggle (Hall et al. 2020; Longarini et al. 2021), Terry et al. (2022) recently found an approximately linear relationship between the amplitude of this wiggle and the disc-to-star mass ratio, thus obtaining an estimate of the mass of gravitationally unstable discs.

An alternative approach to tracers and dust-to-H₂ conversion has recently been proposed (Veronesi et al. 2021; Lodato et al. 2023), which consists in dynamically measuring the mass of the disc. Indeed, when the disc is massive enough compared to the star, the gravity of the disc contributes considerably to the azimuthal gas velocity, making it super-Keplerian. In particular, the azimuthally averaged rotation curve is expected to transition from an inner Keplerian profile to an outer flatter one, indicative of the disc contribution. If high-resolution kinematic data are available, the observed rotation curve can be fitted, allowing a dynamical mass estimate to be obtained. This method has been widely used in galaxies (e.g. Barbieri et al. 2005) and has also been applied to AGN discs (e.g. Lodato & Bertin 2003). It is also worth noting that it does not require precision on previous stellar mass measurements, because it simultaneously fits M_* and the disc mass, M_d, from the rotation curve. To date, in the context of protoplanetary discs, the method has been applied to three systems, Elias 2-27 (Veronesi et al. 2021), IM Lup, and GM Aur (Lodato et al. 2023), providing a disc-to-star mass ratio (~0.1–0.4) for the first two that is consistent with a scenario where the detected spirals are induced by GI.

This method of dynamical measurement, which is efficient for massive discs, relies on a simple model for the structure and the dynamics of the disc. As such, it is important to determine the range of disc-to-stellar mass within which it can provide sufficiently accurate results. In this study, we provide a benchmark of the method by performing a set of smoothedparticle hydrodynamics (SPH) and radiative transfer simulations (PHANTOM+MCFOST, Price et al. 2018; Pinte et al. 2009), taking into account the gravity of the disc.

This article is organised as follows: in Sec. 2 we describe the physical model use for the gas. In Sec. 3, we present and benchmark our set of numerical simulations. In Sect. 4, we detail the analysis workflow and present our results. In Sec. 5, we present our analysis of the uncertainties. In Sect. 6, we discuss these uncertainties and highlight some possible limitations and future outlooks. We finally draw our conclusions in Sec. 7.

2 Physical model

We consider a smooth axisymmetric circular disc, and we neglect the effect of the magnetic field and dust feedback. Indeed, the former brings only a small correction to the turbulence and magnetic pressure, while the latter may only significantly modify the effect of the pressure gradient of the gas for high dust-to-gas ratios, which is assumed not to be the case for the systems studied here. We also do not consider discs with embedded planets. Low-mass planets would introduce a negligible correction of the order of ∝M_p/M_*. Massive planets would create gaps and modify the gas surface densities, but this is beyond the scope of this first exploration. Within these assumptions, the azimuthal velocity v_ϕ of a gas particle at a cylindrical distance R and height z from a central object of mass M_* is related to the gas pressure P and the gravitational potential of the disc Φ_d according to $$v_{\varphi }^2=\frac{\mathcal{G}{M}_{\star }{R}^2}{{(R^2+z^2)}^{3/2}}+\frac{R}{\rho }\frac{\mathrm{\partial }P(R,z)}{\mathrm{\partial }R}+R\frac{\mathrm{\partial }{\mathrm{\Phi }}_{\mathrm{d}}}{\mathrm{\partial }R}.$$(1)

The three contributions of the right hand-side of Eq. (1) correspond to the gravity of the central star, the radial pressure gradient, and the self-gravity of the disc. We assume that the surface density of the gas follows a self-similar profile set by turbulent viscosity (Lynden-Bell & Pringle 1974): $$\mathrm{\Sigma }(R)=\frac{M_{\mathrm{d}}(2-\gamma )}{2\pi R_{\mathrm{c}}^2}{\left(\frac{R}{R_{\mathrm{c}}}\right)}^{-\gamma }\exp [-{\left(\frac{R}{R_{\mathrm{c}}}\right)}^{2-\gamma }],$$(2)

where R_c is the tapering radius, M_d the disc mass, and γ = 1 is the power-law index. Although this profile is commonly used for non-self-gravitating discs, it can also be chosen in the self-gravitating case since the Keplerian shear is not affected by the disc SG. We also assume the disc to be vertically isothermal with the density profile given by $$\rho (R,z)={\rho }_0(R)\exp [-\frac{R^2}{H^2}(1-\frac{1}{\sqrt{1+z^2/R^2}})]\text{, }$$(3)

and the temperature profile to be a power law with the radius $$T(R)=T_{100}{\left(\frac{R}{100\mathrm{a}\mathrm{u}}\right)}^{-q},$$(4)

where H is the disc thickness, ρ₀ is the density at the mid-plane, and T₁₀₀ is a normalisation factor at 100 au. We stress that the discs considered are not only vertically but also locally isothermal (see Sec. 3.1).

Following, for example, Lodato et al. (2023), the rotation profile is $$\begin{array}{rl}{v}_{\varphi }^2=& v_{\mathrm{K}}^2\{1-[{\gamma }^{\mathrm{\prime }}+(2-\gamma ){\left(\frac{R}{R_{\mathrm{c}}}\right)}^{2-\gamma }]{\left(\frac{H}{R}\right)}^2\\ & -q(1-\frac{1}{\sqrt{1+(z/R{)}^2}})\}+v_d^2,\end{array}$$(5)

where ${\gamma }^{\mathrm{\prime }}=\gamma +3/2+q/2,v_{\mathrm{K}}^2=\mathcal{G}{M}_{\star }/R$ and $$\begin{array}{rl}{v}_{\mathrm{d}}^2=& \mathcal{G}{\int }_0^{\mathrm{\infty }}[K(k)-\frac{1}{4}\left(\frac{k^2}{1-k^2}\right)\\ & \times (\frac{r}{R}-\frac{R}{r}+\frac{z^2}{Rr})E(k)]\sqrt{\frac{r}{R}}k\mathrm{\Sigma }(r)\mathrm{d}r,\end{array}$$(6)

where K(k) and E(k) are complete elliptic integrals and k² = 4Rr/[(R + r)² + z²] (Binney & Tremaine 1987; Bertin & Lodato 1999). Here the boundaries of the integral are chosen to correspond to the inner and outer radius of the disc. The terms have been rearranged in four contributions: (1) a purely circular Keplerian term, $v_{\mathrm{K}}\equiv \sqrt{\mathcal{G}{M}_{\star }/R}$; (2) a term related to the radial pressure gradient in the midplane, which gives a sub-Keplerian contribution; (3) another sub-Keplerian term originating from the pressure gradient and the stellar contribution in the vertical direction; and finally (4) the self-contribution of the disc (see Eq. (6)). To be able to compare some of these contributions, we present some orders of magnitude. The correction introduced by the vertical dependence of the stellar potential is ~(z/R)² and is therefore of the same order as the correction due to the pressure gradient, and therefore should not be neglected (we derived the disc mass of Elias 2-27 with the updated model, obtaining a result that is consistent with Veronesi et al. 2021; see Appendix A). In Appendix B, these equations are provided in a dimensionless form for order-of-magnitude estimates.

Fig. 1 Rotation curves extracted from SPH simulations. Left: rotation curves in the midplane of the disc extracted from the SPH simulations with H/R100au = 0.075 after 20 orbits at the outer disc radius, for different disc-to-star mass ratios (solid lines). The black dash-dotted line shows the rotation curve of the model given by Eq. (5), including the pressure gradient term with zero disc mass. Except for profiles obtained with a disc-to-star mass ratio of 0.01, the rotation curves are distinct from a model without the disc self-gravity. Right: same as left panel but with Mdisc = 0.1 M⊙, with rotation curves extracted for different heights (blue solid lines), and compared to the one obtained in the midplane (z = 0, black dashed-dot line).

3 Numerical simulations

3.1 Hydrodynamical simulations

We performed a suite of 3D SPH simulations of gaseous protoplanetary discs using the code PHANTOM (Price et al. 2018). The system consists of a central star of mass M_* = 1 M_⊙ surrounded by a gas disc with mass M_disc = [0.01, 0.025, 0.05, 0.1, 0.15, 0.2] M_⊙. These simulations can be rescaled in terms of disc-to-star mass ratio (see the derivation in Appendix B). The disc extends from R_in = 10 au to R_out = 300 au, and is simulated with 10⁶ SPH particles. The initial profile of surface density is a relaxed exponential tapered power law of Eq. (2), with R_c = 100 au and γ = 1. Simulations include the disc self-gravity as described in Price et al. (2018), adopting a locally isothermal equation of state $P=c_{\mathrm{s}}^2{\rho }_{\mathrm{g}}$ with q = 0.5 for the power-law index of the temperature radial profile (see Eq. (4)). As a consequence of the absence of cooling, the morphology of the discs will be smooth and axisymmetric, without spiral density waves due to gravitational instability. This is consistent with the fact that, to date, most of the observed discs are axisymmetric (with or without substructures; Ohashi et al. 2023). Effects of self-gravity can be important even if the disc does not host spirals. The choice of the index for radial profile of temperature would remain valid even if we were to take into account dust in the hydrodynamical modelling, because variations in the temperature q-index are not expected to significantly impact the evolution of dust (Pinte & Laibe 2014).

The disc is vertically extended by initially setting up a disc aspect ratio of (H/R)_c = 0.075 with a Gaussian profile for the volume density (see Eq. (3)), ensuring nearly vertical hydrostatic equilibrium. These locally isothermal simulations do not require an extra cooling term. As such, we model the angular momentum transport throughout the disc using the SPH shock to capture viscosity (Price et al. 2018, see Sect. 2.6) with α_AV ≃ 0.19, which results in a Shakura & Sunyaev (1973) viscous parameter of α_ss ≈ 0.005. The parameters used are given in Appendix C.

3.2 Preliminary test: Hydrodynamical rotation curves

As a first benchmark, we first extracted the rotation curves directly from the hydrodynamical simulations after having relaxed the disc to centrifugal equilibrium (t ~ 20Ω_out, where Ω_out is the orbital period of the disc outer radius). The rotation curves are obtained by interpolating azimuthal velocities of the SPH particles along the desired coordinates (R, z), and then azimuthally averaging the results. This step provides an uncertainty for the rotational velocity that is estimated as a standard deviation. Figure 1 shows the rotation curves obtained by this procedure. The left panel shows the rotation profiles obtained in the midplane (R, z = 0), with disc masses ranging from 0.01 to 0.2 (blue lines). As a reference, the black line provides the rotation curve of the model given by Eq. (5)), including the pressure gradient term, but setting the mass of the disc to zero. As expected, the contribution of the disc self-gravity increases with the disc-to-star mass ratio. The right panel shows the rotation curves of the simulation with M_disc = 0.1 at different heights ranging from z = 0 to z = 50 au. As expected from the model, the difference in azimuthal velocity between different heights decreases with radius.

To benchmark the validity of the formulation given by Eq. (5), we fit the extracted rotation curves with the code dysc¹, which implements the model described in Sec. 2. The results are shown in Table 1. The fitting procedure relies on a Markov Chain Monte Carlo using the Python package EMCEE (Foreman-Mackey et al. 2013). We used a Gaussian likelihood with flat priors. The code simultaneously fits both the masses of the star and the disc, as well as the truncation radius R_c, similarly to the procedure described in Lodato et al. (2023). After this step, we recovered the mass of the disc with an uncertainty of ~5–12% (see Table 1) for disc-to-star mass ratios ranging between 0.025 and 0.2. In these cases, discs are sufficiently massive for the self-gravity to have a measurable contribution. For a disc-to-star mass ratio of 0.01, the disc mass cannot be extracted accurately even from the numerical simulations.

3.3 Radiative transfer and synthetic observations

We then computed the thermal structure and synthetic observations of our numerical discs by performing 3D radiative transfer simulations with the code MCFOST (Pinte et al. 2006, 2009). We used a Voronoi tessellation, where each MCFOST cell corresponds to a gas SPH particle. The goal is to generate mock ¹²CO and ¹³CO isotopologue channel maps, from which we can extract disc rotation curves.

The main inputs for the radiative transfer modelling are the source of luminosity, the gas density profile extracted from the simulations, and models for dust opacities and densities. Self-gravitating isothermal discs do not present any kind of substructure and the underlying dust density profile can therefore be assumed to be smooth. We do not account for an eventual dust drift, since Stokes numbers are expected to be small in these discs, and we consider sufficiently short disc evolutions. The dust contribution to the gravitational potential of the disc is also negligible. As such, we adopt a constant dust-to-gas ratio of ~10⁻², which corresponds to a standard averaged ISM (Draine 2003).

The thermal structure of the disc is computed with the following assumptions. At first, emission is at local thermal equilibrium and T_gas = T_dust. Dust is treated as a mixture of silicates (70%) and carbon (30%) (Draine & Lee 1984), and the optical properties are calculated using Mie theory for spheres (Andrews et al. 2009). Opacities are computed following the procedure described in the DIANA model (Woitke et al. 2016; Min et al. 2016). We assume an ISM-like grain size distribution (dn/da ∝ a^−3.5), with a_min = 0.01 μm and a_max = 1 mm. The disc is heated passively, that is, the source of radiation is only the central star, which is assumed to radiate isotropically with a Kurucz spectrum at 4000 K. The expected channel maps are computed via ray tracing, using 10⁷ photon packets to sample the radiation field. The disc inclination with respect to the line of sight is i = 30°, and the system is simulated to be located at 140 pc, which corresponds to a typical protostellar disc in a star-forming region such as Taurus. For ¹²CO and ¹³CO, we consider the J = 2–1 transition and assume abundances of 10⁻⁴ and 1.4 × 10⁻⁶ respectively. MCFOST simulations are post-processed with PYMCFOST² by simulating a velocity resolution of 0.1 km s⁻¹. We then spatially convolve the channels with a Gaussian beam of 0.1 arcsec, similarly to the value of the MAPS survey (Öberg et al. 2021). We finally introduce Gaussian noise with an RMS of 5mJy beam⁻¹.

Figure 2 shows an example of the channel maps obtained for a simulation with a disc mass of 0.1 M_⊙. The ¹²CO channels appear more radially extended compared to the ¹³CO channels. The ¹³CO channels are also less patchy, because the τ = 1 region is located closer to the midplane, a region with higher resolution. As a result, the ¹³CO disc rotation curves will be less noisy than the curves extracted from the ¹²CO channels (see Fig. 5). We note that while this is the case for simulations, in real observations we do not expect ¹³CO rotation curves to be less noisy than the ¹²CO. Indeed, ¹³CO is much less abundant than ¹²CO, and at the same sensitivity, the signal-to-noise ratio (S/N) in ¹³CO should be lower than that for ¹²CO.

4 Method and results

4.1 Workflow description

From the cubes of synthetic data, we first determined the height of the emitting layer with DISKSURF (Teague et al. 2021). The emitting layer is defined as the region from which the emission we observe originates. Isotopologues have specific optical depth and column density, and as such, the emitting layer can be located farther from or closer to the midplane of the disc. Having a precise estimate of the height of the emitting layer is crucial in order to correctly deproject azimuthal velocities and compute the model of Eq. (5) at the correct height.

We then used EDDY (Teague 2019) to extract the rotation curves. Two methods have been developed to retrieve the observed velocity (specifically to fit the line centroids): one quadratic and the other Gaussian. The aim of the former is to fit a quadratic curve to the brightest pixel on either side of the maximum of the pixel. The fit depends on the velocity sampling, and is sensitive to the channel correlation. On the other hand, with the latter method, which finds the line centre by fitting a Gaussian profile to it, the selected velocity range may affect the result in cases with skewed profiles.

To perform our analysis, we chose the quadratic method; we briefly discuss this choice in Sec. 4.3 (for an extended discussion of the methods, see e.g. Lodato et al. 2023 and Teague & Foreman-Mackey 2018). Finally, we fit the rotation curves with the model of Eq. (5) using the code DYSC (see Sec. 2).

4.2 Retrieving the height of the emitting layer

We first used the code DISKSURF (Teague et al. 2021) to extract the height of the emitting layer zⁿ (r_i) of the ith radial bin for the nth channel. DISKSURF implements the method presented in Pinte et al. (2018). Briefly, this is a geometrical method that allows the user to trace the emission maxima of each velocity channel and reconstruct the position of the CO layers. For each bin, we determined the heights $z_{16}^n\left(r_i\right),z_{50}^n\left(r_i\right)\text{, and }{z}_{84}^n\left(r_i\right)$, which correspond to the 16th, 50th, and 84th percentiles of the distribution, respectively. We then used the tapered power law, $$z(R)=z_0[\mathrm{arcsec}]{\left(\frac{R}{1\mathrm{arcsec}}\right)}^{\psi }\exp [-{\left(\frac{R}{R_t}\right)}^{q_t}],$$(7)

to continuously parametrise the emitting surfaces z₁₆ (r), z_em (r) = z₅₀ (r), and z₈₄ (r). The parameters obtained were used to retrieve the height of the molecular emission surface z_em (r) as well as to estimate the uncertainties associated with the procedure.

Table 2 presents the results obtained for z_em(r), which we subsequently used as a reference surface in the following analysis. As an example, Fig. 3 shows the emission surface (grey markers) obtained with DISKSURF for the ¹²CO datacube from a mock observation generated by a simulation with M_disc = 0.1 M_⊙. The three lines are the fits of the data points, namely z_em = z₅₀ (solid line), z₁₆ (dotted line), and z₈₄ (dashed line). The scatter between the three fitted surfaces is important with this procedure: it is therefore crucial to refer to this dispersion when analysing uncertainties on the rotation curve.

As expected, we also find that there is a clear trend between the height of the emitting layer and the disc mass: more massive discs have a higher emitting layer, because their surface density is higher and they are optically thicker. Also, the height of the ¹³CO is lower with respect to the ¹²CO because of its reduced abundance, tracing a region closer to the midplane. This is shown in more detail in Figs. D.1 and D.2.

Finally, for consistency, we decided not to include the simulation with M_d/M_* = 0.2 in our analysis at this stage. Indeed, in this case, the ¹²CO emission is poorly reproduced because the SPH resolution at the emitting surface is not sufficient (see Sec. 6.1 for a resolution analysis).

Fig. 2 Example of channel maps obtained with MCFOST for the 12CO (top row) and 13CO (bottom row) isotopologues, for a simulation with disc mass of Md = 0.1 M⊙. Different velocity channels are displayed, from 0.25 to 0.75 km s−1 (going from left to right). The chosen Gaussian beam used to convolve the image matches the value of the MAPS survey (0.1″), and is shown with a grey circle in the left bottom corner of each channel map. The velocity resolution is chosen to be 0.1 km s−1 (as in the MAPS survey).

Fig. 3 Example (12CO, Mdisc = 0.1 M⊙) of one of the emitting layers derived with DISKSURF (grey dots) and then fitted with an exponentially tapered power law (solid line). We also show the 16th (dotted line) and the 84th (dashed line) percentiles, which we take into account to compute the rotation curve errors.

Fig. 4 Differences between rotation curves obtained with the Gaussian and quadratic methods. Left panel: comparison between the rotation curves from the model (black line, Eq. (5) with a disc mass of Md = 0.1 M⊙), the Gaussian method (blue line), and the quadratic method (orange line). We observe that the Gaussian curve is systematically shifted with respect to the model. Right panel: difference between the extracted rotation curve (Gaussian method in blue, quadratic method in orange) and the model vϕ from Eq. (5).

4.3 Extracting the rotation curve

We extract the azimuthal velocity with EDDY (Teague 2019) using the harmonic oscillator method, following Teague et al. (2018a,b). Further details of the extraction methods can be found in Sect. 3.3.1 of Lodato et al. (2023).

For each emitting layer, that is, z₁₆, z_em = z₅₀, and z₈₄, we computed three rotation curves, v₁₆, v₅₀, and v₈₄. We then assumed the azimuthal velocity of the system to be v;_em = v₅₀, with an uncertainty σ_v estimated to be $${\sigma }_v=\sqrt{{|v_{84}-v_{16}|}^2+{\sigma }_{\text{eddy }}^2}\text{, }$$(8)

where σ_eddy is the numerical error obtained with EDDY. This procedure refines the approach of Lodato et al. (2023) by including uncertainties related to the estimate of the height of the emitting layer.

As we fixed the inclination of the disc to i = 30°, we used the quadratic method to extract the rotation curve; this method fits only the peak of emission, which is less sensitive to the lower surface. Indeed, using the Gaussian method instead provides the following bias: the emission coming from the lower surface systematically shifts the position of the line centroids, resulting in a systematic error for the velocity estimate. Figure 4 shows a comparison between the quadratic and the Gaussian methods. Interestingly, despite being smoother, the curve obtained with the Gaussian method underestimates the velocity in the inner disc and overestimates it in the outer disc. This happens because the method tries to fit a double-peaked spectrum with a single Gaussian. We note that in this work we tested both methods, and despite both being biased, we chose the quadratic method (although noisier) because the Gaussian one presents an overall larger systematic shift of the velocity caused by the lower surface. Moreover, with the quadratic method we already manage to recover the correct disc and star masses. For completeness, we also tested our workflow using a rotation curve retrieved with a double Gaussian fit. However, the double-Gaussian method fails to retrieve the correct velocity in the outer region of the disc (R > 300 au), where the S/N is too low (it fits the noise but not the signal; see Appendix E).

Figure 5 shows the rotation curves obtained with the quadratic method for both CO isotopologues following the procedure described above (¹²CO in the left panel, ¹³CO in the right panel). We compare these with the analytical rotation curve obtained assuming a disc of zero mass (black dash-dotted line). As expected, the rotation curve derived with this approximation fits only the disc with the lower mass (M_d = 0.01 M_⊙), because this is the only case where the disc contribution can be considered as negligible. We show how we confirmed this in Sec. 4.5.

Fig. 5 Quadratic rotation curves obtained for all the disc mass simulations (from Mdisc = 0.01 M⊙ to Mdisc = 0.2 M⊙) for the 12CO isotopologue (solid lines in the left panel) and the 13CO isotopologue (dashed lines in the right panel). The black dash-dotted line shows the analytical rotation curve obtained by only considering the star and pressure gradient contribution.

4.4 Fitting procedure

For every simulation, we obtain two rotation curves: one for ¹²CO and another for ¹³CO. These curves are then fitted using the self-gravitating model of Eq. (5). The free parameters are the star and disc mass, and the disc truncation radius. The fits were performed with an MCMC algorithm as implemented in EMCEE (Foreman-Mackey et al. 2013). We chose a simple Gaussian likelihood, and flat priors of [0,5] M_⊙ for the star mass³, [50, 300] au for the truncation radius, and [0, 0.5] M_⊙ for the disc mass. We chose 250 walkers and 1000 steps (having verified that convergence is reached). We fit the two isotopologues both individually and then simultaneously.

4.5 Fit results

Table 3 shows the values retrieved for the star mass M_*, the disc mass M_d, and the disc truncation radius R_c for simulations with disc masses ranging within M_d = 0.01–0.2 M_⊙, and a disc aspect ratio of 0.075. The disc-to-star mass ratio threshold below which the expected value cannot be recovered is 0.05. For reasons related to resolution, we also exclude the simulation with M_d = 0.2 M_⊙ (see Sec. 4.2). From our procedure, we are able to measure a non-zero value for the disc mass (M_disc ≥ 0.05 M_⊙), which indicates that the model for the rotation curve should include the term from the gravity of the disc.

5 Uncertainties

5.1 Physical parameters

In the treatment of real data, the three principal sources of uncertainty are the height of the emitting layer, the aspect ratio of the disc, and its inclination. Uncertainties on z(r) are estimated with the values of z₁₆ and z₈₄ found previously. To estimate uncertainties associated with H/r and i, we first note that in the fitting procedure of Sec. 4, we enforce the values of H/r and i that were set in the numerical simulation. We therefore estimated uncertainties associated to those parameters by carrying the fit over the same synthetic channel maps while varying H/r and i over a range of values. For H/r, we simply performed new fits over the previously extracted rotation curves. For a setup value of H/R = 0.075, we tested H/R = [0.05, 0.07, 0.08, 0.1]. For inclinations, uncertainty estimates require the extraction of new disc emitting layers, necessitating that we obtain novel rotation curves. For a setup value of i = 30°, we tested i = [27, 29, 31, 33]°. The values of the uncertainties obtained by this procedure are summarised in Table 3. Figure 6 shows these uncertainties for discs whose mass can be accurately recovered by the fitting procedure. This excludes the md0.01 and md0.025 (small masses), and md0.2 (lack of resolution). A key result of this study is that disc masses of self-gravitating discs can be estimated from channel maps with averaged systematic uncertainties of the order of 25%. The three parameters H/r, i, and z make similar contributions as sources of error. No clear trend appears when varying H/r and z. Values that differ from the expected one still provide mass estimates with the same level of uncertainty. On the other hand, precise values of i provide uncertainties of the order of 5–10%, whereas an error of a few degrees provides uncertainties of the order of 20–30%. Hence, our procedure is relatively reliable in its capacity to accurately determine disc masses. For instance, when recovering a disc mass of 0.1 M_⊙, the estimated range spans from 0.075 to 0.125. We note that one obtains simultaneously ΔM_*/M_* ~ 8% for an inclination error of ±3° and ΔR_c/R_c ~ 15% for a reasonable disc aspect ratio range of ±0.025 (see Appendix F for a discussion).

5.2 Spatial resolutions

In order to test how the minimum detectable mass is affected by the spatial resolution, we produced new synthetic images by convolving the channel maps previously obtained in Sec. 3.3 with different Gaussian beam sizes [0.2, 0.3]″. We then extracted new rotation curves with the procedure described in Sects. 4.2 and 4.3, and we fit them with the DYSC code (Sec. 4.4). In Fig. 7, we report the results obtained for the disc mass (different markers and lines), showing uncertainties obtained as a function of the spatial resolution of the observations. When M_d = 0.05 M_⊙, the uncertainty ranges from ~50 to 100%, while for the other two cases, M_d = [0.1–0.15] M_⊙, it ranges within ~30–40%ο, which is larger compared to the 0.1″ case discussed above (~25%; see Sec. 5.1). This analysis shows that a resolution of 0.1″ is recommended in order to obtain a reliable measurement.

Fig. 6 Uncertainties related to the fitting procedure for the disc mass Md. Estimates are given as functions of the aspect ratio of the disc H/R (left column), the disc inclination (middle column), and the disc emitting layer z(R) (right column). Different markers and line styles represent results obtained for simulations with different disc masses.

6 Discussion

6.1 Numerical resolution

As the height of the emitting layer increases with the disc mass, in the md0.2 simulation, we notice that the emitting layer is higher than the typical height of SPH particles, while this is not true for the case with low disc mass (for an example, see Fig. D.3). Indeed, the number density of SPH particles decreases with z, because most of the disc mass is concentrated in the midplane. This implies that the size of the Voronoi grid cells used in MCFOST is increasing with z, which results in artificially low-resolution ¹²CO data cubes for the high-disc-mass simulation. To test this, we performed additional hydrodynamical simulations with 2, 4, and 6 × 10⁶ SPH particles. Figure 8 shows that by increasing the number of SPH particles, the uncertainty we obtain for the disc mass decreases until we reach the ‘true’ disc mass value for the 6M particles simulation. The blue dashed line represents the semi-log linear fit of the data points.

Fig. 7 Uncertainties obtained from our fitting procedure for the disc mass ΔMd/Md and for the combined fit procedure with different spatial resolutions [0.1,0.2,0.3]″. Different markers (and lines) represent results obtained for simulations with different disc masses.

Fig. 8 Uncertainties obtained from our fitting procedure for the disc mass ΔMd/Md for the combined fit procedure with different numerical resolutions: [1, 2, 4, 6] × 106 SPH particles. The blue dashed line is the semi-log linear fit of the data points with slope m = −0.58. As expected, accuracy increases with numerical resolution.

6.2 Limitations and further work

The main assumptions adopted in this study are as follows:

• Model: axisymmetric discs, no magnetic field, no planets;

• Hydrodynamic simulations: gas only, as the dust-to-gas ratio is assumed to be so low that it does not affect the result. The disc is assumed to be locally and vertically isothermal;

• Radiative transfer: T_gas = T_dust (LTE), Mie theory for optical properties, ISM-like grain size distribution, dust-to-gas ratio of ~0.01.

• Post-processing of synthetic images: the simulated ALMA images are obtained by convolving with a Gaussian beam and the noise has been added directly to the resulting images. Thus, the filtering of an interferometer and the noise properties of ALMA observations are not fully captured with these simple assumptions. As a test, we also applied the method presented in this work to a simulated ALMA image obtained with CASA, finding no difference in the fitted disc parameters.

Our study is restricted to stable self-gravitating discs, excluding those with peculiar substructures such as GI spirals. It would be important to reproduce the benchmark performed in this study with more general cases; for example, to confirm masses estimated for Elias 2-27, IM Lup, and GM Aur. This study highlights the need to investigate how the extraction of the rotation curve is affected by the presence of non-axisymmetric structures (e.g. by including cooling and allowing the disc to develop GI in the shape of spiral structures). Indeed, non-axisymmetric structures might be important in young objects that are still accreting from the molecular clouds (e.g. streamers). It may also be worth investigating whether or not the presence of an embedded planet or a central binary affects disc mass estimates. Finally, the case of non-isothermal discs is being investigated by Martire et al. (2024).

7 Conclusions

In this study, we benchmarked the determination of disc masses using rotation curves extracted from channel maps as introduced in Veronesi et al. (2021) and Lodato et al. (2023). We generated controlled mock observations from PHANTOM SPH simulations with different disc masses that are post-processed with the MCFOST radiative transfer code. We adopted spatial and spectral resolutions that are the same as those of the MAPS survey (0.1″ and 0.1 km s⁻¹).

We find that for M_d/M_* ≳ 0.05, we recover the correct disc-to-star mass ratio with a typical accuracy of ~25%, with equal uncertainty contributions from determinations of the aspect ratio and the inclination of the disc, as well as the height of the emitting layer (the truncation radius and the stellar mass are also obtained with uncertainties of ~15% and ~8%, respectively). For lower resolutions (0.2–0.3″), disc masses recovered shift towards larger values (0.1 M_⊙ for 0.2″). This constraint should be accounted for when proposing ALMA kinematics observations in order to dynamically measure the disc mass. This method can be used to calibrate existing methods such as those based on flux.

Acknowledgements

The authors want to thank the referee for constructive comments that improved this manuscript. This work has received funding from the ERC CoG project PODCAST No 864965, and from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 823823 (RISE DUSTBUSTERS project). S.F. is funded by the European Union (ERC, UNVEIL, 101076613). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. We used the EMCEE algorithm Foreman-Mackey et al. (2013), the CORNER package Foreman-Mackey (2016) to produce corner plots, and PYTHON-based MATPLOTLIB package Hunter (2007) for all the other figures. Computational resources were provided by the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon, and the INDACO Platform, a project of High Performance Computing at the Università degli Studi di Milano. The authors thank Pietro Curone, Andrés F. Izquierdo, Sean Andrews, Richard Teague, Giovanni Rosotti and Enrico Ragusa for useful discussions.

Appendix A Elias 2-27 updated fit results

We derived with the updated model the disc mass of Elias 2-27, obtaining a result that is consistent with the findings of Veronesi et al. (2021). In particular, when fitting the west and the east side of the disc simultaneously, and the two isotopologues (¹³CO and C¹⁸O), we obtain a star mass of ~0.40 M_⊙ and a disc mass of −0.13 M_⊙ (~0.46 M_⊙ and ~0.08 M_⊙ in Veronesi et al. 2021).

Appendix B Dimensionless model

Equation 5 can be expressed in terms of x = R/R_c and ζ = z/R_c as $$v_{\varphi }^2/\frac{G{M}_{\star }}{R_c}=f_1(x,\zeta )+\frac{1}{2\pi }\frac{M_{\mathrm{d}}}{M_{\star }}{f}_2(x,\zeta ),$$(B.1)

where $$f_1(x,\zeta )=\frac{1}{x}\{1-[{\gamma }^{\mathrm{\prime }}+(2-\gamma )x^{2-\gamma }]{\left(\frac{H_{\mathrm{c}}}{R_{\mathrm{c}}}\right)}^2{x}^{1-q}-q(1-\frac{1}{\sqrt{1+{\zeta }^2}})\}\text{, }$$

and $$\begin{array}{r}{f}_2(x,\zeta )={\int }_0^{\mathrm{\infty }}[K(y)-\frac{1}{4}\left(\frac{y^2}{1-y^2}\right)(\frac{u}{x}-\frac{x}{u}+\frac{{\zeta }^2}{xu})E(y)]\\ \sqrt{\frac{u}{x}}y\tilde{\mathrm{\Sigma }}(u)\mathrm{d}u,\end{array}$$

where H_c denotes the scale height at the tapered radius $R_{\mathrm{c}},y\equiv \frac{4u}{(u+x{)}^2+{\zeta }^2},\tilde{\mathrm{\Sigma }}\left(r/R_{\mathrm{c}}\right)\equiv \mathrm{\Sigma }(r)/\frac{{\tilde{M}}_{\mathrm{d}}}{2\pi R_{\mathrm{c}}^2}$ denotes the dimensionless surface density, and ${\tilde{M}}_{\mathrm{d}}\equiv 2\pi {\int }_{R_{\text{in }}}^{R_{\text{out }}}R\mathrm{\Sigma }(R)\mathrm{d}R\simeq M_{\mathrm{d}}$ denotes the true mass of the disc. Hence, the azimuthal velocity can be rescaled in terms of the disc-to-star-mass ratio M_d/M_*.

Appendix C Simulation parameters

The parameters used in the hydrodynamic simulations described in Sect. 3.1 are given in Table B.1.

Appendix D Properties of the emitting layer

Figure D.1 shows the height of emitting layers as a function of the radius for different disc masses. As expected, the emitting layer is higher for large disc masses, because the column density is larger. Fig. D.2 highlights the difference between the height of layers derived for the ¹²CO (orange line) and ¹³CO (blue line) isotopologues. ¹³CO is optically thinner than ¹²CO and therefore traces a region closer to the midplane; its emitting layer is also lower compared to that of ¹²CO.

Fig. D.1 Emitting layer derived with DISKSURF for different disc masses for the 12CO isotopologue.

Fig. D.2 Comparison between the emitting layer of 12CO (solid line) and 13CO (dashed line) for a disc mass of 0.1 M⊙.

Fig. D.3 shows a comparison between the emitting layers derived for 0.2 M_⊙ (top panel) and 0.025 M_⊙ (bottom panel). Emitting layers are plotted with orange solid lines, while the SPH particles particles are represented with cyan markers. The number density of SPH particles decreases with z, as most of the disc mass tends to be concentrated in the midplane. This implies that the size of the Voronoi grid cells increases with z, resulting in an artificially low resolution in the upper layer of the ¹²CO data cubes. This becomes problematic in the case of the high-disc-mass simulations (e.g. here 0.2 M_⊙ for 10⁶ SPH particles), as the disc is optically thicker and the emitting layer is higher with respect to the case of lighter discs.

Fig. D.3 Distribution of SPH particles (cyan dots) and the emitting layer obtained with DISKSURF (orange line) for two simulations, one with high mass (upper panel) and one with low mass (lower panel).

Appendix E Double Gaussian method

Fig. E.1 shows a comparison between the analytical rotation curve obtained with the self-gravitating model (Eq. 5), and the velocity values obtained with the quadratic and the double Gaussian methods, computed at two different radii. We note that, in the outer region, the value extracted with the double Gaussian method deviates further from the expected value. Moreover, in the outer region the S/N is too low, and the double Gaussian method fits the noise rather than the signal.

Fig. E.1 Comparison among the analytical rotation curve obtained with the self-graviting model (Eq. 5, black solid line), and the velocity obtained with the quadratic (red circle) and the double Gaussian method (blue star), computed at two different radii (70 au and 560 au).

Appendix F Other uncertainties

Fig. F.1 shows the uncertainties for the mass of the star (top row) and the disc truncation radius (bottom row) as function of the aspect ratio of the disc H/R (left), the disc inclination (centre), and the disc emitting layer z(R) (right). While for the uncertainties for the disc mass (see Fig. 6), the three parameters H/R, i, and z(r) make similar contributions as sources of error, now we can identify two main contributions. For the mass of the star and the disc truncation radius, the principal sources of uncertainty are the inclination of the disc (panel e in Fig. F.1) and its aspect ratio (panel g in Fig. F.1), respectively. Indeed, varying the inclination with respect to the expected value implies that the annulus along which the observed velocity is computed in order to extract the rotation curve is not at a constant radius. Another source of error is related to the deprojection of the rotation curve. This strongly affects the fitted value of the mass of the star. As for the aspect ratio, it mostly influences the truncation radius, as it scales with (H/R)² (Eq. 5).

Fig. F.1 Uncertainties related to the fitting procedure for the mass of the star M* (top row) and the disc truncation radius Rc (bottom row). Estimates are given as functions of the aspect ratio of the disc H/R (left column), the disc inclination (middle column), and the disc emitting layer z(R) (right column). Different markers and lines styles represent results obtained for simulations with different disc masses.

References

All Tables

All Figures

Fig. 1 Rotation curves extracted from SPH simulations. Left: rotation curves in the midplane of the disc extracted from the SPH simulations with H/R100au = 0.075 after 20 orbits at the outer disc radius, for different disc-to-star mass ratios (solid lines). The black dash-dotted line shows the rotation curve of the model given by Eq. (5), including the pressure gradient term with zero disc mass. Except for profiles obtained with a disc-to-star mass ratio of 0.01, the rotation curves are distinct from a model without the disc self-gravity. Right: same as left panel but with Mdisc = 0.1 M⊙, with rotation curves extracted for different heights (blue solid lines), and compared to the one obtained in the midplane (z = 0, black dashed-dot line).

In the text

Fig. 2 Example of channel maps obtained with MCFOST for the 12CO (top row) and 13CO (bottom row) isotopologues, for a simulation with disc mass of Md = 0.1 M⊙. Different velocity channels are displayed, from 0.25 to 0.75 km s−1 (going from left to right). The chosen Gaussian beam used to convolve the image matches the value of the MAPS survey (0.1″), and is shown with a grey circle in the left bottom corner of each channel map. The velocity resolution is chosen to be 0.1 km s−1 (as in the MAPS survey).

In the text

	Fig. 3 Example (12CO, Mdisc = 0.1 M⊙) of one of the emitting layers derived with DISKSURF (grey dots) and then fitted with an exponentially tapered power law (solid line). We also show the 16th (dotted line) and the 84th (dashed line) percentiles, which we take into account to compute the rotation curve errors.
In the text

Fig. 4 Differences between rotation curves obtained with the Gaussian and quadratic methods. Left panel: comparison between the rotation curves from the model (black line, Eq. (5) with a disc mass of Md = 0.1 M⊙), the Gaussian method (blue line), and the quadratic method (orange line). We observe that the Gaussian curve is systematically shifted with respect to the model. Right panel: difference between the extracted rotation curve (Gaussian method in blue, quadratic method in orange) and the model vϕ from Eq. (5).

In the text

	Fig. 5 Quadratic rotation curves obtained for all the disc mass simulations (from Mdisc = 0.01 M⊙ to Mdisc = 0.2 M⊙) for the 12CO isotopologue (solid lines in the left panel) and the 13CO isotopologue (dashed lines in the right panel). The black dash-dotted line shows the analytical rotation curve obtained by only considering the star and pressure gradient contribution.
In the text

	Fig. 6 Uncertainties related to the fitting procedure for the disc mass Md. Estimates are given as functions of the aspect ratio of the disc H/R (left column), the disc inclination (middle column), and the disc emitting layer z(R) (right column). Different markers and line styles represent results obtained for simulations with different disc masses.
In the text

	Fig. 7 Uncertainties obtained from our fitting procedure for the disc mass ΔMd/Md and for the combined fit procedure with different spatial resolutions [0.1,0.2,0.3]″. Different markers (and lines) represent results obtained for simulations with different disc masses.
In the text

	Fig. 8 Uncertainties obtained from our fitting procedure for the disc mass ΔMd/Md for the combined fit procedure with different numerical resolutions: [1, 2, 4, 6] × 106 SPH particles. The blue dashed line is the semi-log linear fit of the data points with slope m = −0.58. As expected, accuracy increases with numerical resolution.
In the text

	Fig. D.1 Emitting layer derived with DISKSURF for different disc masses for the 12CO isotopologue.
In the text

	Fig. D.2 Comparison between the emitting layer of 12CO (solid line) and 13CO (dashed line) for a disc mass of 0.1 M⊙.
In the text

	Fig. D.3 Distribution of SPH particles (cyan dots) and the emitting layer obtained with DISKSURF (orange line) for two simulations, one with high mass (upper panel) and one with low mass (lower panel).
In the text

	Fig. E.1 Comparison among the analytical rotation curve obtained with the self-graviting model (Eq. 5, black solid line), and the velocity obtained with the quadratic (red circle) and the double Gaussian method (blue star), computed at two different radii (70 au and 560 au).
In the text

Fig. F.1 Uncertainties related to the fitting procedure for the mass of the star M* (top row) and the disc truncation radius Rc (bottom row). Estimates are given as functions of the aspect ratio of the disc H/R (left column), the disc inclination (middle column), and the disc emitting layer z(R) (right column). Different markers and lines styles represent results obtained for simulations with different disc masses.

In the text