© The Authors 2025

1 Introduction

In recent years, high spatial and spectral resolution gas observations from the Atacama Large Millimeter/submillimeter Array (ALMA) have provided detailed insights into the kinematics and velocity fields of planet-forming disks (Pinte et al. 2023). Rotation curves (i.e., the azimuthally averaged radial profiles of the gas rotation velocity) are a direct probe of gas kinematics in protoplanetary disks, and they are crucial to improving our knowledge of several disk properties. First, the gas rotation velocity at different heights is determined by the balance between the stellar gravity, the disk self-gravity, and the pressure gradients. As a consequence, accurate rotation curves allow us to dynamically measure the fundamental properties of young star-disk systems, as stellar mass M_★, disk mass M_d, and disk scale radius R_c. Dynamical masses of protostars were previously estimated by Guilloteau et al. (2014); Simon et al. (2017); Yen et al. (2018); Braun et al. (2021); dynamical measurements of M_d were instead provided through rotation curves fitting by Veronesi et al. (2021) for Elias 2-27, by Lodato et al. (2023) for IM Lup and GM Aur, by Martire et al. (2024) for the whole MAPS sample, and by Andrews et al. (2024) for MWC 480. In addition, this methodology was benchmarked and further investigated with models and simulations in recent works by Veronesi et al. (2024) and Andrews et al. (2024). These quantities, especially M_d and R_c, are not trivial to measure; rotation curves provide complementary accurate estimates via a methodology that is completely independent of line fluxes or continuum flux densities (see Miotello et al. 2023, for a review). Moreover, disks are known to host both radial and vertical gradients in their temperature and density structure (see, e.g., Öberg et al. 2023, for a recent review), which jointly determine the pressure structure. With sets of molecular lines originating from different vertical layers, due to a range of abundances, optical depths, and excitation conditions (e.g., van Zadelhoff et al. 2001; Aikawa et al. 2002; Dartois et al. 2003; Semenov et al. 2008; Henning & Semenov 2013; Dutrey et al. 2014), the rotation curves of these same transitions can be leveraged to extract precious information on the pressure structure of disks at different heights.

A first simple approach that we consider in this paper is to verify that the vertical and radial (2D) temperature structure estimated by other means is directly imprinted in the rotation curves of an ensemble of molecular lines spanning a wide range of vertical layers. We focus on temperature rather than density since the former is more directly determined observationally. Historically, observations of both spectral energy distributions and spatially unresolved gas lines in the far-infrared (FIR) and millimeter regimes have been successfully forward-modeled by radiative transfer and sophisticated thermo-chemical codes to infer the thermal structure of the outer regions (>5–10 au) of protoplanetary disks (see reviews by Miotello et al. 2023; Öberg et al. 2023, and references therein). In recent years, the unique capabilities of ALMA have provided direct stringent constraints on the 2D temperature structure of disks. In particular, with high angular and spectral resolution observations, molecular emitting layers can be reconstructed using geometric arguments from channel maps (Pinte et al. 2018a). In the assumption of optically thick and geometrically thin emission, thermalized lines can thus be used to measure the temperature along their emission layer (e.g., Pinte et al. 2018a; Law et al. 2021b; Paneque-Carreño et al. 2023). Using a set of lines originating from different layers, it is thus possible to fit for parameterized 2D (T(r, z)) temperature distributions (Law et al. 2021b, 2023, 2024; Paneque-Carreño et al. 2023).

Recently, the effect of vertical temperature gradients in rotation curves has been analyzed by Martire et al. (2024). The authors showed that differences in the rotation velocities of ¹²CO and ¹³CO are not explainable through classical vertically isothermal disk models, but instead must be traced back to vertical stratification in temperature.

For this work we took a step forward in the analysis of gas kinematics, following a multi-molecule approach. We simultaneously analyzed the rotation velocities of seven molecular tracers emitting from a wide range of disk layers. In this way, our aim is to provide an accurate estimate of the disk mass and scale radius, based on a multi-molecule fit of the curves with a thermally stratified model. From this fit we also investigated the impact of thermal structure on gas kinematics, testing whether vertical stratification can be detected in the kinematic structure of the disk, probed through multiple tracers located at different heights above the midplane. We focused on the very bright disk of HD 163296 and selected exquisite quality multi-molecule data from the Molecules with ALMA at Planet-forming Scales (MAPS) Large Program (Öberg et al. 2021).

This paper is organized as follows. In Sect. 2, we illustrate in detail the selected data and targeted molecular lines. In Sect. 3 we explain the methodology we followed for our analysis. First, we discuss the extracted molecular optical depth profiles that we used as a complementary test to verify the location of the emitting layers. Then, we present the adopted routine for the rotation curves extraction; we describe the new multi-Gaussian methods we implemented to fit the hyperfine structured spectra of C₂H and HCN, and the corresponding rotation curves. We also highlight a possible observational bias in the rotation curves, arising from beam smearing in the presence of emission gaps. Finally, we describe the model including thermal stratification that we used for the multi-molecule fit of the data. In Sect. 4, we show our results from the multi-molecule fit of multiple rotation curves we performed with the thermally stratified model, and we discuss its impact on gas kinematics. From the fit, we derive dynamical estimates of stellar mass M_★, disk mass M_d, and disk scale radius R_c , and we map the 2D pressure structure of the disk. We also discuss the dependence of our results on the assumed emitting surfaces. We finally draw the main conclusions of this work in Sect. 5.

2 Observations

2.1 The source

In this work, we focus our analysis on HD 163296: thanks to its relative proximity and massive disk, this is one of the most studied Herbig Ae star-disk systems. This protoplanetary disk presents multiple features pointing at ongoing planet formation, such as dust rings, azimuthal asymmetries, deviations from Keplerian velocities, kinks in CO emission, and meridional flows (Isella et al. 2016; Andrews et al. 2018; Isella et al. 2018; Teague et al. 2019; Pinte et al. 2018b; Izquierdo et al. 2022). Observations in scattered light at optical and infrared wavelengths also show radial substructures (Rich et al. 2019). HD 163296 has been observed in CO and other molecular tracers at millimeter wavelengths in several single-dish and interferometric projects (Thi et al. 2004; de Gregorio-Monsalvo et al. 2013; Pegues et al. 2020; Öberg et al. 2021). The HD 163296 disk is thought to be hosting two planets in formation, identified by kinematic detections, at 94 au (Izquierdo et al. 2022) and 260 au (Teague et al. 2018a, 2021; Pinte et al. 2018b). HD 163296 presents extremely bright and radially extended molecular emission, which is pivotal to extract precise and accurate gas kinematics by leveraging the high signal-to-noise ratio (S/N) of specific lines.

2.2 Targeted emission lines

To extract information on molecular gas dynamics, we based our work on the publicly available high spatial and spectral resolution data provided by the MAPS ALMA Large Program (Öberg et al. 2021). The survey investigated more than 50 emission lines from different chemical species: in Table 1, we report the molecules and the corresponding emission lines that we considered for this analysis. An initial visual datacube exploration was performed for each molecule in order to select suitable lines: the candidates were chosen according to their brightness and high signal-to-noise ratio.

3 Methods

In this section, we give an overview of the general workflow of our analysis and describe in detail the methodology we followed.

Our main aim is to explore the kinematics of the HD 163296 disk through a multi-molecule approach: as different molecules emit at slightly different physical and thermal conditions, considering different emission lines is essential to trace the gas dynamics along the vertical extent of the disk. So far, HD 163296 has been kinematically analyzed through its brightest gas emission (Pinte et al. 2018b; Teague et al. 2018a, 2021; Casassus & Pérez 2019; Law et al. 2021b; Izquierdo et al. 2023) – ¹²CO and ¹³CO lines; as for abundant tracers with a high signal-to-noise ratio, it is more straightforward to perform a reliable reconstruction of the emitting layer and an exhaustive kinematic analysis. In this work, we expand the extraction of rotation curves and the following kinematic considerations also to different tracers, to obtain a more complete picture of gas dynamics in the disk. We show a schematic summary of the workflow:

• To extract and fit the rotation curves with a thermally stratified model, we assume the 2D thermal structure of the disk from Law et al. (2021b) and the molecular emitting layers ɀ(R) from Paneque-Carreño et al. (2023). We then perform a complementary test to verify the solidity of the extracted emitting layers, given the assumed thermal structure. The validation test consists of the extraction of molecular optical depth profiles τ(R) from the data: as optical depth is easily related to the height ɀ(R), checking that we obtain reasonable τ profiles is a good way to validate the solidity of the assumed emitting layers. To extract τ profiles, we proceed as follows: we consider the temperature structure T(R, ɀ) computed in Law et al. (2021b) through CO isotopologs and the molecular emitting layers from Paneque-Carreño et al. (2023); we extract the molecular optical depth profiles by comparing the surface brightness profile from the data with the expected profile in case of optical thick emission, considering the full Planck law. We assess that our test provides the expected results, in agreement with optical depth estimates obtained with independent methods;

