Rotation curves in protoplanetary disks with thermal stratification. Physical model and observational evidence in MAPS disks

P. Martire, C. Longarini, G. Lodato, G. P. Rosotti, A. Winter, S. Facchini, C. Hardiman, M. Benisty, J. Stadler, A. F. Izquierdo, et al.

► To cite this version:

P. Martire, C. Longarini, G. Lodato, G. P. Rosotti, A. Winter, et al. Rotation curves in protoplanetary disks with thermal stratification. Physical model and observational evidence in MAPS disks. Astronomy & Astrophysics - A&A, 2024, 686, 10.1051/0004-6361/202348546. insu-04836867

HAL Id: insu-04836867

https://insu.hal.science/insu-04836867v1

Submitted on 13 Dec 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Rotation curves in protoplanetary disks with thermal stratification

## Physical model and observational evidence in MAPS disks

P. Martire<sup>1,2</sup>①, C. Longarini<sup>1,8</sup>, G. Lodato<sup>1</sup>, G. P. Rosotti<sup>1</sup>①, A. Winter<sup>3</sup>, S. Facchini<sup>1</sup>①, C. Hardiman<sup>4</sup>①, M. Benisty<sup>3,5</sup>, J. Stadler<sup>3,5</sup>①, A. F. Izquierdo<sup>2</sup>①, and Leonardo Testi<sup>6,7</sup>①

Received 10 November 2023 / Accepted 16 February 2024

## ABSTRACT

**Context.** In recent years, the gas kinematics probed by molecular lines detected with ALMA has opened a new window into the of study protoplanetary disks. High spatial and spectral resolution observations have revealed the complexity of protoplanetary disk structure. Drawing accurate interpretations of these data allows us to better comprehend planet formation.

**Aims.** We investigate the impact of thermal stratification on the azimuthal velocity of protoplanetary disks. High-resolution gas observations reveal velocity differences between CO isotopologues, which cannot be adequately explained with vertically isothermal models. The aim of this work is to determine whether a stratified model can explain this discrepancy.

**Methods.** We analytically solved the hydrostatic equilibrium for a stratified disk and we derived the azimuthal velocity. We tested the model with SPH numerical simulations and then we used it to fit for the star mass, disk mass, and scale radius of the sources in the MAPS sample. In particular, we used <sup>12</sup>CO and <sup>13</sup>CO datacubes.

**Results.** When thermal stratification is taken into account, it is possible to reconcile most of the inconsistencies between rotation curves of different isotopologues. A more accurate description of the CO rotation curves offers a deeper understanding of the disk structure. The best-fit values of star mass, disk mass, and scale radius become more realistic and more in line with previous studies. In particular, the quality of the scale radius estimate significantly increases when adopting a stratified model. In light of our results, we computed the gas-to-dust ratio and the Toomre Q parameter. Within our hypothesis, for all the sources, the gas-to-dust ratio appears higher but still close to the standard value of 100 (within a factor of 2). The Toomre Q parameter suggests that the disks are gravitationally stable (Q > 1). However, the systems that show spirals presence are closer to the conditions of gravitational instability (Q ~ 5).

**Key words.** accretion, accretion disks – hydrodynamics – molecular data – protoplanetary disks

## 1. Introduction

Our understanding of the physical properties of protoplanetary disks has improved in recent years thanks to the Atacama Large Millimeter Array (ALMA; ALMA Partnership et al. 2015). High spectral and spatial resolution gas observations enable us to probe density, temperature, and velocity fields of protostellar disks, gaining unique information about their structure (Law et al. 2021; Calahan et al. 2021; Teague et al. 2022; Miotello et al. 2023; Pinte et al. 2023; Lodato et al. 2023). More recently, the large program Molecules with ALMA at Planet-forming Scales (MAPS; Öberg et al. 2021) targeted five protoplanetary disks (MWC 480, IM Lup, GM Aur, HD 163296, and AS 209) in several molecular lines. For optically thick line emission, the gas temperature can be measured along the emission surface directly

from the peak surface brightness of the channel maps (Law et al. 2021). Given the varying heights of these emitting layers surfaces, it is possible to infer the thermal structure in disks, proving the existence of a vertical thermal stratification in them (Dartois et al. 2003; Rosenfeld et al. 2013; Pinte et al. 2018), as expected from basic radiative transfer arguments (Chiang & Goldreich 1997; D’Alessio et al. 1998, 1999). Although disk models have usually been considered as vertically isothermal, the vertical gradient of temperature leads to considerable corrections in the calculation of density structure and azimuthal velocity, which results in several percent deviations from the Keplerian velocity (Rosenfeld et al. 2013). Accounting for such differences is important not only to infer stellar masses, but also to accurately constrain the disk pressure structure and disk mass. Such parameters are of great importance with respect

<footnote n="1">¹ Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, Milano, 20133, Italy</footnote>

<footnote n="2">² Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
e-mail: martire@strw.leidenuniv.nl</footnote>

<footnote n="3">³ Laboratoire Lagrange, Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, 06304 Nice, France</footnote>

<footnote n="4">⁴ School of Physics and Astronomy, Monash University, Clayton, VIC 3800, Australia</footnote>

<footnote n="5">⁵ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="6">⁶ Alma Mater Studiorum Università di Bologna, Dipartimento di Fisica e Astronomia (DIFA), Via Gobetti 93/2, 40129 Bologna, Italy</footnote>

<footnote n="7">⁷ INAF – Osservatorio Astrofisico di Arcetri, L.go E. Fermi 5, 50125 Firenze, Italy</footnote>

<footnote n="8">⁸ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
e-mail: cl2000@cam.ac.uk</footnote>
to interpreting velocity deviations, which may serve as sign-
posts for planets (Pinte et al. 2018; Rabago & Zhu 2021; Bollati
et al. 2021; Izquierdo et al. 2021; Bae et al. 2021), dust trap-
ping (Teague et al. 2018; Rosotti et al. 2020) or disk instabilities
(Hall et al. 2020; Terry et al. 2022; Longarini et al. 2021;
Barraza-Alfaro et al. 2021).

In this paper, we analytically derive the density and velocity
field of protostellar disks with thermal stratification, generaliz-
ing the work of Takeuchi & Lin (2002). A similar analysis of
MAPS data with vertically isothermal disks was performed by
Lodato et al. (2023). We test the model against hydrodynami-
cal simulations and we apply it to the whole MAPS sample for
¹²CO and ¹³CO data. In Sect. 2, we present the model, solving
the vertical hydrostatic equilibrium and obtaining an expression
for the azimuthal velocity. In Sect. 3, we present the numerical
setup and a comparison between the model and simulations. In
Sect. 4, we apply the model and we discuss our findings. Finally,
in Sect. 5, we compute the gas-to-dust ratio and the Toomre Q
parameter, and we draw our conclusions.

## 2. Model

### 2.1. Assumptions

In our analytical calculations, we did not make any assumption
on the surface density Σ, considering it as arbitrary. However,
to apply the model to observations, we were forced to choose
a parameterization for the surface density and we assume that
it is described by the self-similar solution from Lynden-Bell &
Pringle (1974):

$$
\Sigma = \frac{(2-\gamma)M_d}{2\pi R_c^2} \left(\frac{R}{R_c}\right)^{-\gamma} \exp\left[-\left(\frac{R}{R_c}\right)^{2-\gamma}\right], \quad (1)
$$

where $M_d$ and $R_c$ are the disk mass and the scale radius respec-
tively; $R$ is the cylindrical radius and $\gamma$ is a free parameter
describing the steepness of the surface density. The disk density
at the midplane $\rho_{\text{mid}}$ is:

$$
\rho_{\text{mid}} = \frac{\Sigma}{\sqrt{2\pi}H_{\text{mid}}} \propto R^{-(\gamma+(3-q)/2)} \exp\left[-\left(\frac{R}{R_c}\right)^{2-\gamma}\right], \quad (2)
$$

where $H_{\text{mid}} = c_{\text{s,mid}}/\Omega_k$ is the typical scale height of the disk at the midplane, $c_{\text{s,mid}} = \sqrt{k_b T_{\text{mid}}/(\mu m_p)} \propto R^{-q/2}$ is the sound speed at the disk midplane, $k_b$ is the Boltzmann constant, $T_{\text{mid}} = T_{\text{mid,100}}(R/100 \text{ au})^{-q}$ is the temperature at midplane, $\mu$ is the mean molecular weight (usually assumed to be 2.1), $m_p$ is the proton mass, and $\Omega_k = \sqrt{GM_*/R^3}$ is the Keplerian angular velocity (G is the gravitational constant and $M_*$ is the stellar mass).