• Once we have verified that the assumed ɀ(R) are reasonable using the optical depth profiles as benchmark, we extract rotation curves from the considered molecules, taking into account the corresponding emitting layers. Considering the large uncertainties on the assumed ɀ(R) from Paneque-Carreño et al. (2023) (for the faintest lines we consider, Δɀ/ɀ can be as high as ≈50%), we also check that rotation curves are not significantly influenced by the exact location of the emitting layers, so that the lack of accuracy in ɀ(R) does not severely impact our rotation curves analysis. For the hyperfine lines in our sample (i.e., HCN and C₂H), we extend existing methodologies to more efficiently fit the spectra and obtain improved rotation curves. We also analyze beam smearing effects on curves, verifying that it causes deviations opposite to pressure gradients;

• We fit the obtained rotation curves all together with an extended model, including contributions from the star, the disk self-gravity, the pressure gradient and vertical thermal stratification; from this fit we obtain a dynamical estimate for stellar mass, disk mass, and scale radius. Knowing these parameters and the 2D thermal structure, we show the retrieved 2D pressure structure of the disk. Finally, we compare the predictions from the stratified model and a simple isothermal disk structure, to recover which one can best describe the differences in rotation velocity between the considered molecules. In this way, we show the impact of the pressure structure on the kinematics of molecules different from ¹²CO and ¹³CO, and the improvement induced by the implementation of the vertical stratification, with respect to the isothermal model.

3.1 Molecular optical depth profiles

Law et al. (2021b) retrieved the emission layers of the three CO isotopologs – ¹²CO, ¹³CO, and C¹⁸O – within the radial extent where they are optically thick. In the optically thick regime, the gas temperature along the emitting layers can be easily measured from the peak surface brightness of the emission lines (Weaver et al. 2018); thus, from the lines profiles in each pixel it is possible to extract the gas kinetic temperature along the CO surfaces, assuming LTE conditions. They fitted the obtained temperature values with the 2D parametric expression for the kinetic temperature of the gas by Dullemond et al. (2020), which depends on the radius R and the height above the midplane ɀ: $$T^4(R,ɀ)=T_{\text{mid }}^4(R)+\frac{1}{2}\left[1+\tanh \left(\frac{ɀ-\alpha {ɀ}_q(R)}{{ɀ}_q(R)}\right)\right]T_{\text{atm}}^4(R).$$(1)

Here, ɀ_q(R) = ɀ₀(R/l00 au)^β, T_atm = T_atm,0(R/100 au)^q_atm , and T_mid = T_mid,0(R/100 au)^q_mid. They fitted seven parameters using data from the three CO optically thick isotopologs: T_atm,₀, q_atm, T_mid,0, q_mid, α, ɀ₀, and β. We report the best-fit parameters obtained by Law et al. (2021b) in Table 2. In this work, we use the kinetic temperature by Law et al. (2021b) and the observed brightness temperature from the data to derive radial profiles for the molecular optical depth.

We also consider H₂CO (3₀₃ − 2₀₂), HCO⁺ (J = 1 − 0), HCN (J = 3 − 2, F = 3 − 2), and C₂H (N = 3 − 2, J = 5/2 − 3/2, F = 3 − 2). We adopt the molecular emitting layers by Paneque- Carreño et al. (2023) into Eq. (1), to find the kinetic temperature of the gas at the specific layers. With the obtained temperature profile, we calculate the expected surface brightness profile I_Planck at the emitting layer, under the assumption of optically thick emission, through the Planck law: $$I_{\text{Planck }}(\nu ,T)=\frac{2h{\nu }^3}{c^2}\frac{1}{e^{\frac{h\nu }{k_{B}T(R,ɀ)}}-1}.$$(2)

Here, ν is the rest frequency of the considered line and T (R, z) is the temperature on the emitting surface computed from Eq. (1). Then, we extract the effective surface brightness I_data from the data, using the GoFish (Teague 2019b) function radial_profile, where we set a 30° wedge along the disk major axis. At this point, we know that the two obtained profiles for the surface brightness are linked by the following relation: $$I_{\text{data }}(R)=I_{\text{Planck }}(R)\left[1-e^{-\tau (R)}\right].$$(3)

By inverting Eq. (3), we can compute the optical depth radial profile τ(R) as $$\tau (R)=-\ln \left(1-\frac{I_{\text{data}}(R)}{I_{\text{Planck }}(R)}\right).$$(4)

In Fig. 1, we show the optical depth profiles of H₂CO, HCO⁺, HCN, and C₂H we obtained following this method. We calculated uncertainties on the τ profiles (shaded profiles in Fig. 1) propagating the errors on the extracted surface brightness. Uncertainties are generally larger for H₂CO and HCO⁺, the two lines with lower resolution data (0.3″).

First of all, we observe that H₂CO and HCO⁺ lines result in an optically thin emission along the whole radial extent of the disk, reaching maximum peaks at τ ∼ 0.11 and τ ∼ 0.45 respectively. C₂H emission becomes marginally thick, with a value for τ approaching to 1, at the radii corresponding to intensity peaks. HCN emission is mostly optically thick within ∼130 au, with higher values in the inner ∼60 au and a local maximum around the radius of the intensity peak (R = 109 au); it then evolves into optically thin emission, with τ decreasing by an order of magnitude. As expected, C₂H and HCN result in a marginally thick emission at the rings locations, as previously observed in other disks (Öberg et al. 2021, 2023). We also find that in the innermost tens au we measure a general decrease of the optical depth for all the molecules. This decrease is an observational artifact and is due to the beam dilution effect (Weaver et al. 2018): at smaller scales, the beam size may be larger than the spatial dimension of the considered emitting region of the disk. As a consequence, the measured surface brightness is underestimated by a beam dilution factor, proportional to the ratio of the emitting area to the beam size (lower than 1 at small radii).

The obtained estimates for τ are compatible with existent optical depth values measured for the same molecules in this (Bergner et al. 2021) and other disks (Bergner et al. 2019; Facchini et al. 2021), and are in line with the expected optical depth values that can be retrieved from column density profiles in Guzmán et al. (2021) and from the brightness temperature profiles computed by Paneque-Carreño et al. (2023) for HD 163296. Therefore, the consistent results of this preliminary test guarantee that, given the thermal structure for the disk by Law et al. (2021b), the emitting layers by Paneque-Carreño et al. (2023) are suitable for the rotation curves extraction. Considered the large errors on the assumed molecular emitting layers, we also verified that the rotation curves obtained with eddy are not significantly affected by the emitting layers choice, varying them within the given uncertainties. Thus the assumed layers are suitable to extract reliable curves, without the need to rely on deeply accurate measurements of ɀ(R), whose extraction procedure shows high systematic uncertainties and is thus intrinsically difficult.

Finally, we note that this is an efficient and straightforward method to easily extract optical depth and column density profiles of the considered line, relying only on thermal structure T(R, ɀ) and emitting height ɀ(R) of one single transition. This is an independent method that does not require complicated chemical studies about excitation temperatures and represents a valid alternative to approaches applied to hyperfine lines (Guzmán et al. 2021) for systems where a precise and accurate 2D temperature structure is known.

Fig. 1 Optical depth radial profiles for H2CO, HCO+, HCN, and C2H found in this work, taking into account the molecular emitting layers. Emission peaks and emission gaps, when present, are highlighted in red and black dashed lines, respectively. We also report the τ = 1 line for easier visualization.

3.2 Rotation curves extraction

Once we have verified with the optical depth complementary test described in the previous section that the assumed emitting layers are consistent, we proceed to the rotation curves extraction. In order to extract the velocity radial profiles of the HD 163296 disk, we used the eddy Python package (Teague 2019a), following the method presented in Teague et al. (2018a) and Teague et al. (2018b). First of all, we have to know the geometry of the disk (i.e., inclination i, position angle PA, possible disk center offsets in declination and right ascension) and the emitting layers of each considered molecule, to deproject the data from projected location in the sky to disk-frame coordinates. We assumed i = 46.7°, PA = 313.3° from Öberg et al. (2021) (the PA is shifted by 180° according to its definition in eddy) and center offsets equal to zero. For the molecular emission layers we assumed an exponentially tapered power law $$ɀ(R)=Z_0{\left(\frac{R}{100\text{au}}\right)}^{\varphi }\exp \left[-{\left(\frac{R}{R_{\text{taper }}}\right)}^{\psi }\right],$$(5)

with best-fit values by Paneque-Carreño et al. (2023), listed in Table 3, for Z₀, ϕ, r_taper, and ψ. In addition, we also performed our analysis adopting nonparametric emitting layers for all the molecules, and we show the results we obtained in Appendix E. In Sect. 4.2, we discuss these results, showing how the choice of using parametric emitting surfaces for the considered lines has no substantial effect on the qualitative results and trends we obtain from our analysis.

We note that taking into account the molecular emitting layers in the data deprojection does not influence the extracted curves in a dramatic way, for molecules emitting close to the midplane. Especially we are interested in testing the molecules with lower S/N emission, for which it is more difficult to extract reliable rotation curves (i.e., H₂CO, HCO⁺, C₂H). To test this, we computed the rotation curves of these molecules assuming flat emission (i.e., originating at the disk midplane) and we obtained that the differences between the two scenarios are compatible with zero within: 0.6σ ξ for HCO⁺ and H₂CO, 4.5σ for C₂H (and 3σ for C¹⁸O). For thoroughness, we did the same test for all the molecules: as expected, for the other lines the discrepancy is above 6σ, as these molecules are well known to emit from very elevated layers. For every line, we show the difference between the extracted rotation curves in case of flat and elevated emission in Fig. A.1. We also note that the curves would be even more compatible if we had more realistic estimates for the uncertainties on the rotation curves, as the errors on the curves are underestimated, as we explain in Sect. 4.1. This test guarantees a reliable extraction of rotation curves for lower S/N lines emitting close to the midplane, even if the available parameterizations for their emitting heights have large uncertainties, without affecting the quality of the results. In Sect. 4.2, we show the results on the stellar and disk parameters we obtain from the multimolecule fit of the rotation curves in case of flat emission for all the considered molecules.

To extract the rotation curves, we used the SHO (Simple Harmonic Oscillator) method in eddy, that assumes axisymmetric disks. Following this approximation, emission line spectra that arise from the same radial distance from the disk center are expected to exhibit the same profile, as they trace the same physical and chemical conditions. The line centroids are expected to be shifted by v_θ cos(θ), an offset with respect to the systemic velocity due to the azimuthal location of the emission in the disk (θ represents the polar angle, measured from the redshifted major axis). With this method, the best-fit rotation velocity at a specific radius is determined as the one that best aligns the deprojected spectra from the same annulus. Mathematically, we fit the set of line centroids of the same annulus with a cosine function. In order to align the spectra, in our analysis we fit each spectrum with a Gaussian curve to identify the intensity peak and the corresponding velocity, and then we proceed to the fit of the line centroids with the harmonic function. For each tracer we repeat the extraction of the rotation curve 20 times, considering different independent samples of spectra inside each annulus for different iterations and we compute the uncertainties on the velocity values as their weighted standard deviation. A similar approach is followed in Teague et al. (2021) to extract the ¹²CO rotation curve for the HD 163296 disk with GoFish: in that case, the best-fit velocity values are selected as those that best align the spectra, as quantified by fitting a Gaussian processes mean-model to the averaged spectra.

If we want to consider also the radial contribution to the line-of-sight velocity v_los , in eddy we can fit the azimuthal and radial velocities together, and fitting every set of line centroids with a double harmonic oscillator. We note that for molecules other than the CO, the MAPS quality of data is not sufficient to provide reliable kinematical estimates of the radial and vertical velocity components. Nevertheless, in our analysis we are interested in the azimuthal velocity profiles. We decided not to focus on v_R and v_ɀ, considering them negligible with respect to the difference in azimuthal velocities between data and an axisymmetric model, as shown by Teague et al. (2019) and Izquierdo et al. (2023).

3.2.1 Limitations of the method

As explained in Sect. 3.2, the extraction of rotation curves is based on the identification and fit of line centroids in each pixel of the data cube. Besides requiring a high signal-to-noise ratio, this methodology relies on a few assumptions. First, we assumed that the disks are azimuthally symmetric and that the spectra along each annulus vary only by a phase cos θ, but maintain the same profile. This hypothesis might be inadequate for highly perturbed line emission morphologies. Second, eddy does not distinguish the upper and lower surfaces of disks, and assumes that the front side dominates over the back side. However, an important contamination from the back side of the disks can complicate the spectra alignment and so the goodness of the fit. In some cases (such as GM Aur in the MAPS sample, see Lodato et al. 2023) the signals arising from the two surfaces are comparable and thus the procedure results in a quite inefficient alignment of the spectra. In these cases, the lower surface can shift the line centroids in the observed intensity spectra, inducing a wrong derivation of the best-fit velocities. We note that the level of contamination from the back side depends on the considered disk and molecule: for example, line emission coming from lower layers is generally less contaminated by the back side with respect to the emission from the upper layers, since its emission surface is closer to the disk midplane (Lodato et al. 2023).

To quantify the effects of these limitations on our analysis, we quantitatively compared the rotation curves we extracted with eddy with the ones obtained by Martire et al. (2024) with the discminer tool (Izquierdo et al. 2021) for ¹²CO and ¹³CO, considering also the back side of the disk. We measure very low differences, which are <4% for ¹²CO and <5% for ¹³CO, confirming that the limitations of this method are not our primary source of systematic effects.

Fig. 2 Teardrop plots for C2H (left) and HCN (right). Each row of the plot represents an averaged spectrum at separate annuli after shifting and stacking. The color scale indicates the intensity of the emission line. In both cases, hyperfine transitions (two for C2H and three for HCN) are clearly visible.

3.2.2 Multi-Gaussian methods for hyperfine transitions kinematics: C₂H and HCN rotation curves

In this work, we developed a new method to extract rotation curves from molecules that have hyperfine transitions. In fact, in our analysis molecules such as C₂H and HCN do exhibit a hyperfine structure (see Table 1 for details about the transitions). Hyperfine structure is typically characterized by energy shifts that are orders of magnitude smaller than fine-structure ones. Therefore, for these molecules, it is possible to detect two (or more) emission lines in the same frequency channel (i.e., in the same channel map). Typically, a spectrum from molecules with hyperfine lines shows two or more different intensity peaks (one for each transition) instead of one, as clearly shown by teardrop plots in Fig. 2, where the color scale represents the averaged intensity of the considered emission line along separate annuli. The intensity of each visible transition depends on its Einstein coefficient: for C₂H, in each spectral channel we see two similar transitions, roughly equal in intensity. This reflects in double-peaked spectra with two maxima at the same intensity. For HCN, instead, we can clearly detect a central dominant transition, which is the strongest one, alongside with two fainter lines nearly of the same intensity. Therefore, the standard extraction methods (that fit a quadratic or Gaussian function to the spectrum) are not adequate to fit these multi-peaked spectra, and they fail to return reliable velocity profiles.

For this reason, we implemented in eddy a new fitting procedure for hyperfine lines using multiple Gaussian components (two or three, depending on the number of hyperfine transitions), where the width of the multiple Gaussians is assumed to be the same, as the hyperfine components are expected to emit under the same conditions of temperature, density, and turbulence.

For C₂H, we fit each spectrum with a double-Gaussian function $$I=I_0\exp \left[{\left(\frac{v-v_0}{\sigma }\right)}^2\right]+I_0\xi \exp \left[{\left(\frac{v-\left(v_0-\text{Δ}v\right)}{\sigma }\right)}^2\right],$$(6)

where ξ is the ratio of the two intensity peaks and Δv is the velocity difference between the two transitions. We fit four free parameters: the center of the Gaussian v₀, which represents the line centroid, the width of the Gaussians σ, the intensity corresponding to the first peak I₀ , and the peak intensity ratio ξ.

For HCN, we fit each spectrum with a triple-Gaussian function $$\begin{array}{c}I=I_0{\xi }_{\text{blue }}\exp \left[{\left(\frac{v-\left(v_0-\text{Δ}{v}_{\text{blue }}\right)}{\sigma }\right)}^2\right]+I_0\exp \left[{\left(\frac{v-v_0}{\sigma }\right)}^2\right]\\ +I_0{\xi }_{\text{red }}\exp \left[{\left(\frac{v-\left(v_0+\text{Δ}{v}_{\text{red }}\right)}{\sigma }\right)}^2\right],\end{array}$$(7)