From the literature (Chiang & Goldreich 1997; Dullemond
et al. 2020) and observational data (Rosenfeld et al. 2013; Pinte
et al. 2018; Law et al. 2021), we know that protoplanetary disks
are thermally stratified. We take this into account by defining a
function, *f*, that describes the dependency of the temperature, *T*,
on height, such that:

$$
T(R, z) = T_{\text{mid}}(R)f(R, z), \quad (3)
$$

$$
c_s^2(R, z) = c_{\text{s,mid}}^2(R) f(R, z). \quad (4)
$$

We underline that the isothermal case can be obtained consider-
ing *f* ≡ 1, thus *T* = *T*<sub>mid</sub>(*R*). As for the density, we assume that:

$$
\rho = \rho(R, z) = \rho_{\text{mid}}(R)g(R, z), \quad (5)
$$

where *g* describes how the density changes vertically. We note
that in order to smoothly connect the functions above to their
value at midplane it is necessary that *f*(*z* = 0) = 1 = *g*(*z* = 0).
Assuming a barotropic fluid, the pressure, *P*, is given by:

$$
P(R, z) = P_{\text{mid}}(R)fg(R, z) = c_{\text{s,mid}}^2(R)\rho_{\text{mid}}(R)fg(R, z). \quad (6)
$$

While the profile of *f* is arbitrary, this does not hold for *g*, whose
value is set by solving the vertical hydrostatic equilibrium.

## 2.2. Hydrostatic equilibrium and rotation curve

To compute the vertical density profile we assume a non-self-
gravitating disk under the condition of hydrostatic equilibrium
in the vertical direction:

$$
\frac{1}{\rho} \frac{dP}{dz} = - \frac{d\Phi_{\star}}{dz}, \qquad (7)
$$

where $\Phi_{\star} = -GM_{\star}/r$ is the stellar potential ($r = \sqrt{R^2 + z^2}$ is the spherical radius). Equation (7) can be written as (for further details, see Appendix A):

$$
c_{s, \text{mid}}^2 f \frac{d \log(fg)}{d z} = - \Omega_k^2 z \left[ 1 + \left( \frac{z}{R} \right)^2 \right]^{-3/2} . \quad (8)
$$

Solving for $\log(fg)$, we find:

$$
\log(fg) = - \frac{1}{H_{\text{mid}}^2} \int_0^z \frac{z'}{f} \left[ 1 + \left( \frac{z'}{R} \right)^2 \right]^{-3/2} dz', \quad (9)
$$

and hence the density is given by:

$$
\rho(R, z) = \frac{\rho_{\text{mid}}(R)}{f(R, z)} \exp \left\{ - \frac{1}{H_{\text{mid}}^2} \int_0^z \frac{z'}{f(R, z')} \left[ 1 + \left( \frac{z'}{R} \right)^2 \right]^{-3/2} dz' \right\}. \quad (10)
$$

Assuming the condition of centrifugal balance, the rotation curve
is given by the radial component of Navier-Stokes equation:

$$
v_{\phi}^{2}(R, z) = \frac{R}{\rho} \frac{dP}{dR}(R, z) + R \frac{d\Phi_{\star}}{dR}(R, z). \quad (11)
$$

The first term in Eq. ((11)) can be written as (for further details, see Appendix A)

$$
\frac{R}{\rho} \frac{dP}{dR} = c_{s, \text{mid}}^2 f \left[ \frac{d \log P_{\text{mid}}}{d \log R} + \frac{d \log(fg)}{d \log R} \right], \quad (12)
$$

and the second one as

$$
R \frac{d\Phi_{\star}}{dR}(R, z) = v_{\star}^2 \left[ 1 + \left( \frac{z}{R} \right)^2 \right]^{-3/2}, \quad (13)
$$

where $v_k = \sqrt{GM_*/R}$ is the Keplerian velocity. Therefore, the
rotation curve is

$$
v_{\phi}^2(R, z) = v_k^2 \left\{ \left[ 1 + \left( \frac{z}{R} \right)^2 \right]^{-3/2} + \left[ \frac{d \log P_{\text{mid}}}{d \log R} + \frac{d \log(fg)}{d \log R} \right] \left( \frac{H}{R} \right)_{\text{mid}}^2 f(R, z) \right\}. \quad (14)
$$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparison between the analytical models and the simulation after eight outer orbits and at z = 20 au. Left panel: density field. Right panel: Pressure gradient term of the rotation curve. The red dots represent the simulation data, while the blue and orange line show the thermally stratified model and the isothermal one, respectively. The model with thermal stratification matches very well the simulation.</figcaption></figure>

In the self-similar case, this becomes:

$$v_{\phi}^{2}(R, z) = v_{k}^{2} \left\{ \left[ 1 + \left( \frac{z}{R} \right)^{2} \right]^{-3/2} - \left[ \gamma' + (2 - \gamma) \left( \frac{R}{R_{c}} \right)^{2-\gamma} - \frac{d \log(fg)}{d \log R} \right] \left( \frac{H}{R} \right)_{\text{mid}}^{2} f(R, z) \right\}, \quad (15)$$

where $\gamma' = \gamma + (3 + q)/2$. Each term of Eq. (15) can be easily interpreted: $[1 + (z/R)^2]^{-3/2}$ is the star contribution at the height z, $\gamma'$ is the effect of the power law scaling of the pressure, $(2 - \gamma)(R/R_c)^{2-\gamma}$ is the effect of the exponential truncation and the logarithmic term is the effect of the vertical stratification. Since the latter is the derivative of a product, we do not know its sign a priori; thus, we also do not know if the rotation is accelerated or slowed down by thermal stratification (see Appendix A). In any case, in all our attempts this term never dominates over the variation of gravity with z. Thus, we found rotation to slow down with z and this effect is more pronounced as compared to the isothermal case when considering the parameters of the MAPS sample.

We underline that for the isothermal case ($f \equiv 1$), this expression reduces to the one derived and analyzed by Lodato et al. (2023), while Eq. (9) is simplified as:

$$\log g = -\frac{1}{H_{\text{mid}}^2} \int_0^z \frac{z' dz'}{[1 + (z'/R)^2]^{3/2}} = -\frac{R^2}{H_{\text{mid}}^2} \left( 1 - \frac{1}{\sqrt{1 + z^2/R^2}} \right). \quad (16)$$

Therefore, the density in the isothermal case is given by:

$$\rho(R, z) = \rho_{\text{mid}}(R) \exp \left[ \frac{R^2}{H_{\text{mid}}^2} \left( \frac{1}{\sqrt{1 + z^2/R^2}} - 1 \right) \right]. \quad (17)$$

If the disk is self-gravitating, we should add to the right-hand side of Eq. (11) the self-gravitating term (Bertin & Lodato 1999):

$$\frac{R}{\rho} \frac{d\Phi_d}{dR}(R, z) = G \int_0^\infty \left[ K(k) - \frac{1}{4} \left( \frac{k^2}{1 - k^2} \right) \times \left( \frac{r}{R} - \frac{R}{r} + \frac{z^2}{Rr} \right) E(k) \right] \sqrt{\frac{r}{R}} k \Sigma(r) dr, \quad (18)$$

where $K(k)$ and $E(k)$ are complete elliptic integrals (Abramowitz & Stegun 1970) and $k^2 = 4Rr/[(R+r)^2 + z^2]$.

## 2.3. Temperature prescriptions

The two parameterizations of the vertical temperature more often used are given by Dartois et al. (2003) and Dullemond et al. (2020). In this work we will use the one by Dullemond et al. (2020), which is given by:

$$T(R, z)^4 = T_{\epsilon}^4(R) + \frac{1}{2} T_{\text{atm}}^4(R) \left[ 1 + \tanh\left(\frac{z}{Z_q(R)} - \alpha\right) \right] \quad (19)$$

and, thus,

$$f(R, z) = \left\{ \left( \frac{T_{\epsilon}}{T_{\text{mid}}} \right)^4 + \frac{1}{2} \left( \frac{T_{\text{atm}}}{T_{\text{mid}}} \right)^4 (R) \left[ 1 + \tanh \left( \frac{z}{Z_q(R)} - \alpha \right) \right] \right\}^{1/4}, \quad (20)$$

where the atmospheric temperature is parameterized as $T_{\text{atm}}(R) = T_{\text{atm},100}(R/100 \text{ au})^{-q_{\text{atm}}}$, $T_{\epsilon}$ is considered as an approximation of the temperature at the midplane $T_{\epsilon} \approx T_{\text{mid}}$. $Z_q(R)$ is defined as $Z_q(R) = \zeta_{100}(R/100 \text{ au})^{\beta}$ and $\alpha$ is a parameter that describes where the transition from midplane to atmospheric temperature occurs in the vertical direction. We note that in this case, $f(R, z = 0) \neq 1$ and, thus, the temperature does not smoothly connect to its value at midplane. We discuss this in Appendix B, but we underline that Eq. (19) is a good approximation for the five disks within the MAPS large program in most of the radial extent of the disk.
Once the function *f* is defined, Eqs. (10) and (14) can be solved semi-analytically and entirely specify the rotation curve. We have implemented this calculation in DYSC<fnref n="1" />.

### 3. Comparison with numerical simulations

In this work, we performed numerical smoothed particle hydrodynamics (SPH) simulations of protostellar disks using the PHANTOM code (Price et al. 2018). This code is widely used in the astrophysical community to study gas and dust dynamics in accretion disks (Dipierro et al. 2015; Ragusa et al. 2017; Curone et al. 2022) and it has recently been employed in kinematical studies (Pinte et al. 2018; Hall et al. 2020; Terry et al. 2022; Verrios et al. 2022). The aim of this simulation is to test the model before applying it to actual data.

To test the analytical model, we simulated a thermally stratified disk using the parameters of MWC 480 presented in Law et al. (2021). The simulation has been performed with $N = 10^6$ gas particles, initially distributed as a tapered power law density profile, smoothed at the inner radius, with $\gamma = 1$ and $R_c = 150$ au, between $R_{in} = 10$ au and $R_{out} = 400$ au. The mass of the star is 2.1 $M_\odot$. For the temperature structure we used the Dullemond prescription given by Eq. (19), with $\zeta_0 = 7$ au, $\alpha = 2.78, \beta = -0.05, T_{mid,100} = 27$ K, $q = 0.23, T_{atm,100} = 69$ K, and $q_{atm} = 0.7$. The $\alpha_{SS}$ Shakura & Sunyaev (Shakura & Sunyaev 1973) viscosity coefficient has been set to 0.005. No self-gravity or dust have been included in the simulation.

We let the system evolve and reach hydrostatic equilibrium. We observed that after a couple of orbits the system reaches a relaxed state. We decided to analyze the output of the simulation after eight outer orbits (~45 kyr). In Fig. 1, we show a comparison between the density and the velocity of the simulations (red dots) at $z = 20$ au and both the isothermal and stratified model predictions. The red dots represent the azimuthal average of the respective quantity computed by averaging over all SPH particles within each of the 50 radial bins and the error bar is the corresponding standard deviation. Since we are plotting quantities at $z = 20$ au, we have excluded the inner points because at those radii ($R < 100$ au) the disk has a smaller hydrostatic height $H$, causing numerical issues in the resolution of our simulation. The stratified model perfectly describes the density and the velocity field of the simulation and is a significant improvement over the isothermal one. In particular, in the right panel of Fig. 1, we can see that the difference between the azimuthal velocity and the Keplerian velocity ($v_k^2 - v_\phi^2)/v_k^2$ reaches the 10–12% and only the stratified model is able to reproduce it.

### 4. Applying the model

#### 4.1. The curves

In this section, we describe how we applied the model to the entire sample of disks from the MAPS ALMA Large Program (Öberg et al. 2021). We performed our fits under the assumption of vertically isothermal or stratified disk to then compare the results. For the vertically isothermal model, the thermal structure is defined by the hydrostatic height of the disk at $R = 100$ au and the power law coefficient of the temperature profile, $q$. These parameters were taken from Zhang et al. (2021). As for the stratified model, Law et al. (2021) obtained the two-dimensional (2D) temperature structure of the MAPS disks, using the Dullemond et al. (2020) prescription given by Eq. (19). We note that the

<table><caption>**Table 1.** Velocity extraction method, orientation parameters, thermal parameters and emitting surfaces for ¹²CO and ¹³CO data of the MAPS disks.</caption><thead><tr><th></th><th>MWC 480</th><th>IM Lup</th><th>GM Aur</th><th>HD 163296</th><th>AS 209</th></tr></thead><tbody><tr><td>Extraction</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>¹²CO</td><td>Gauss</td><td>Dbell</td><td>Dbell</td><td>Dbell</td><td>Gauss</td></tr><tr><td>¹³CO</td><td>Gauss</td><td>Dbell</td><td>Dbell</td><td>Dbell</td><td>Gauss</td></tr><tr><td>Orientation</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>i (deg)</td><td>37.00</td><td>47.50</td><td>53.20</td><td>46.69</td><td>35.00</td></tr><tr><td>PA (deg)</td><td>328.15</td><td>144.50</td><td>53.98</td><td>312.75</td><td>85.20</td></tr><tr><td>Isothermal</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>H<sub>100</sub> (au)</td><td>10</td><td>10</td><td>7.5</td><td>8.4</td><td>6</td></tr><tr><td>q</td><td>0.82</td><td>0.66</td><td>0.3</td><td>0.84</td><td>0.5</td></tr><tr><td>Stratified</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>T<sub>mid</sub> (K)</td><td>27</td><td>25</td><td>20</td><td>24</td><td>25</td></tr><tr><td>T<sub>atm</sub> (K)</td><td>69</td><td>36</td><td>48</td><td>63</td><td>37</td></tr><tr><td>q</td><td>0.23</td><td>0.02</td><td>0.01</td><td>0.18</td><td>0.18</td></tr><tr><td>q<sub>atm</sub></td><td>0.7</td><td>-0.03</td><td>0.55</td><td>0.61</td><td>0.59</td></tr><tr><td>ζ<sub>0</sub> (au)</td><td>7</td><td>3</td><td>13</td><td>9</td><td>5</td></tr><tr><td>α</td><td>2.78</td><td>4.91</td><td>2.57</td><td>3.01</td><td>3.31</td></tr><tr><td>β</td><td>-0.05</td><td>2.07</td><td>0.54</td><td>0.42</td><td>0.02</td></tr><tr><td>¹²CO Surface</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>z<sub>0</sub> (au)</td><td>17.04</td><td>34.13</td><td>32.00</td><td>27.14</td><td>16.47</td></tr><tr><td>ψ</td><td>1.35</td><td>0.99</td><td>0.97</td><td>1.07</td><td>1.24</td></tr><tr><td>R<sub>t</sub> (au)</td><td>579.43</td><td>889.40</td><td>729.91</td><td>534.00</td><td>327.52</td></tr><tr><td>q<sub>t</sub></td><td>1.63</td><td>3.18</td><td>3.22</td><td>2.99</td><td>3.01</td></tr><tr><td>¹³CO Surface</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>z<sub>0</sub> (au)</td><td>11.52</td><td>22.84</td><td>18.21</td><td>16.09</td><td>4.13</td></tr><tr><td>ψ</td><td>1.09</td><td>1.27</td><td>1.14</td><td>1.12</td><td>0.96</td></tr><tr><td>R<sub>t</sub> (au)</td><td>402.77</td><td>529.06</td><td>512.13</td><td>392.75</td><td>180.22</td></tr><tr><td>q<sub>t</sub></td><td>1.87</td><td>1.65</td><td>2.73</td><td>3.43</td><td>3.59</td></tr></tbody></table>

**Notes.** The orientation parameters and the emitting surfaces are taken from Izquierdo et al. (2023), the thermal parameters for the isothermal model are taken from Zhang et al. (2021), for the stratified model from Law et al. (2021).

rotation curve traced by a specific molecule is defined by:

$$v_{\text{rot}}^2(R) = v_{\phi}^2(R, z(R)), \qquad (21)$$

where $z(R)$ is the height of the emitting layer of the considered molecule. For the emitting layer, we use:

$$z(R) = z_0 \left( \frac{R}{100 \text{ au}} \right)^{\psi} \exp \left[ - \left( \frac{R}{R_t} \right)^{q_t} \right], \qquad (22)$$

where the best-fit parameters have been obtained by Izquierdo et al. (2023). All the parameters used are summarized in Table 1.

The rotation curves (Fig. 2) can be obtained through different moment maps, according to the disk emission. We underline that with respect to the rotation curve extraction, we are only interested in measuring velocities from the frontside. Since three of the sources have strong contributions from the disk backside, we used a double-Bell decomposition to distinguish between these two components as introduced (Izquierdo et al. 2022). In this work, we have used an improved algorithm that performs this decomposition based on velocity priors obtained from the DISCMINER models (Izquierdo et al., in prep.).

#### 4.2. Results

The ¹²CO and ¹³CO data were simultaneously fitted with both the isothermal and stratified model including the self-gravitating

<footnote n="1">¹ The code is publicly available at https://github.com/crislong/DySc</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. $^{12}$CO and $^{13}$CO rotation curves of the MAPS disks extracted with DISCMINER.</figcaption></figure>

contribution. The results are shown in Figs. 3–7, and the best-fitting parameters reported in Table 2 and Appendix C.

To quantify the importance of thermal stratification, we computed the relative difference between the squares of $^{12}$CO and $^{13}$CO rotation curves, as shown in Fig. 8. According to the vertical isothermal model, this quantity is:

$$ (v_{13}^2 - v_{12}^2)_{\text{iso}} = v_k^2 q \frac{\sqrt{1 + z_{12}^2/R^2} - \sqrt{1 + z_{13}^2/R^2}}{\sqrt{(1 + z_{13}^2/R^2)(1 + z_{12}^2/R^2)}}, \quad (23) $$

which solely depends on the different height of the tracer, since it is assumed that the temperature does not change vertically. As for the stratified model, the expression is more complex, since it involves the evaluation of the term given by Eq. (9) at different heights. In this case, we expect to observe larger differences between the velocity of the two isotopologues, since there is an additional shift caused by the different emission temperature. To determine the importance of vertical stratification, we quantified the maximum value of the velocity shift between $^{12}$CO and $^{13}$CO, which can be predicted in the isothermal case:

$$ \frac{(v_{13}^2 - v_{12}^2)_{\text{iso}}}{v_k^2} \approx q \frac{\Delta z^2}{2R^2} < 5\%, \quad (24) $$

where we assumed that typically $z/R < 0.5$. Hence, if the quantity $(v_{13}^2 - v_{12}^2)/v_k^2$ is higher than 5%, the system cannot be

described by an isothermal model, while it is likely that vertical stratification plays a significant role. It is important to note that the Eq. (24) depends on the star mass through $v_k^2$. We normalized the squared differences of the velocities by the square of the velocity for $^{13}$CO since this quantity is independent of the stellar mass. Figure 8 shows this quantity for the studied systems.

In the following, we will present the results of each disk, along with a discussion of the importance of thermal stratification. To compare the results, we performed our fits for both the vertically isothermal and stratified case. In addition, we computed the dust mass from millimetric emission at 283 GHz, using (Hildebrand 1983):

$$ M_{\text{dust}} = \frac{d^2 F_v}{\kappa_v B_v(T)}, \quad (25) $$

where d is the distance, $F_v$ is the flux density in Jy, $\kappa_v = 2.3(\nu/230 \text{ GHz})^{0.4} \text{cm}^2 \text{g}^{-1}$ is the dust opacity, and $B_v$ is the blackbody spectrum. In our analysis, we assumed $T = 20 \text{ K}$ and $\nu = 283 \text{ GHz}$, while the flux densities were extracted from MAPS data. We recall that this equation implies that dust emission is optically thin. The results are reported in Table 3.

### 4.3. MWC 480

MWC 480 is a ~7 Myr Herbig Ae star located in the Taurus-Aurigae star-forming region at a distance of d = 162 pc
<figure><img src="image_4.png" /><figcaption>Fig. 3. Rotation curves of MWC 480 obtained from CO isotopologues. Left panel: Rotation curve obtained from ¹²CO data (red points), along with our best-fitting curve for the stratified model (blue line) and the isothermal one (orange line). Right panel: same for the ¹³CO data.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Same as Fig. 3 but for IM Lup.</figcaption></figure>

(Montesinos et al. 2009). The most recent value of the stellar mass has been derived dynamically by Izquierdo et al. (2023) to be $M_{\star} = 1.97 M_{\odot}$. Zhang et al. (2021) through 2D thermo-chemical models computed disk mass and scale radius of the MAPS disks. For MWC 480, these values are $M_{\mathrm{d}} = 0.16 M_{\odot}$ and $R_{\mathrm{c}} = 200$ au.

By inspecting the ¹²CO and ¹³CO rotation curves (Fig. 2), there is no evident sign of thermal stratification, since the two curves do not differ significantly. Figure 3 shows that the two models are nearly indistinguishable, but in Fig. 8 we see that the stratified model is better at reproducing the data. When we assume an isothermal model, we obtain $M_{\star} = 1.969 \pm 0.002 M_{\odot}$, $M_{\mathrm{d}} = 0.201 \pm 0.002 M_{\odot}$ and $R_{\mathrm{c}} = 80 \pm 1$ au, while for the stratified model, it is $M_{\star} = 2.027 \pm 0.002 M_{\odot}$, $M_{\mathrm{d}} = 0.150 \pm 0.002 M_{\odot}$.

and $R_{\mathrm{c}} = 128 \pm 1$ au. The disk mass obtained with the stratified model is in agreement with the literature value (Zhang et al. 2021). Since the reduced chi-squared $\chi^2_{\mathrm{red}}$ is smaller in the stratified case (see Table 2), we adopted it as the best-fit model.

### 4.4. IM Lup

IM Lup is a young pre-main sequence star (~1 Myr) located in the Lupus star-forming region at a distance of 158pc (Gaia Collaboration 2018). The dynamical stellar mass is estimated to be 1.1 $M_{\odot}$ (Teague et al. 2021), and it hosts an unusually large disk, extending out to ≈ 300 au in the dust continuum and out to ≈1000 au in the gas (Cleeves et al. 2016). The dust continuum emission shows clear evidence of a spiral morphology, which
<figure><img src="image_6.png" /><figcaption>Fig. 5. Same as Fig. 3 but for GM Aur.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Same as Fig. 3, but for HD 163296.</figcaption></figure>

may be triggered by gravitational instability (Huang et al. 2018). Cleeves et al. (2016) first estimated the disk mass from mm visibilities and found a massive disk of $0.2 M_{\odot}$. Verrios et al. (2022) claimed that the spiral structure of IM Lup could be generated by an embedded protoplanet. They performed numerical SPH simulations of planet-disk interaction and then post-processed them to compare their results with CO, dust, and scattered light emission. Interestingly, a high disk mass (~0.1 $M_{\odot}$) is required to match the scattered light image, so that the sub-micron sized grains could remain well coupled in the top layers of the disk. Cleeves et al. (2016) first estimated the disk scale radius $R_c = 100$ au by comparing SED to a simple tapered power-law density profile. Afterwards, Pinte et al. (2018) analyzed CO data and found that a tapered power law density profile with $R_c = 284$ au more optimally reproduces the data. They also analyzed the rotation

curve of the disk and found that while the inner disk is in good agreement with Keplerian rotation around a $1 \pm 0.1 M_{\odot}$ star, both the $^{12}$CO and the $^{13}$CO rotation curves become sub-Keplerian in the outer disk. The authors attributed this effect to the pressure gradient. Lodato et al. (2023) analyzed $^{12}$CO and $^{13}$CO rotation curves and fitted for star mass, disk mass and scale radius with an isothermal model. In particular, the authors found that for the rotation curves extracted with EDDY the best-fit are $M_{\star} = 1.012 \pm 0.003 M_{\odot}, M_{\mathrm{d}} = 0.096 \pm 0.003 M_{\odot}, R_{\mathrm{c}} = 89 \pm 1$ au and DISMINER are $M_{\star} = 1.02 \pm 0.02, M_{\mathrm{d}} = 0.10 \pm 0.01 M_{\odot}, R_{\mathrm{c}} = 66 \pm 1$ au. We underline that in this work the rotation curves have been obtained again, and they are different from those of Lodato et al. (2023). This is also true for the case of GM Aur.

Figure 4 shows both the isothermal and stratified fit. While for $^{12}$CO, both models describe the rotation curve well, for $^{13}$CO
<figure><img src="image_8.png" /><figcaption>Fig. 7. Same as Fig. 3, but for AS 209.</figcaption></figure>

<table><caption>**Table 2.** Results of the fitting procedure and reduced chi-squared for the two different models: isothermal and stratified.</caption><thead><tr><th></th><th>M<sub>★</sub> (M<sub>☉</sub>)</th><th>M<sub>d</sub> (M<sub>☉</sub>)</th><th>R<sub>c</sub> (au)</th><th>χ<sup>2</sup><sub>red</sub></th></tr></thead><tbody><tr><td>MWC 480</td><td></td><td></td><td></td><td></td></tr><tr><td>Isothermal</td><td>1.969 ± 0.002</td><td>0.201 ± 0.002</td><td>80 ± 1</td><td>11.21</td></tr><tr><td>Stratified</td><td>2.027 ± 0.002</td><td>0.150 ± 0.002</td><td>128 ± 1</td><td>6.14</td></tr><tr><td>IM Lup</td><td></td><td></td><td></td><td></td></tr><tr><td>Isothermal</td><td>1.055 ± 0.002</td><td>0.200 ± 0.003</td><td>55 ± 1</td><td>35.68</td></tr><tr><td>Stratified</td><td>1.194 ± 0.002</td><td>0.106 ± 0.002</td><td>115 ± 1</td><td>6.29</td></tr><tr><td>GM Aur</td><td></td><td></td><td></td><td></td></tr><tr><td>Isothermal</td><td>0.872 ± 0.003</td><td>0.312 ± 0.003</td><td>56 ± 1</td><td>90.84</td></tr><tr><td>Stratified</td><td>1.128 ± 0.002</td><td>0.118 ± 0.002</td><td>96 ± 1</td><td>8.48</td></tr><tr><td>HD 163296</td><td></td><td></td><td></td><td></td></tr><tr><td>Isothermal</td><td>1.842 ± 0.002</td><td>0.124 ± 0.001</td><td>38 ± 1</td><td>29.60</td></tr><tr><td>Stratified</td><td>1.948 ± 0.002</td><td>0.134 ± 0.001</td><td>91 ± 1</td><td>19.74</td></tr><tr><td>AS 209</td><td></td><td></td><td></td><td></td></tr><tr><td>Isothermal</td><td>1.272 ± 0.003</td><td>0.042 ± 0.003</td><td>45 ± 1</td><td>25.13</td></tr><tr><td>Stratified</td><td>1.311 ± 0.001</td><td>0.0002 ± 0.0002</td><td>126 ± 2</td><td>10.55</td></tr></tbody></table>

the isothermal model fails, since the velocity shift is so high that it cannot be explained just in terms of emitting surface. This difference is clearly visible when considering the $\chi^2_{red}$, which for the stratified model is considerably smaller. The best-fit parameters for the isothermal model are $M_{★} = 1.055 \pm 0.002 M_{☉}$, $M_{d} = 0.200 \pm 0.003 M_{☉}$, and $R_{c} = 55 \pm 1$ au; while for the stratified model are $M_{★} = 1.1994 \pm 0.002 M_{☉}$, $M_{d} = 0.106 \pm 0.002 M_{☉}$, and $R_{c} = 115 \pm 1$ au. The effects of thermal stratification are visible in Fig. 8. At $R \sim 250$ au, the difference in the data between $^{12}$CO and $^{13}$CO is on the order of ~10% and it significantly increases in the outer part. In that case, the stratified model is not able to explain that difference either. Izquierdo et al. (2023) pointed out that the emission from the outer disk is so diffuse that the retrieval of the emitting surface, as well as the velocity extraction, needs to be taken with care. This is possibly an effect of external photoevaporation. Indeed, despite the

very weak external radiation field irradiating IM Lup, Haworth et al. (2017) showed that the disk is sufficiently large that the outer part, which is weakly gravitationally bound, can undergo photoevaporation.

### 4.5. GM Aur

GM Aur is a T-Tauri star in the Taurus-Auriga star-forming region that hosts a transition disk. The stellar mass has been estimated dynamically to be $M_{★} = 1.1 M_{☉}$ by Teague et al. (2021), in agreement with previous measurements (Macías et al. 2018). Its CO morphology is very complex, showing spiral arms, tails, and interactions with the environments (Huang et al. 2021). From thermochemical models of MAPS data, Schwarz et al. (2021) obtained a disk mass of $M_{d} = 0.2 M_{☉}$ and a scale radius of $R_{c} = 111$ au, making GM Aur a possibly gravitationally unstable disk. Lodato et al. (2023) performed a fitting procedure for the star mass, disk mass, and scale radius, using an isothermal model, finding that for GM Aur, the two CO lines provide inconsistent rotation curves, which cannot be attributed only to a difference in the height of the emitting layer. In addition, the authors provided a simple order-of-magnitude estimate of the expected velocity shift due to thermal stratification, concluding that the difference between the two rotation curves could not be explained by this effect. They drew this conclusion by taking into account the different temperature of the two molecules at their emission height, $z_{i}(R)$, given by Law et al. (2021). However, as shown in Appendix A, in the azimuthal velocity, it is not only essential to know the temperature at $(R, z)$, but also its radial and vertical gradient at that location.

By analyzing the rotation curves of the two CO isotopolgues (Fig. 2), a systematic shift between $^{12}$CO and $^{13}$CO curves is clearly visible, which may possibly be attributed to thermal stratification. When we perform the fitting with the isothermal model, we obtain (as the best-fit parameters) $M_{★} = 0.872 \pm 0.003 M_{☉}$, $M_{d} = 0.312 \pm 0.003 M_{☉}$, and $R_{c} = 56 \pm 1$ au. This is in agreement with Lodato et al. (2023) and leads to a high $\chi^2_{red}$ (see Table 2). As a matter of fact, Fig. 5 shows that an
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Relative difference between the squares of observed ¹²CO and ¹³CO rotation curves predicted by the thermally stratified model (blue line), the isothermal model (orange line), and the data (red dots). Except for AS 209, where this quantity is negative in the inner part, it is clearly visible that data are well reproduced by the stratified model. Indeed, the difference of speed between the two curves cannot be explained just in terms of different height.</figcaption></figure>
isothermal model is not able to reproduce both ¹²CO and ¹³CO rotation curves. Conversely, when thermal stratification is taken into account, the two rotation curves are compatible and are in agreement with data, especially for $R > 180$ au. In this case, the best-fit value for the star mass is $M_★ = 1.128 ± 0.002 M_⊙$, which is in line with the literature values (Teague et al. 2021; Macías et al. 2018). As for the disk mass, the best-fit value is $M_d = 0.118 ± 0.002 M_⊙$. Finally, the best-fit value for the scale radius is $R_c = 96 ± 1$ au, almost twice the value obtained with the isothermal model and in good agreement with Schwarz et al. (2021). A stratified model reproduces the difference between ¹²CO and ¹³CO rotation curves very well, as shown in Fig. 8, which leads to a significant decrease in the $\chi^2_{red}$ value.

## 4.6. HD 163296

HD 163296 is one of the most well-studied Herbig Ae star system at millimeter wavelengths thanks to its relative close distance ($d = 101$ pc) and bright disk. The disk presents several features that suggest ongoing planet formation, such as dust rings, deviations from Keplerian velocities due to gas pressure variations, “kinks” in the CO emission, and meridional flows (Isella et al. 2016, 2018; Pinte et al. 2018, 2023; Teague et al. 2018; Izquierdo et al. 2022, 2023; Calcino et al. 2022).

This system has also been extensively studied because there is considerable evidence to support a massive disk. By modeling the dust lines, Powell et al. (2019) found that the disk mass is $M_d = 0.21 M_⊙$. As for the scale radius, de Gregorio-Monsalvo et al. (2013) through radiative transfer modeling found that $R_c = 125$ au is the value that better reproduces dust and CO ALMA observations. Guidi et al. (2016) presented a multiwavelength ALMA and VLA study of the disk and, via a modeling of the visibilities, they found that the best-fit value of the scale radius is $R_c = 118$ au, in agreement with de Gregorio-Monsalvo et al. (2013).

When we fit data with a vertical isothermal model, we obtain (as the best-fit parameters) $M_★ = 1.842 ± 0.002 M_⊙$, $M_d = 0.124 ± 0.001 M_⊙$ and $R_c = 38 ± 1$ au. While the star mass is realistic, the scale radius is unrealistically small compared to the gas emission extent on the order of 400 au (Law et al. 2021). Additionally, the isothermal model is not able to reproduce the difference between the rotation curves of the two CO isotopologues (see Fig. 8), resulting in a relatively poor fit with a large $\chi^2_{red}$. If we include the 2D thermal structure, the quality of the fit increases (see $\chi^2_{red}$ in Table 2). In this case, the best-fit for stellar mass and disk mass does not change significantly ($M_★ = 1.948 ± 0.002 M_⊙$, $M_d = 0.134 ± 0.001 M_⊙$), while the scale radius does shift to $R_c = 91 ± 1$ au. Comparing our result for the disk mass to the literature values, we observe that our fit gives a value that is roughly half. Figure 6 shows that both the isothermal and the stratified model aptly describe the rotation curve of ¹²CO and ¹³CO. However, the shift between them, presented in Fig. 8, is well recovered only by the stratified model, which offers only a partial explanation for the significant increase of the plotted quantity. The presence of pressure modulated substructures in the rotation curves (Izquierdo et al. 2023) impacts the quality of the fit (clearly visible in Fig. 8). One possible solution would be to include them in the fitting model.

## 4.7. AS 209

AS 209 is a young T-Tauri star in the Ophiucus star forming region ($d \sim 121$ pc). The most recent stellar mass estimate is

$M_★ = 1.14 M_⊙$ (Izquierdo et al. 2023). Fedele et al. (2018) gave an estimate for the scale radius $R_c = 80$ au through the modeling of the mm visibilities. Afterwards, through thermochemical modeling, they found a dust mass of $M_{dust} = 3.5 × 10^{-4} M_⊙$; with a gas-to-dust ratio of 100, this translates into $M_d = 0.0035 M_⊙$, which is in agreement with the recent value $M_d = 0.0045 M_⊙$ of Zhang et al. (2021). Interestingly, when inspecting the rotation curves of AS 209 (Fig. 2), the ¹³CO is slower compared to the ¹²CO, despite it being closer to the midplane. This trend is observed up to ~125 au. A possible explanation for this is the compactness of the disk, which makes more difficult to extract a precise emitting surface due to beam smearing. Indeed, line centroids from pixels near the center of the disk are an averaged composition of multiple surrounding velocities because of the limited resolution and the steepness of $v(r)$. Since AS209 is the smallest disk in the sample, it is more prone to be affected by this in the largest fraction of its total extent compared to the other sources. When we fit with the isothermal model, we obtained, as the best-fit parameters, $M_★ = 1.272 ± 0.003 M_⊙$, $M_d = 0.042 ± 0.003 M_⊙$, and $R_c = 45 ± 1$ au. When we fit with the stratified model, we obtain as the best-fit parameters $M_★ = 1.311 ± 0.001 M_⊙$, and $R_c = 126 ± 2$ au, while for the disk mass we report a $3-\sigma$ upper limit of $M_d = 0.00025 ± 0.00025 M_⊙$, since the best-fit parameter is compatible with zero. Both models are shown in Fig. 7. As for ¹²CO, the two models behave in the same way, showing little difference in the outer edge. Conversely, for ¹³CO the isothermal model works better in the inner part, where ¹³CO is slower, while in the outer part the stratified model describes well the rotation curve. According to the $\chi^2_{red}$, the stratified model describes the data better (see Table 2).

# 5. Discussion

## 5.1. Thermal stratification in MAPS disks

Table 2 presents a summary of the findings of this study, comparing the isothermal model with the stratified one. It is evident from the results that the reduced $\chi^2$ value consistently decreases when employing the stratified model. This indicates that the inclusion of thermal stratification provides a more effective way of describing the observed data. In this context, MWC 480 is particularly interesting. Despite the small kinematic signatures of thermal stratification, as depicted in Fig. 8, the quality of the stratified fit is higher and it yields more reliable values for star mass, disk mass, and scale radius. On the opposite side, GM Aur is the system that shows the strongest effects of thermal stratification, given that the ¹²CO and ¹³CO systematically shifted over all the radial extent of the disk. The introduction of thermal stratification is able to reconcile these differences, reducing the $\chi^2_{red}$ by an order of magnitude. The only case where the stratified model encounters challenges in accurately describing both curves is in AS 209. This system is peculiar because the compactness of the disk influences the extraction of emission surfaces. Consequently, contrary to what expected, we observe that the ¹³CO rotates slower than the ¹²CO in the inner part. Despite that, the $\chi^2_{red}$ is smaller when thermal stratification is taken into account.

## 5.2. Disk masses

In this paragraph, we aim to contextualize our work within the broader framework of disk mass estimation.

One solid tracer of the disk mass is the carbon dioxide HD, which is a good tracer of the disk gas because it follows the
<table><caption>**Table 3.** Continuum fluxes at 283 GHz, dust masses from Eq. (25) and gas-to-dust ratio using the best-fit value of the disk mass of the stratified model.</caption><thead><tr><th></th><th><i>F</i><sub>283</sub> (mJy)</th><th><i>M</i><sub>dust</sub> (<i>M</i><sub>☉</sub>)</th><th>Gas-to-dust ratio</th></tr></thead><tbody><tr><td>MWC 480</td><td>943.51</td><td>0.00138</td><td>108</td></tr><tr><td>IM Lup</td><td>536.25</td><td>0.00075</td><td>134</td></tr><tr><td>GM Aur</td><td>347.95</td><td>0.00049</td><td>240</td></tr><tr><td>HD 163296</td><td>1127.97</td><td>0.00064</td><td>202</td></tr><tr><td>AS 209</td><td>414.83</td><td>0.00034</td><td>&lt; 192</td></tr></tbody></table>

distribution of molecular hydrogen and its emission is sensitive to the total mass. The first detection of HD emission in a protoplanetary disk comes from Bergin et al. (2013) for TW Hya. Afterwards, the detection of HD *J* = 1–0 line was used to estimate disk mass of GM Aur. The HD based disk mass is 2.5–20.4×10⁻² *M*☉ (McClure et al. 2016), in line with our estimate of 0.118 *M*☉. Finally, the non-detection of HD in HD163296 (Kama et al. 2020) translates into an upper limit for the disk mass of 0.067 *M*☉, which is almost half of the value we obtained in this work.

Another reliable method to trace the disk mass uses the N₂H⁺. This molecule is a chemical tracer of CO-poor gas and can be used to measure the CO–H₂ ratio and calibrate CO-based gas masses. By combining N₂H⁺ with C¹⁸O, Trapman et al. (2022) estimated disk masses of three protoplanetary disks, including GM Aur. The value they obtained, 1.5–9.6 × 10⁻² *M*☉, is slightly higher compared to our estimate, but it is in an overall good agreement. This method has also been used to probe disk masses of protoplanetary disks in the Lupus star-forming region (Anderson et al. 2022).

In this context, it is worth mentioning observations of the ¹³C¹⁷O, a very rare CO isotopologue. Booth et al. (2019) observed this molecule in HD 163296 and this allows for a precise disk mass measurement to be obtained. These authors found that the disk mass that is better at reproducing observations is *M*<sub>d</sub> = 0.31 *M*☉, which is discrepant with our inferred value.

As for the dust, its ability to trace the mass is discussed in the next subsection.

### 5.3. Gas-to-dust ratio

With the knowledge of the disk mass, it is possible to evaluate the gas-to-dust ratio, using Eq. (25) for the dust mass. The results are shown in Table 3. We found values between 100–250, within a factor of 2 from the usually assumed value of 100. This is surprisingly, due to the assumptions we made to obtain the dust mass. Indeed, as we have already mentioned, the optically thin hypothesis for dust emission could lead to a difference of a more than a factor of 2 in the dust mass calculation (Guidi et al. 2016), underestimating it. In addition, the dust opacity could also vary of a factor of ~10, depending on the grain size and composition. Hence, it is significant overall that the inferred gas-to-dust ratio is so close to the standard value. As for AS 209, we estimated an upper limit for this quantity. Indeed, according to Veronesi et al. (in prep.), the minimum measurable mass with the rotation curve is 5% of the star mass. Taking this value as an upper limit for AS 209 disk mass, it is possible to give an upper limit for the gas-to-dust ratio.

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Toomre Q parameter of the MAPS disks with the best-fit parameters of the stratified model. We excluded AS 209 because its best-fit disk mass is compatible with zero.</figcaption></figure>

### 5.4. Toomre Q

To investigate the presence of gravitational instability, we used our best-fit parameters for the stratified model to compute the Toomre parameter (Toomre 1964) which, based on the hypothesis of nearly Keplerian disk ($\kappa \approx \Omega$), is

$$Q \approx \frac{c_s \Omega}{\pi G \Sigma} = 2 \left[ \frac{H}{R} \right]_{\text{mid}} \frac{M_{\star}}{M_{\text{d}}} \left( \frac{R}{R_{\text{c}}} \right)^{-1} \exp \left[ \frac{R}{R_{\text{c}}} \right], \quad (26)$$

where we used Eq. (1) for the surface density. According to the WKB quadratic dispersion relation (Lin & Shu 1964; Toomre 1964), the onset of the instability happens when *Q* ~ 1. Figure 9 shows the profile of the Q parameter for the MAPS sample, except for AS209, since its disk mass estimate is compatible with zero. Every disk is gravitationally stable, according to the Toomre criterion, since *Q* > 1. Interestingly, the two disks that showed spiral structures (IM Lup and GM Aur) have a Toomre profile that is lower than the others, with a minimum value of ~4 for GM Aur and ~6 for IM Lup. Lau & Bertin (1978) showed that a WKB description of gravitational instability can still be obtained under less restrictive conditions compared to the quadratic relation. They showed that disks that are locally stable according to the Q criterion might still generate large scale spiral waves. In general, other mechanisms could increase the critical value of the Toomre parameter, such as external irradiation (Lin & Kratter 2016; Löhnert et al. 2020) or dust-driven gravitational instability (Longarini et al. 2023b,a). Hence, we do not exclude the possibility that gravitational instability is at play in GM Aur and IM Lup.

## 6. Conclusions

The kinematic data of protoplanetary disks display velocity differences between ¹²CO and ¹³CO that cannot be explained through a vertically isothermal model, given the systematic shift between rotation curves of CO isotopologues. In this work, we predict how thermal stratification affects the density and the velocity field of a protoplanetary disk.
We used SPH simulations to test our model, finding excel-
lent agreement, and then we applied it to the MAPS sample.
We extracted the rotation curves of CO isotopologues (¹²CO and
¹³CO) and we carried out fitting for the star mass, disk mass, and
scale radius – both with a vertically isothermal and a stratified
model. The quality of the fit significantly improves when ther-
mal stratification is taken into account and the best-fit parameters
are more realistic and aligned with literature. All our results are
summarized in Table 2.

Typically, when thermal stratification is considered, the best-
fit value for the star mass tends to rise. This can be intuitively
understood, as an isothermal model would favor a star mass that
lies between that of ¹³CO and ¹²CO; its underestimation is due
to the slower rotation of ¹²CO. Conversely, the stratified model
encapsulates the difference between the two curves, mitigating
the underestimation issue and resulting in a more accurate mass
estimate. While an isothermal model provides a satisfactory fit
at small radii, the fit worsens at large radii where the difference
between ¹²CO and ¹³CO is larger. The fit tries to compensate
for this by increasing the disk mass, most of which resides at
large radii, thereby changing the predicted curve only in the outer
parts of the disk. Ultimately, a more accurate description of the
thermal structure through a stratified model leads to a realistic
estimate of the scale radius.

We note that the inclusion of the vertical gradient of tem-
perature into our model results in improved χ² values across all
systems under examination. This work demonstrates the impact
of thermal stratification in the disk dynamics and highlights the
importance of having a precise knowledge of the disk tempera-
ture structure to infer physical quantities such as the stellar mass,
the disk mass and the disk scale radius in a meaningful manner.
Thanks to the generality of our calculations, we can study the
density and velocity profile of thermally stratified disk according
different prescription of temperature simply.

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), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. This work has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement # 823823 (RISE DUSTBUSTERS project). G.R. is funded by the European Union under the European Union's Horizon Europe Research & Innovation Programme No. 101039651 (discEvol) and by the Fondazione Cariplo, grant no. 2022-1217. 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. S.F. acknowledges financial contribution from PRIN-MUR 2022YP5ACE. CH is funded by a Research Training Program Scholarship from the Australian Government, and acknowledges funding from the Australian Research Council via DP220103767. M.B. and J.S. have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (PROTOPLANETS, grant agreement No. 101002188). The authors thank Pietro Curone and Claudia Toci for useful discussions.

References

Abramowitz, M., & Stegun, I. A. 1970, Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables

ALMA Partnership (Brogan, C. L., et al.) 2015, ApJ, 808, L3

Anderson, D. E., Cleeves, L. I., Blake, G. A., et al. 2022, ApJ, 927, 229

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41

Bae, J., Teague, R., Andrews, S. M., et al. 2022, ApJ, 934, L20

Bae, J., Teague, R., & Zhu, Z. 2021, ApJ, 912, 56

Barraza-Alfaro, M., Flock, M., Marino, S., & Pérez, S. 2021, A&A, 653, A113

Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature, 493, 644

Bertin, G., & Lodato, G. 1999, A&A, 350, 694

Bollati, F., Lodato, G., Price, D. J., & Pinte, C. 2021, MNRAS, 504, 5444

Booth, A. S., Walsh, C., Ilee, J. D., et al. 2019, ApJ, 882, L31

Calahan, J. K., Bergin, E. A., Zhang, K., et al. 2021, ApJS, 257, 17

Calcino, J., Hilder, T., Price, D. J., et al. 2022, ApJ, 929, L25

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, ApJ, 832, 110

Curone, P., Izquierdo, A. F., Testi, L., et al. 2022, A&A, 665, A25

D'Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411

D'Alessio, P., Cantó, J., Hartmann, L., Calvet, N., & Lizano, S. 1999, ApJ, 511, 896

Dartois, E., Dutrey, A., & Guilloteau, S. 2003, A&A, 399, 773

de Gregorio-Monsalvo, I., Ménard, F., Dent, W., et al. 2013, A&A, 557, A133

Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73

Dullemond, C. P., Isella, A., Andrews, S. M., Skobleva, I., & Dzyurkevich, N. 2020, A&A, 633, A137

Dutrey, A., Guilloteau, S., Piétu, V., et al. 2008, A&A, 490, L15

Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 605, A16

Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Galloway-Sprietsma, M., Bae, J., Teague, R., et al. 2023, ApJ, 950, 147

Guidi, G., Tazzari, M., Testi, L., et al. 2016, A&A, 588, A112

Hall, C., Dong, R., Teague, R., et al. 2020, ApJ, 904, 148

Haworth, T. J., Facchini, S., Clarke, C. J., & Cleeves, L. I. 2017, MNRAS, 468, L108

Hildebrand, R. H. 1983, QJRAS, 24, 267

Huang, J., Andrews, S. M., Pérez, L. M., et al. 2018, ApJ, 869, L43

Huang, J., Bergin, E. A., Öberg, K. I., et al. 2021, ApJS, 257, 19

Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101

Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L49

Izquierdo, A. F., Testi, L., Facchini, S., Rosotti, G. P., & van Dishoeck, E. F. 2021, A&A, 650, A179

Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022, ApJ, 928, 2

Izquierdo, A. F., Testi, L., Facchini, S., et al. 2023, A&A, 674, A113

Kama, M., Trapman, L., Fedele, D., et al. 2020, A&A, 634, A88

Lau, Y. Y., & Bertin, G. 1978, ApJ, 226, 508

Law, C. J., Teague, R., Loomis, R. A., et al. 2021, ApJS, 257, 4

Lin, M.-K., & Kratter, K. M. 2016, ApJ, 824, 91

Lin, C. C., & Shu, F. H. 1964, ApJ, 140, 646

Lodato, G., Rampinelli, L., Viscardi, E., et al. 2023, MNRAS, 518, 4481

Löhnert, L., Krätschmer, S., & Peeters, A. G. 2020, A&A, 640, A53

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17

Longarini, C., Lodato, G., Toci, C., et al. 2021, ApJ, 920, L41

Longarini, C., Armitage, P. J., Lodato, G., Price, D. J., & Ceppi, S. 2023a, MNRAS, 522, 6217

Longarini, C., Lodato, G., Bertin, G., & Armitage, P. J. 2023b, MNRAS, 519, 2017

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Macías, E., Espaillat, C. C., Ribas, Á., et al. 2018, ApJ, 865, 37

McClure, M. K., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 831, 167

Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, in Astronomical Society of the Pacific Conference Series, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 534, 501

Molyarova, T., Akimkin, V., Semenov, D., et al. 2017, ApJ, 849, 130

Montesinos, B., Eiroa, C., Mora, A., & Merín, B. 2009, A&A, 495, 901

Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Piétu, V., Dutrey, A., Guilloteau, S., Chapillon, E., & Pety, J. 2006, A&A, 460, L43

Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13

Pinte, C., Teague, R., Flaherty, K., et al. 2023, in Astronomical Society of the Pacific Conference Series, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 534, 645

Powell, D., Murray-Clay, R., Pérez, L. M., Schlichting, H. E., & Rosenthal, M. 2019, ApJ, 878, 116

Price, D. J., Wurster, J., Tricco, T. S., et al. 2018, PASA, 35, e031

Rabago, I., & Zhu, Z. 2021, MNRAS, 502, 5325

Ragusa, E., Dipierro, G., Lodato, G., Laibe, G., & Price, D. J. 2017, MNRAS, 464, 1449

Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16
Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173

Schwarz, K. R., Calahan, J. K., Zhang, K., et al. 2021, ApJS, 257, 20

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14

Simon, M., Guilloteau, S., Beck, T. L., et al. 2019, ApJ, 884, 42

Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344

Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12

Teague, R., Bae, J., Aikawa, Y., et al. 2021, ApJS, 257, 18

Teague, R., Bae, J., Andrews, S. M., et al. 2022, ApJ, 936, 163

Terry, J. P., Hall, C., Longarini, C., et al. 2022, MNRAS, 510, 1671

Toci, C., Rosotti, G., Lodato, G., Testi, L., & Trapman, L. 2021, MNRAS, 507, 818

Toomre, A. 1964, ApJ, 139, 1217

Trapman, L., Rosotti, G., Bosman, A. D., Hogerheijde, M. R., & van Dishoeck, E. F. 2020, A&A, 640, A5

Trapman, L., Zhang, K., van't Hoff, M. L. R., Hogerheijde, M. R., & Bergin, E. A. 2022, ApJ, 926, L2

Trapman, L., Rosotti, G., Zhang, K., & Tabone, B. 2023, ApJ, 954, 41

Veronesi, B., Ragusa, E., Lodato, G., et al. 2020, MNRAS, 495, 1913

Verrios, H. J., Price, D. J., Pinte, C., Hilder, T., & Calcino, J. 2022, ApJ, 934, L11

Zhang, K., Booth, A. S., Law, C. J., et al. 2021, ApJS, 257, 5
## Appendix A: Computing the pressure gradient

Assuming a barotropic fluid, the pressure contribution to the rotation curve is

$$ \frac{R}{\rho} \frac{dP}{dR} \Big|_z = \frac{R}{\rho} \left( c_s^2 \frac{d\rho}{dR} \Big|_z + \rho \frac{dc_s^2}{dR} \Big|_z \right), \quad (A.1) $$

which can also be expressed as:

$$ \begin{aligned} \frac{R}{\rho} \frac{dP}{dR} \Big|_z &= c_s^2 \frac{R}{\rho_{\text{mid}} g} \left( \rho_{\text{mid}} \frac{dg}{dR} \Big|_z + g \frac{d\rho_{\text{mid}}}{dR} \Big|_z \right) + \\ &R \left( c_{s, \text{mid}}^2 \frac{df}{dR} \Big|_z + f \frac{dc_{s, \text{mid}}^2}{dR} \Big|_z \right) = \\ &c_s^2 \left[ -\gamma' - (2-\gamma) \left( \frac{R}{R_c} \right)^{2-\gamma} + R \frac{d \log(fg)}{dR} \Big|_z \right]. \end{aligned} \quad (A.2) $$

Assuming hydrostatic equilibrium in the vertical direction, we can write an explicit expression for the last term in Eq. (A.2):

$$ \begin{aligned} -\Omega_k^2 z \left[ 1 + \left(\frac{z}{R}\right)^2 \right]^{-3/2} &= -\frac{d\Phi_*}{dz} = \frac{1}{\rho} \frac{dP}{dz} = \\ \frac{1}{\rho} P_{\text{mid}} \frac{d(fg)}{dz} &= c_{s, \text{mid}}^2 f \frac{d \log(fg)}{dz} \end{aligned} \quad (A.3) $$

and, assuming $f(z=0) = 1$ (therefore $T(R, z=0) = T_{\text{mid}}(R)$),

$$ \log(fg) = -\frac{1}{H_{\text{mid}}^2} \int_0^z \frac{z'}{f} \left[ 1 + \left(\frac{z'}{R}\right)^2 \right]^{-3/2} dz'. \quad (A.4) $$

Deriving with respect to R we obtain:

$$ \begin{aligned} \frac{R}{\rho} \frac{d \log(fg)}{dR} \Big|_z = R & \left[ \frac{2}{H_{\text{mid}}^3} \frac{dH_{\text{mid}}}{dR} \int_0^z \frac{1}{f} \frac{z' dz'}{(1 + (z'/R)^2)^{3/2}} + \right. \\ & \left. \frac{1}{H_{\text{mid}}^2} \int_0^z \frac{1}{f^2} \frac{df}{dR} \frac{z' dz'}{(1 + (z'/R)^2)^{3/2}} - \frac{1}{H_{\text{mid}}^2} \int_0^z \frac{1}{f} \frac{3z'^3/R^3 dz'}{(1 + (z'/R)^2)^{5/2}} \right] = \\ & - \frac{q}{H_{\text{mid}}^2} \int_0^z \frac{1}{f} \frac{z' dz'}{(1 + z'^2/R^2)^{3/2}} + \frac{R}{H_{\text{mid}}^2} \int_0^z \frac{\dot{f}}{f^2} \frac{z' dz'}{(1 + z'^2/R^2)^{3/2}} + \\ & \frac{3}{H_{\text{mid}}^2} \int_0^z \frac{1}{f} \frac{z' dz'}{(1 + z'^2/R^2)^{5/2}}, \end{aligned} \quad (A.5) $$

where $\cdot = d/dR$. Expanding the logarithmic term is useful since it is now clear that the azimuthal velocity depends both on the temperature ($f$) and on the temperature gradient ($\dot{f}$) along the radial direction. Moreover, it is easy to see the contribution of the vertical thermal stratification: in our model $f$ and $\dot{f}$ are functions of $z$ and thus they have to be integrated to compute azimuthal velocity; whereas in vertically isothermal models, this does not happen, since $f=1$ and $\dot{f}=0$. We note that it is not analytically possible to determine generally the sign of Eq. (A.5) and, thus, the effect that we see on the pressure gradient. However, for all the considered cases it never overcomes the gravitational contribution and it does lead to a faster deceleration in the rotation.

<table><caption>Table B.1. Values of parameters used in the estimate (B.1).</caption><thead><tr><th></th><th>(T<sub>atm</sub>/T<sub>mid</sub>)<sub>100</sub></th><th>q<sub>atm</sub> - q</th><th>α</th><th>|Δ(R = 100 au)|</th></tr></thead><tbody><tr><td>IM Lup</td><td>1.44</td><td>0.05</td><td>4.91</td><td>5.8 · 10<sup>-5</sup></td></tr><tr><td>GM Aur</td><td>2.4</td><td>0.54</td><td>2.57</td><td>4.8 · 10<sup>-2</sup></td></tr><tr><td>AS 209</td><td>1.48</td><td>0.41</td><td>3.31</td><td>1.6 · 10<sup>-3</sup></td></tr><tr><td>HD 163296</td><td>2.63</td><td>0.43</td><td>3.01</td><td>2.9 · 10<sup>-2</sup></td></tr><tr><td>MWC 480</td><td>2.56</td><td>0.47</td><td>2.78</td><td>4.1 · 10<sup>-2</sup></td></tr></tbody></table>

## Appendix B: Dullemond prescription

According to Law et al. (2021), the temperature prescription given by Dullemond et al. (2020) (Eq.(19)) fits the data well, but we note that in this case, $f(R, z=0) \neq 1$ and the temperature does not smoothly connect to its value at midplane since $T(R, z=0) \neq T_{\text{mid}}(R)$. To evaluate this discrepancy, considering $T_{\epsilon}(R) \approx T_{\text{mid}}(R)$, we compute:

$$ \begin{aligned} \Delta(R) = 1 - \frac{T(R, 0)}{T_{\text{mid}}} &= 1 - f(R, z=0) \approx \\ & - \frac{1}{8} \left( \frac{T_{\text{atm}}}{T_{\text{mid}}} \right)^4 (R)(1 - \tanh \alpha) = \\ & - \frac{1}{8} \left( \frac{T_{\text{atm},100}}{T_{\text{mid},100}} \right)^4 \left( \frac{R}{100\text{au}} \right)^{-4(q_{\text{atm}}-q)} (1 - \tanh \alpha), \end{aligned} \quad (B.1) $$

which is strongly dependent on $T_{\text{mid}}$, $T_{\text{atm}}$ and $\alpha$. As shown in Table B.1, the deviation of the actual midplane temperature from $T(R, z=0)$ as computed from the prescription given by Dullemond et al. (2020) is $\lesssim 10\%$ in our regions of interest ($R > 100$ au). This discrepancy could have relevance for systems such as GM Aur, HD 163296 and MWC 480. To examine this further, we conducted fits for these systems using both $T_{\text{mid}}(R)$ and $T(R, z=0)$ as the midplane temperature. It was observed that this choice does not significantly alter the results. Thus, $T_{\text{mid}} \approx T(R, z=0)$ and $f(R, z=0) \approx 1$. We can consider Eqs. (19) and (20) as a good parameterization for the temperature.

## Appendix C: Corner plots

In Figs C.1 and C.2, we present the corner plots of the MCMC fitting procedure for the studied parameters under, respectively, the vertically isothermal and the stratified model.
<figure><img src="image_11.png" /><figcaption>**Fig. C.1.** Corner plots of the MCMC fitting procedure according to the vertically isothermal model. They show the distribution of the three relevant fitting parameters for the five disks of the MAPS large program. From top left to bottom: MWC 480, IM Lup, GM Aur, HD 163296, and AS 209.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. C.2.** Corner plots of the MCMC fitting procedure according to the stratified model. They show the distribution of the three relevant fitting parameters for the five disks of the MAPS large program. From top left to bottom: MWC 480, IM Lup, GM Aur, HD 163296, and AS 209.</figcaption></figure>