where ξ_blue (ξ_red) is the ratio of the blueshifted (redshifted) intensity peak to the central dominant peak and Δv_blue (Δv_red) is the velocity difference between the central and the blueshifted (red- shifted) peak. In this case we fit five free parameters: the center of the middle Gaussian v₀, which represents the line centroid, the width of the Gaussians σ, the intensity corresponding to the central peak I₀ , and the peak intensity ratios ξ_blue and ξ_red .

The velocity shifts Δv (between the two C₂H hyperfine lines), Δv_blue, and Δv_red (between the weaker hyperfine lines and the central one of HCN) are fixed parameters, given the known frequency shifts Δν = 2.5 MHz (between the two C₂H transitions), Δν_blue = 2.09 MHz, and Δν_red = 1.54 MHz (between the HCN transitions) as taken from the Cologne Database for Molecular Spectroscopy (CDMS) (Sastry et al. 1981; Müller et al. 2000 for C₂H and Ahrens et al. 2002 for HCN). We set priors for some parameters:

• Intensities must be positive: I₀ , ξ, ξ_blue , ξ_red > 0;

• ξ ∈ [0.2, 1.5] for C₂H and ξ_blue, ξ_red ∈ [0.1, 1.] for HCN: we choose a reasonably large interval centered on the expected value for the ratio between the intensity peaks, considering the Einstein coefficients and upper state degeneracies of the considered transitions $(\xi ~A_{ij}^1{𝑔}_1/A_{ij}^2{𝑔}_2\approx 0.65$ for C₂H and $\xi ~A_{ij}^1{g}_1/A_{ij}^2{g}_2\approx 0.13$ for HCN);

• σ ∈ [0, Δv/2] for C₂H and σ ∈ [0, Δv_blue/2] for HCN: the width of the Gaussians is not larger than half of the distance in velocity between the two considered transitions, as the expected thermal broadening of the line is below ∼1 km/s (Izquierdo et al. 2023) and thus much smaller than Δv.

For C₂H, we set reasonable starting positions for the parameter space exploration, to ease the convergence of the fit. In particular, we computed the Keplerian velocity at each radius R and azimuthal angle θ as $$v_{\text{rot},\text{K}}={\left(\frac{G{M}_{\star }}{R}\right)}^{1/2}\sin i\cos \theta +v_{\text{sys}},$$(8)

where M_★ is the stellar mass computed fitting the ¹²CO curve with a Keplerian model, and used it as starting position for the line centroid v₀ fit.

In Fig. 3, we show the fit of the hyperfine spectra of both C₂H and HCN with, respectively, the double- and tripleGaussian methods, compared to the single-Gaussian one originally implemented in eddy. These methods lead to a remarkable improvement in the extraction of the rotation curves of hyperfine structured lines. In Fig. 4, for both molecules we show the comparison between the rotation curve extracted with the multiGaussian method and the curve previously obtained with the single-Gaussian method already implemented in eddy.

These methods prove very efficient in fitting the multipeaked spectra in both cases. For C₂H, the remarkable upgrade in the peak fitting procedure leads to the most significant improvement in the azimuthal velocity profile: fitting the spectra with a double-Gaussian, we obtain a smooth curve with small errors, a more reliable result than the one obtained using a singleGaussian. In that case, since the two peaks in the C₂H spectra are mostly equal in intensity, a single-Gaussian model leads to a fragmented curve affected by high variance. For HCN, by fitting three Gaussian components to each spectrum we maximize the exploitation of the available data, taking into account also the two minor intensity peaks. Nevertheless, the two curves do not differ substantially and the observed differences in the selected radial range are small: as the central peak is dominant, also the standard single-Gaussian method, even if less accurate than the triple-Gaussian, could still fit the central prominent peak providing the correct value of the line centroids, neglecting the hyperfine splitting visible in the two minor peaks.

Fig. 3 Improvement in the fitting of the spectra using multi-Gaussian methods. Top panel: comparison between the single-Gaussian (light blue) and the double-Gaussian (pink) methods used to fit the intensity peaks of C2H at 50 au. Bottom panel: same as top panel, but for the three hyperfine components of the HCN transition.

Fig. 4 Comparison between the rotation curves for C2H (top) and HCN (bottom) using the double- or triple-Gaussian methods (pink) and the single-Gaussian method already implemented in eddy (dashed light blue) for the fitting of the spectra.

3.2.3 Beam smearing effect

Among the considered molecules, C₂H and HCN show strong intensity gradients in their emission. Since the extraction of velocity profiles comes from molecular line emission, the presence of gradients in the emission intensity of these lines can affect the corresponding rotation curves. In Fig. 5, we show these observational effects for C₂H; we show the same effect for HCN in Appendix B.

We find a strong anticorrelation between the two following measurable quantities as a function of radius (evaluated at the same radii), in agreement with the results from Keppler et al. (2019):

• d ( ∫ Idv)/ dR: the intensity gradient of C₂H molecular emission line, computed from the radial intensity profile¹;

• $v_{{\text{C}}_2\text{H}}-v_{{}^{12}\text{CO}}$ : the difference between the rotation velocities of C₂H and ¹²CO (shown separately in panel b of Fig. 5). In this case, the ¹²CO curve was chosen as reference because it exhibits a typical smooth Keplerian profile, and it does not show any strong gradient in the intensity radial profile.

These two quantities display an opposite modulation, as shown in panel c of Fig. 5. In particular, the difference in rotation velocity $v_{{\text{C}}_2\text{H}}-v_{{}^{12}\text{CO}}$ is positive at the inner edge of an emission gap, while it becomes negative at its outer edge. This effect is shown in the comparison between the curves reported in panel b of Fig. 5. This particular modulation can be explained by a beam smearing effect. When a source is observed at a given angular resolution, every extracted spectrum is the intensity weighted average of all the spectra that fall inside the beam area. For this reason, a strong intensity gradient across the beam would bias the recovered spectrum and, consequently, the extracted velocity. In our results we are biased to see velocities corresponding to higher intensity values inside the beam, that do not correspond to the beam center. For instance, in correspondence of the inner edge of the emission gap we are affected by a strong negative intensity gradient. So in this region, the extracted spectrum weighs more the velocities corresponding to the higher intensities, which are the ones from the inner radii, and we measure an artificially higher rotation velocity. On the contrary, at the outer edge of the emission gap the strong intensity gradient is positive, and velocities corresponding to the outer radii are weighted more and we measure an artificially lower rotation.

Beam smearing can also explain the flattening at small radii we see in the C₂H rotation curve in panel b of Fig. 5, since this molecule shows a central cavity in emission, followed by an intensity peak. In the internal cavity, the signal inside the beam is dominated by the higher intensities found at larger radii: as a consequence, we measure lower velocities in correspondence of the cavity, leading to a flattening of the curve.

This observational effect has important consequences: we have observed that beam smearing triggers apparent velocity deviations in rotation curves, that are linked to strong intensity gradients of the emission line. Nevertheless, if any substructure (such as rings or gaps) is present in the gas surface density, then also the pressure profile undergoes alterations induced by the change in density, causing similar velocity deviations. Therefore, it is important to distinguish rotation curve modulations that are due to this observational bias from effective deviations linked to pressure gradients, that could be signs of the presence of gas density substructures. However, we note that beam smearing works in an opposite way to pressure gradients, as we show in Fig. 6: assuming that intensity peaks trace pressure peaks (and analogously with gaps), pressure gradients are expected to induce velocity deviations that are of opposite sign with respect to the ones deriving from intensity gradients that we observe. This means that the apparent difference in velocity induced by beam smearing in correspondence of an intensity gap could be incorrectly interpreted as a pressure maximum, not as a pressure gap as expected if intensity and pressure peaks correspond. Conversely, if we observe kinematic evidence of a pressure gap in correspondence of an intensity gap, we can immediately conclude that it cannot have been artificially generated by beam smearing, since this effect would have induced the opposite behavior (and analogously for a pressure peak). This different impact that beam smearing and pressure gradients have on rotation curves may help distinguish these effects; nevertheless, this may be not straightforward, since both of them can be present and influence data in different ways.

Fig. 5 Beam smearing effect for C2H. Panel a: azimuthally averaged integrated intensity profile of the C2H line. Emission peaks are marked with dashed red lines and the emission gap is highlighted with a dashed grey line. Panel b: rotation curve of C2H (purple line) compared to that of 12CO (green line). The rotation velocity of C2H is super-Keplerian just before the gap and sub-Keplerian after. It shows the opposite behavior in correspondence of the emission peaks. Panel c: anticorrelation between the C2H integrated intensity gradient and its difference in rotation velocity compared to 12CO. The two quantities show an opposite modulation, and they cross at the radii corresponding to the emission gap or peaks.

Fig. 6 Sketch of the different effects that beam smearing due to line intensity gradients (pink dashed line) and additional pressure gradient effects due to the presence of gas density substructures (light blue solid line) have on a standard rotation curve (black solid line). We assume that intensity gradients trace gas density gradients; a peak or gap in line intensity corresponds to a peak or gap in gas density. Line intensity and gas density peaks are marked with dashed red lines, and line intensity and gas density gaps are highlighted with a dashed grey line. At the bottom we show a colored map of the radial intensity. The solid black line represents the rotation curve of a molecule showing no gas density or line intensity substructures. The typical curve distortion due to beam smearing effect is shown by the pink dashed line, and is in agreement with our observations. The expected distortion of the curve due to the presence of gas density substructures is shown by the light blue solid line. The rings and gaps in the gas density induce additional pressure gradients that influence the rotation curve, according to the $v_{\varphi }^2(R,ɀ)$ expression shown in the upper right corner. These pressure gradients are positive before rings and beyond gaps locations, while they are negative before gaps and beyond rings locations. Thus, in this case, the deviations in rotational velocity are expected to be in the opposite direction of the beam smearing-induced variations.

3.3 Vertically stratified model

We extract and analyze multiple molecular lines emitting from different layers to investigate how disk kinematics is affected by height above the midplane across the vertical extent of the disk. Recent high resolution gas observations have shown differences in the rotation velocity of the CO isotopologs, that could not be correctly explained with simple vertically isothermal disk models. In Martire et al. (2024), the authors show that a thermal stratified model can better explain the velocity discrepancies between the ¹²CO and ¹³CO rotation curves, and that vertical stratification should be considered when building a solid model to adequately describe them. Therefore, we choose to simultaneously fit our seven rotation curves including the contribution from vertical thermal stratification, in order to provide a complete picture of the several contributions that can affect the disk kinematics. From this fit we derive dynamical estimates for the disk and stellar masses, and the disk scale radius; we also compare our results with a simple vertically isothermal model, showing that it is not adequate to describe the complexity of the data.

We consider a disk surface density as described by the selfsimilar solution of Lynden-Bell & Pringle (1974): $$\text{Σ}=\frac{(2-\gamma )M_{\text{d}}}{2\pi R_{\text{c}}^2}{\left(\frac{R}{R_{\text{c}}}\right)}^{-\gamma }\exp \left[-{\left(\frac{R}{R_{\text{c}}}\right)}^{2-\gamma }\right].$$(9)

Here M_d and R_c are the disk mass and the scale radius, R is the cylindrical radius, and γ is the index that describes the steepness of the surface density (we adopt γ = 1). On the midplane we have that the temperature can be described by a power law as $T_{\text{mid }}(R)\propto R^{q_{\text{mid }}}$ and consequently the sound speed is $c_{\text{s, mid}}\propto \sqrt{T_{\text{mid}}}\propto R^{q_{\text{mid}}/2}$. To account for the thermal vertical stratification of the disk, following the method and parameterization described in Martire et al. (2024), we can define two functions f and 𝑔 that determine the vertical dependence of the disk temperature, sound speed, and density: $$T(R,ɀ)=T_{\text{mid }}(R)f(R,ɀ),$$(10) $$c_{\text{s}}^2(R,ɀ)=c_{\text{s},\text{mid}}^2(R)f(R,ɀ),$$(11) $$\rho (R,ɀ)={\rho }_{\text{mid }}(R)𝑔(R,ɀ).$$(12)

As a consequence, considering a barotropic fluid, the 2D pressure in the disk is given by $$P(R,ɀ)=P_{\text{mid }}(R)f(R,ɀ)𝑔(R,ɀ)=c_{\text{s},\text{ mid }}^2(R){\rho }_{\text{mid }}(R)f𝑔(R,ɀ).$$(13)

We note that 𝑔(R, ɀ) is a function that only depends on f (R, ɀ) as a consequence of the hydrostatic equilibrium (see Martire et al. 2024).

We adopt the model for our rotation curves fit from Martire et al. (2024), given by $$\begin{array}{c}{v}_{\theta }^2(R,ɀ)=v_{\text{k}}^2\{{\left[1+{\left(\frac{ɀ}{R}\right)}^2\right]}^{-3/2}-[{\gamma }^{\prime }+(2-\gamma ){\left(\frac{R}{R_{\text{c}}}\right)}^{2-\gamma }\\ -\frac{\text{d\hspace{0.17em}}\log (f𝑔)}{\text{d\hspace{0.17em}}\log R}]{\left(\frac{H}{R}\right)}_{\text{mid }}^{2}f(R,ɀ)\}\\ +G{{{\displaystyle \int }}^{\text{}}}_0^{\infty }[K(k)-\frac{1}{4}\left(\frac{k^2}{1-k^2}\right)\\ \times \left(\frac{r}{R}-\frac{R}{r}+\frac{{ɀ}^2}{Rr}\right)E(k)]\sqrt{\frac{r}{R}}k\text{Σ}(r)dr,\end{array}$$(14)

where $v_{\text{k}}=\sqrt{G{M}_{\star }/R}$ is the Keplerian velocity, ${\gamma }^{\prime }=\gamma \text{+(3+ }q\text{) }/2,r=\sqrt{R^2+{ɀ}^2}$ is the spherical radius, K(k) and E(k) are complete elliptical integrals (Abramowitz & Stegun 1970), and k² = 4Rr/[(R + r)² + ɀ²]. In this expression we can distinguish the different contributions of the stellar gravity, the pressure gradient, the thermal stratification (the term depending on f and 𝑔), and the self-gravity (Bertin & Lodato 1999). For further details on this model, see Martire et al. (2024). The adopted expression for the isothermal model is explained in Lodato et al. (2023), Eqs. (12)–(14).

To perform the simultaneous fit of the rotation curves, we use the code DySc²: given the data, the thermal structure of the disk, and the geometry of the emitting layers, the code performs the fit of the data to the model by maximizing the likelihood of the model, with the posterior distribution explored by Monte Carlo Markov Chains, as implemented in emcee³. As output, we obtain the probability distributions of the fitted parameters, their best-fit values, and corresponding uncertainties. We implemented the possibility to fit multiple (>2) rotation curves in the code. We fit our seven rotation curves obtained from data with a single model for the rotation velocity described in Eq. (14). From this multi-molecule fit, we can extract accurate best-fit values for some fundamental parameters: stellar mass M_✶, disk mass M_d, and scale radius R_c.

4 Results and discussion

4.1 Multi-molecule fit with thermal stratification: dynamical estimates and pressure structure

We extended the code DySc to perform a multi-molecule fit with seven different molecules and we fitted the obtained molecular rotation curves simultaneously with one single model, including contributions from the stellar and the disk gravity, the pressure gradient and thermal stratification, as described by Eq. (14). We ran the MCMC exploration using 200 walkers, 1000 burn-in steps and 2000 steps to determine the final posterior distribution of the fitted parameters. To evaluate the contributions of vertical structure, we took into account the molecular emitting layers found by Paneque-Carreño et al. (2023). For the thermal structure to use in our fit, we considered the parameterization for the temperature presented by Dullemond et al. (2020) and reported in Eq. (1), with the parameters obtained by Law et al. (2021b) reported in Table 2. In this work we did not account for the uncertainties in the temperature structure of the disk. Andrews et al. (2024) analyzed the systematic uncertainties driven by the disk thermal structure, with related uncertainties that are of the same order of surfaces mis-specifications.

We report the multi-molecule fit in Fig. 7: we show the fit of each molecular rotation curve with the recovered bestfit model. It is clear that the model well reproduces the seven curves, with the exception of a short range of radii in H₂CO around R ~ 150 au: in a short interval around this radius, which represents an emission gap for H₂CO, the reconstruction of the rotation curve with eddy is not accurate, since the signal-tonoise ratio is too low to allow an efficient extraction of the rotation velocity (see Appendix C for further details). As for HCN and C₂H, the differences with respect to the model are due to the apparent deviations in the rotation curves caused by beam smearing, as explained in Sect. 3.2.3. It is important to note that the inclusion of several molecules different from the CO isotopologs, with their corresponding emitting heights, does not hinder the successful outcome of the fit, but rather it contributes to its accuracy.

The uncertainties on the rotation curves are very small: this is an underestimate as we are taking into account only the statistical errors on the values, but no systematic errors. To provide a more exhaustive estimate of the uncertainties, it would be necessary to implement in eddy the systematic component of the errors on the derived rotation curves associated with the uncertainty on the geometrical parameters of the disk (inclination, position angle, and centering), the thermal structure, and the molecular emitting layers. Andrews et al. (2024) discuss an efficient way forward in the determination of systematic uncertainties in kinematic measurements.

From the multi-molecule fit, we are able to find the bestfit values for the stellar mass M_⋆ , the disk mass M_d, and scale radius R_c. We note the potential of detailed kinematical analysis in extracting dynamical estimates of these crucial quantities, whose estimate may be very difficult with other methods. For each parameter, we compute a 95% credibility interval (that is the Bayesian counterpart of the confidence level) as the range of values that contains 95% of the probability density in the posterior distribution and we consider the extremes of these intervals as estimates of the best-fit parameters uncertainties. We summarize the best-fit values we obtained from the multimolecule fit, including the thermal stratification contribution, in Table 4; for completeness, we report the corresponding corner plots in Appendix D. We compare these estimates with the values extracted from a multi-molecule fit, this time assuming an isothermal temperature profile T(R) ∝ R ^q with H/R = 0.084 and q = −0.84 (Zhang et al. 2021). We also compare our results with the ones obtained by Martire et al. (2024): the authors retrieve the rotation curves using discminer⁴, and use only ¹²CO and ¹³CO for the fit, including thermal stratification contribution.

We note that in our work the estimate of the stellar mass does not change between the isothermal and stratified fits, while the disk mass appears slightly lower in the second scenario. What changes the most is the disk scale radius, as expected: in the isothermal case, we obtain R_c = 48.1 au, a very low value for the scale radius and a clear sign that we are probably underestimating it. On the contrary, with the thermally stratified model we provide a more accurate characterization of the thermal structure of the disk and this leads to R_c = 143.1 au. This value is more consistent with the best-fit estimates of R_c obtained with complementary methods by de Gregorio-Monsalvo et al. (2013) and Zhang et al. (2021). When comparing our results for the thermally stratified scenario with the ones from Martire et al. (2024), we find that they obtain slightly higher values for both the stellar and the disk masses, and a lower estimate of the disk scale radius. We point out that the fitting procedure in the two works is the same; the main differences are represented by the rotation curves to fit (the authors use discminer by Izquierdo et al. 2021 to extract the velocity profiles, while we use eddy) and the number of considered emission lines (they use ¹²CO and ¹³CO while we use seven different molecules). We also observe that we obtain uncertainties on the best-fit parameters that are an order of magnitude smaller than the ones obtained by Martire et al. (2024). This is understandable considering that eddy and discminer are based on different methods to provide rotation curves, so we do not expect to obtain the same order of magnitude for the uncertainties a priori. Error bars on the determined velocity values affect the simultaneous fit performed on the curves with DySc: as we are underestimating the uncertainties on the rotation curves with eddy, consequently we obtain low estimates of the errors on the best-fit parameters.

To emphasize the importance of the multi-molecule approach, we fitted the rotation curves starting with two molecules and then progressively adding one molecule at a time. All the fits were performed with the same number of walkers, burn-in steps, and steps (200, 1000, 2000). To quantify the goodness of the fit, we show in Fig. 8 the ${\chi }_{\text{red }}^2$ obtained from all the fits, increasing one by one the number of fitted molecules. We normalized the ${\chi }_{\text{red }}^2$ to the maximum value, since the uncertainties on the extracted rotation curves are systematically underestimated (as explained above in this section).

As we can see from Fig. 8, the most significant improvement is achieved with the inclusion of C¹⁸O in the fit, going from two to three molecules. This addition also induces the largest improvement in the dynamical estimates of the disk mass and scale radius: the values we obtain by fitting two and three molecules differ by ≲10% and ≲6% respectively for M_d and R_c. Moreover, the ${\chi }_{\text{red }}^2$ keeps decreasing with the further inclusion of lines in the multi-molecule fit. We added HCN and C₂H as last, since their rotation curves are the most affected by the beam smearing effect, as shown in Sect. 3.2.3. We also found that the width of the posterior distribution shrinks as expected when fitting more rotation curves together, but our error bars on the velocity profiles are so small that they are still dominated by the systematics.

Knowing the disk mass and scale radius from the multimolecule fit, and considering the thermal structure of the disk by Law et al. (2021b), we can show the normalized 2D pressure structure of the disk in Fig. 9, following Eq. (13). We note that the shape of the P(R, ɀ) profile is given by the assumed disk temperature T(R, ɀ), while the normalization term P_mid(R) is determined by the values for the disk mass and scale radius obtained by the fit of the rotation curves. We also show the location of the emitting layers of the tracers considered in the disk: the selected molecules probe the 2D pressure structure of the disk across its (R, ɀ) extent, from the upper surface down to the midplane.

Martire et al. (2024) showed that thermal stratification plays an important role in explaining the discrepancies between ¹²CO and ¹³CO rotation curves. To do so, they compute the difference in the squared rotation velocities of the two molecules, with a proper normalization; this quantity strongly depends on the different heights of the tracers, and it is thus very useful to assess the relevance of thermal structure in disks. Hence, they show the expected behavior of this quantity according to both the isothermal and the stratified models, to see which one can better reproduce the observed data. In this work, we perform the same calculation for all the considered molecules at different emitting heights, using our multi-molecule stratified model, to establish if thermal stratification has a strong impact on gas kinematics across the vertical extent of the disk. We compare squared rotation velocities of each molecule with the ¹²CO one, and we illustrate the expected behavior of their normalized difference in the isothermal and in the vertically stratified scenarios. We show in Fig. 10 the results we obtained for C¹⁸O, while in Appendix F we report the results for the other five molecules.

From our results, it is clear that the thermally stratified prediction is more appropriate to describe the observed differences between the two rotation curves, while the isothermal model cannot reproduce the observed behavior. We observe that the effect of thermal stratification is important also in the lower layers located closer to the disk midplane, as shown in Figs. 10 and F.1. Temperature leaves a clear imprint in the kinematics of molecular lines originating from the whole vertical extent of the disk, in the upper layers (as proved by Martire et al. 2024) and in the layers close to the midplane. Once this is included in the model, this further large-scale modification to the pressure can be distinguished and measured efficiently.

Considering Fig. 10, we point out that the difference in rotation velocity between the two isotopologs increases with radius and then starts decreasing at ~200 au; this trend is well reproduced by the stratified model, contrarily to the isothermal one. This discrepancy between the two models we considered for the rotation curves lies in the pressure gradient term, which includes an additional contribution depending on the height above the midplane when accounting for thermal vertical structure: $$\begin{array}{c}\frac{R}{\rho }\frac{dP}{dR}(R,ɀ)=-v_{\text{k}}^2[{\gamma }^{\prime }+(2-\gamma ){\left(\frac{R}{R_{\text{c}}}\right)}^{2-\gamma }\\ -\frac{\text{d\hspace{0.17em}}\log (f𝑔)}{\text{d\hspace{0.17em}}\log R}]{\left(\frac{H}{R}\right)}_{\text{mid }}^{2}f(R,ɀ).\end{array}$$(15)

In particular it can be verified that, for the thermal structure we assumed (i.e., the selected functions f and 𝑔 in Eqs. (10), (11)), this term remains negative across the whole radial extent of the disk and it decreases faster for ¹²CO than for C¹⁸O as a function of radius, up to ~200 au. From this point onward, the pressure contribution keeps decreasing only for C¹⁸O while it remains approximately constant for ¹²CO. Physically, this translates into a stronger effect of the pressure gradient in the higher layers with respect to the lower ones, when vertical thermal stratification is taken into account. For this reason, the gas rotation is slowed down at elevated heights more than near the midplane, because of the larger negative contribution of the pressure gradient. This behavior depends on the thermal structure of the disk and we only tested it for the temperature prescription we considered; nevertheless, we expect to observe the same result when describing the disk thermal structure with a similar function, both radially and vertically monotonic. The difference between the predictions of the two models becomes less significant when the emitting layer of the considered molecule is closer to the ¹²CO one: this is clear from the results for HCN shown in Appendix F.

In Fig. 10, we also see an oscillating behavior of the difference of the velocities squared. These modulations are likely correlated to the radial pressure substructure exhibiting peaks and troughs, as imprinted in molecular rotation curves (Teague et al. 2019; Izquierdo et al. 2023). As we are not accounting for pressure modulated substructures in our model, we are not able to recover the oscillating behavior of the measured velocity difference. An important improvement would be to include these small-scale pressure oscillations in order to perform a more exhaustive fit of the curves. It is worth highlighting that these additional pressure effects are seen in the difference of the velocities squared, extracted from two different emission heights, which thus suggests different amplitudes of the pressure substructure at distinct heights in the disk. This indicates that rotation curves can be a powerful tool to characterize the 3D structure of gaps, which should be further investigated in the future.

In summary, we can use multi-molecule emission to investigate and directly probe the 2D pressure structure of disks, one of the main physical properties that shape disks dynamics; we can measure the pressure effects on a global scale (i.e., extended pressure gradients along the radial and vertical range) and on a local scale (i.e., localized pressure-modulated substructures).

Fig. 7 Multi-molecule fit of the seven considered rotation curves with the thermally stratified model. For each tracer, the red dots with error bars represent the data, while the blue lines are the best-fit models of each velocity profile.

Fig. 8 Measured ${\chi }_{\text{red }}^2$ (normalized to the maximum value) of the multimolecule fits. We changed the number of molecules considered in the fit, adding the lines one by one.

Fig. 9 Retrieved 2D pressure structure P(R, ɀ) as a function of radius and height across the disk. The location of the selected molecular emitting layers are marked with colored solid lines, and highlight that we are tracing the pressure along an extended radial and vertical range in the disk.

Fig. 10 Normalized difference between the squared rotation velocities of C18O and 12CO: comparison between the data (light blue dots), the prediction from isothermal model (green line), and the prediction from thermally stratified model (purple line).

4.2 Dependence on the emitting layers

In this section, we analyze how changing the considered molecular emitting layers can impact our results on the dynamical estimates of M_⋆, M_d, and R_c.

As a first test, we repeated the whole identical routine of rotation curves extraction and multi-molecule fit of the curves, considering flat emission for all the molecules. We note that this scenario is extreme and does not well represent the real system, as we are assuming that all the molecules emit from the midplane, which is clearly not realistic as ¹²CO, ¹³CO, and HCN are known to emit from very elevated layers. In this case, we obtain the following best-fit estimates: M_⋆ = 1.7827 ± 0.0003 M_⊙, M_d = 0.1396 ± 0.0003 M_⊙, and R_c = 174.8 ±0.3 au. Thus, the differences in the final dynamical estimates of the parameters are within 10% for M⋆, 18% for M_d, and 20% for R_c. We note that these differences are mainly driven by the overestimation of the rotation curve of ¹²CO in the outer regions of the disk; in fact, we fit a lower value for M_⋆ and a higher value for M_d, which is strongly related to the disk outer regions. However, in this case we are assuming ¹²CO is emitting from the midplane and we know that this scenario is unrealistic. However, a much more realistic scenario where we consider flat emission only for lower S/N lines and/or molecules emitting close to the midplane (HCO⁺, H₂CO, C₂H, and C ¹⁸O) would impact the results far less, as the extracted curves in the flat and nonflat scenarios are compatible (as explained in Sect. 3.2 and shown in Fig. A.1).

As a second test, we repeated the whole analysis adopting nonparametric emitting layers for all the molecules, as retrieved by Paneque-Carreño et al. (2023). In Appendix E we show the corner plot obtained from the multi-molecule fit, with the dynamical estimates of the stellar mass, disk mass, and scale radius. The best-fit values we obtain from the fit in the nonparametric case differ from the parametric case results by ≲1% for M_⋆, ≲30% for M_d, and ≲13% for R_c. We see that the shape of the emitting surface has an impact mainly on the disk mass, that is affected by the systematics described by Veronesi et al. (2024) and Andrews et al. (2024). The difference we measure for the disk mass estimate is in agreement with the results by Veronesi et al. (2024) and Andrews et al. (2024), predicting an uncertainty close to ~25%. Also, in Appendix E we show the same plot as Fig. 10 for the nonparametric case, in which we recover exactly the same trends as the parametric case; thus, the same conclusions we drew in Sect. 4.1 also apply to the nonparametric scenario.

5 Conclusions

In this work, we explored the potential of a multi-molecule approach to gas kinematics to accurately characterize the 2D (R, ɀ) pressure structure of disks. We focused our analysis on the HD 163296 protoplanetary disk, and investigated the 2D structure by analyzing seven molecular tracers originating at different heights across the disk vertical extent to gain insights into the physical properties across the distinct disk layers. Considering the disk 2D temperature structure as estimated through the optically thick emission layering by Law et al. (2021b) and the emission heights by Paneque-Carreño et al. (2023), we extended existent tools to extract the rotation curves along these layers. We then simultaneously fitted all the curves for the first time with a multi-molecule (≥2) model of a disk including thermal stratification. Our results show the following:

• The 2D pressure structure leaves a measurable kinematical imprint in the gas rotation curves that is important across the whole disk vertical extent, also close to the midplane. Therefore, thermal stratification must be included in models to correctly reproduce the measured velocity differences at distinct heights, extending the approach by Martire et al. (2024) closer to the midplane;

• There are pressure substructures linked to the existence of dust gaps (not considered in our model) that can be seen from the rotation curves analysis and clearly exhibit a vertical dependence;

• A multi-molecule fit (≥2) of the rotation curves including the contribution from vertical stratification can be performed efficiently and dynamical estimates can be constrained more accurately (within the limits of our parametric method) for the stellar mass (M_★ = 1.89 M_⊙), disk mass (M_d = 0.12 M_⊙), and disk scale radius (R_c = 143 au), in agreement with Martire et al. (2024);

• The main results we obtain do not depend significantly on the chosen parameterization for the molecular emitting layers. The disk mass is the parameter that is more affected by these systematics (≲30%);

• We can retrieve the 2D pressure structure of the disk by knowing the 2D temperature profile, the disk mass, and the scale radius from the rotation curves fit;

• Spectra from molecular lines with hyperfine splitting can be efficiently leveraged by extending standard methodologies to enhance precision and accuracy in the retrieved rotation curves;

• Optical depth profiles can be easily constrained knowing the 2D temperature and the molecular emitting layers only, with reasonable accuracy. These results are comparable to estimates from other independent and more complicated methods.

A first interesting development of this work would be to extract the temperature field directly from the rotation curves instead of assuming it from Law et al. (2021a), by fitting for the temperature structure as imprinted in the vertically separated rotation curves. With this method it would be possible to characterize the 2D temperature field T(r, ɀ) directly from molecular rotation curves without the restriction to optically thick emitters. We leave this analysis to subsequent work.

We also emphasize that the simultaneous characterization of the 2D pressure structure and temperature scalar field has the potential to reveal the 2D gas density structure in disks. A multi-molecule kinematical analysis could then represent a solid further step toward a more comprehensive physical and chemical view of the 2D structure of disks.

Acknowledgements

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01055.L ALMA 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. V.P. and S.F. are 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. S.F. acknowledges financial contribution from PRIN-MUR 2022YP5ACE. CL has been supported by the UK Science and Technology research Council (STFC) via the consolidated grant ST/W000997/1. GL acknowledges the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement No 823823 (DUSTBUSTERS). The authors thank the anonymous referee for valuable comments that improved the scientific value of this work.

Appendix A Flat emission scenario

For every considered line, in Fig. A.1 we show the difference between the extracted rotation curves in case of flat and elevated emission, where the error on each point is evaluated as the sum in quadrature of the errors on the two curves. The errors are generally small, as the uncertainties on the curves are in the first place underestimated, as explained in Sect. 4.1. The largest discrepancies not compatible within the error bars are obtained for ¹²CO, ¹³CO, and HCN, for which the flat assumption is clearly incorrect, as these molecules emit from elevated layers in the disk. This was expected, and these molecules do not represent an issue as their emission has a high S/N and their emitting layer is well-constrained. As expected, we obtain a large discrepancy for ¹²CO, especially at larger radii: in case of flat emission, the curve is overestimated in the outer regions of the disk (R > 200 au), leading in turn to an overestimate of the disk mass, as described in Sect. 4.2.

Fig. A.1 Difference between the rotation curves assuming flat and elevated emission, for each of the considered lines.

Appendix B Sample extension: Beam smearing effect

We show in Fig. B.1 the beam smearing effects for HCN, as previously shown in Fig. 5 for C₂H in Sect. 3.2.3.

Fig. B.1 Beam smearing effect for HCN.

Panel (a): Azimuthally averaged integrated intensity profile of the HCN line. Emission peaks are marked with dashed red lines and the emission gap is highlighted with a dashed grey line. Panel (b): Rotation curve of HCN (purple line) compared to the one of ¹²CO (green line). The rotation velocity of HCN is super-Keplerian just before the gap and sub- Keplerian after. It shows the opposite behavior in correspondence of the emission peaks. Panel (c): Anticorrelation between the HCN integrated intensity gradient and its difference in rotation velocity compared to¹² CO. The two quantities show an opposite modulation and they cross at the radii corresponding to the emission gap or peaks.

Appendix C H₂CO: A complicated case

As mentioned in Sect. 4.1, the rotation curve extraction with eddy does not provide reliable results for a small range of radii around R ~ 150 au, as the signal-to-noise in correspondence of this emission gap is too low to allow an efficient reconstruction of the profile. We note that the extraction in this region is not reliable by looking at the river plots at the corresponding radii and at how efficiently they are aligned. We show in Fig. C.1 the river plots, before and after the alignment with the SHO method, at two different radii: R ~ 240 au and R ~ 160 au. We see that in the first case the alignment procedure works perfectly and the aligned spectra form a straight line, while in the second case (corresponding to the emission gap radius) the method fails in aligning the intensity peaks.

Fig. C.1 River plots before and after the alignment procedure with the SHO method from eddy, atR ~ 240 au (two upper panels) andR ~ 150 au (two bottom panels).

Appendix D Thermally stratified fit: Corner plot

We show in Fig. D.1 the corner plot for the simultaneous multi-molecule fit performed with DySc on our seven extracted rotation curves.

Fig. D.1 Corner plot of the multi-molecule stratified fit performed with DySc on the seven considered rotation curves. The uncertainties reported in the plot represent the default errors we obtain from the MCMC process, not the one we consider for our analysis, as explained in Sect. 4.1.

Appendix E Nonparametric emitting layers

We show in Fig. E.2 the corner plot obtained from the multimolecule fit assuming nonparametric emitting layers for all the molecules, with the dynamical estimates of the stellar mass, disk mass, and scale radius. We also show in Fig. E.1 the same plot as Fig. 10 for the nonparametric case.

Fig. E.1 Normalized difference between the squared rotation velocities of C18O and 12CO, assuming nonparametric emitting layers for all the molecules: comparison between the data (light blue dots), the prediction from isothermal model (green line) and from thermally stratified model (purple line).

Fig. E.2 Corner plot of the multi-molecule stratified fit performed with DySc on the seven considered rotation curves, assuming nonparametric emitting layers for all the molecules.

Appendix F Complete sample: Thermally stratified versus isothermal model

Here we show in Fig. F.1 the normalized difference between the squared rotation velocities of respectively ¹³CO, HCO⁺, H₂CO, HCN, and C₂H with respect to the ¹²CO. We do not discuss in detail H₂CO, since the rotation curve extraction proved to be quite difficult and we do not have the necessary accuracy to draw a reasonable conclusion. We also do not recover any significant difference between the two considered models for what concerns HCO⁺.

On the contrary, we find a significant improvement in the fit for ¹³CO, when considering thermal structure: in this case, as for C¹⁸O, the stratified model better reproduces the measured velocity difference. This is the case also for HCN and C₂H, even if the difference between the two models is less appreciable here because in the considered range of radii their emitting layers are closer to the ¹²CO with respect to the ¹³CO. Moreover, for these two molecules we can clearly see a global oscillation in the observations that can be traced back to the distortion of the rotation curves due to the beam smearing effect, as explained in Sect. 3.2.3. Finally, for ¹³CO we find small-scale pressure modulated oscillations similar to what we see for C¹⁸O.

We observe that, especially in the ¹³CO case, for some radial ranges (130–210 au, >330 au) even the stratified model struggles to reproduce the data. This can be traced back to the inflexibility of the model in this case, caused by the not complete suitability of the temperature parameterization from Dullemond et al. (2020) we are assuming, as the thermal structure might be more complex.

Fig. F.1 Normalized difference between the squared rotation velocities of each molecule and 12CO: comparison between the data (light blue dots), the prediction from isothermal model (green line) and from thermally stratified model (purple line).

References

All Tables

All Figures

	Fig. 1 Optical depth radial profiles for H2CO, HCO+, HCN, and C2H found in this work, taking into account the molecular emitting layers. Emission peaks and emission gaps, when present, are highlighted in red and black dashed lines, respectively. We also report the τ = 1 line for easier visualization.
In the text

	Fig. 2 Teardrop plots for C2H (left) and HCN (right). Each row of the plot represents an averaged spectrum at separate annuli after shifting and stacking. The color scale indicates the intensity of the emission line. In both cases, hyperfine transitions (two for C2H and three for HCN) are clearly visible.
In the text

	Fig. 3 Improvement in the fitting of the spectra using multi-Gaussian methods. Top panel: comparison between the single-Gaussian (light blue) and the double-Gaussian (pink) methods used to fit the intensity peaks of C2H at 50 au. Bottom panel: same as top panel, but for the three hyperfine components of the HCN transition.
In the text

	Fig. 4 Comparison between the rotation curves for C2H (top) and HCN (bottom) using the double- or triple-Gaussian methods (pink) and the single-Gaussian method already implemented in eddy (dashed light blue) for the fitting of the spectra.
In the text

Fig. 5 Beam smearing effect for C2H. Panel a: azimuthally averaged integrated intensity profile of the C2H line. Emission peaks are marked with dashed red lines and the emission gap is highlighted with a dashed grey line. Panel b: rotation curve of C2H (purple line) compared to that of 12CO (green line). The rotation velocity of C2H is super-Keplerian just before the gap and sub-Keplerian after. It shows the opposite behavior in correspondence of the emission peaks. Panel c: anticorrelation between the C2H integrated intensity gradient and its difference in rotation velocity compared to 12CO. The two quantities show an opposite modulation, and they cross at the radii corresponding to the emission gap or peaks.

In the text

Fig. 6 Sketch of the different effects that beam smearing due to line intensity gradients (pink dashed line) and additional pressure gradient effects due to the presence of gas density substructures (light blue solid line) have on a standard rotation curve (black solid line). We assume that intensity gradients trace gas density gradients; a peak or gap in line intensity corresponds to a peak or gap in gas density. Line intensity and gas density peaks are marked with dashed red lines, and line intensity and gas density gaps are highlighted with a dashed grey line. At the bottom we show a colored map of the radial intensity. The solid black line represents the rotation curve of a molecule showing no gas density or line intensity substructures. The typical curve distortion due to beam smearing effect is shown by the pink dashed line, and is in agreement with our observations. The expected distortion of the curve due to the presence of gas density substructures is shown by the light blue solid line. The rings and gaps in the gas density induce additional pressure gradients that influence the rotation curve, according to the $v_{\varphi }^2(R,ɀ)$ expression shown in the upper right corner. These pressure gradients are positive before rings and beyond gaps locations, while they are negative before gaps and beyond rings locations. Thus, in this case, the deviations in rotational velocity are expected to be in the opposite direction of the beam smearing-induced variations.

In the text

	Fig. 7 Multi-molecule fit of the seven considered rotation curves with the thermally stratified model. For each tracer, the red dots with error bars represent the data, while the blue lines are the best-fit models of each velocity profile.
In the text

	Fig. 8 Measured ${\chi }_{\text{red }}^2$ (normalized to the maximum value) of the multimolecule fits. We changed the number of molecules considered in the fit, adding the lines one by one.
In the text

	Fig. 9 Retrieved 2D pressure structure P(R, ɀ) as a function of radius and height across the disk. The location of the selected molecular emitting layers are marked with colored solid lines, and highlight that we are tracing the pressure along an extended radial and vertical range in the disk.
In the text

	Fig. 10 Normalized difference between the squared rotation velocities of C18O and 12CO: comparison between the data (light blue dots), the prediction from isothermal model (green line), and the prediction from thermally stratified model (purple line).
In the text

	Fig. A.1 Difference between the rotation curves assuming flat and elevated emission, for each of the considered lines.
In the text

	Fig. B.1 Beam smearing effect for HCN.
In the text

	Fig. C.1 River plots before and after the alignment procedure with the SHO method from eddy, atR ~ 240 au (two upper panels) andR ~ 150 au (two bottom panels).
In the text

	Fig. D.1 Corner plot of the multi-molecule stratified fit performed with DySc on the seven considered rotation curves. The uncertainties reported in the plot represent the default errors we obtain from the MCMC process, not the one we consider for our analysis, as explained in Sect. 4.1.
In the text

	Fig. E.1 Normalized difference between the squared rotation velocities of C18O and 12CO, assuming nonparametric emitting layers for all the molecules: comparison between the data (light blue dots), the prediction from isothermal model (green line) and from thermally stratified model (purple line).
In the text

	Fig. E.2 Corner plot of the multi-molecule stratified fit performed with DySc on the seven considered rotation curves, assuming nonparametric emitting layers for all the molecules.
In the text

	Fig. F.1 Normalized difference between the squared rotation velocities of each molecule and 12CO: comparison between the data (light blue dots), the prediction from isothermal model (green line) and from thermally stratified model (purple line).
In the text