Molecules with ALMA at Planet-forming Scales
(MAPS). XII. Inferring the C/O and S/H Ratios in
Protoplanetary Disks with Sulfur Molecules

Romane Le Gal, Karin I. Öberg, Richard Teague, Ryan A. Loomis, Charles J.

Law, Catherine Walsh, Edwin A. Bergin, François Ménard, David J. Wilner,
Sean M. Andrews, et al.

► To cite this version:

Romane Le Gal, Karin I. Öberg, Richard Teague, Ryan A. Loomis, Charles J. Law, et al. Molecules with ALMA at Planet-forming Scales (MAPS). XII. Inferring the C/O and S/H Ratios in Protoplanetary Disks with Sulfur Molecules. The Astrophysical Journal Supplement Series, 2021, 257, 10.3847/1538-4365/ac2583. insu-03672368

HAL Id: insu-03672368

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

Submitted on 29 Jul 2025

**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
Molecules with ALMA at Planet-forming Scales (MAPS) XII:
Inferring the C/O and S/H ratios in Protoplanetary Disks with Sulfur Molecules

ROMANE LE GAL,¹,²,³,⁴,* KARIN I. ÖBERG,¹ RICHARD TEAGUE,¹ RYAN A. LOOMIS,⁵ CHARLES J. LAW,¹

CATHERINE WALSH,⁶ EDWIN A. BERGIN,⁷ FRANÇOIS MÉNARD,³ DAVID J. WILNER,¹ SEAN M. ANDREWS,¹

YURI AIKAWA,⁸ ALICE S. BOOTH,⁶,⁹ GIANNI CATALDI,¹⁰,⁸ JENNIFER B. BERGNER,¹¹,† ARTHUR D. BOSMAN,⁷

L. ILSE CLEEVES,¹² IAN CZEKALA,¹³,¹⁴,¹⁵,¹⁶,¹⁷,† KENJI FURUYA,¹⁰ VIVIANA V. GUZMÁN,¹⁸ JANE HUANG,¹,⁷

JOHN D. ILEE,⁶ HIDEKO NOMURA,¹⁰ CHUNHUA QI,¹ KAMBER R. SCHWARZ,¹⁹,† TAKASHI TSUKAGOSHI,¹⁰

YOSHIHIDE YAMATO,⁸ AND KE ZHANG,²⁰,⁷,‡

¹Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

²IRAP, Université de Toulouse, CNRS, CNES, UT3, F-31000 Toulouse, France

³Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁴IRAM, 300 rue de la Piscine, F-38406 Saint-Martin d'Hères, France

⁵National Radio Astronomy Observatory, Charlottesville, VA 22903, USA

⁶School of Physics and Astronomy, University of Leeds, Leeds, UK, LS2 9JT

⁷Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, USA

⁸Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

⁹Leiden Observatory, Leiden University, 2300 RA Leiden, the Netherlands

¹⁰National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

¹¹University of Chicago Department of the Geophysical Sciences, Chicago, IL 60637, USA

¹²University of Virginia, Charlottesville, VA 22903, USA

¹³Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

¹⁴Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

¹⁵Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

¹⁶Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA

¹⁷Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA

¹⁸Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile

¹⁹Lunar and Planetary Laboratory, University of Arizona, 1629 E. University Blvd, Tucson, AZ 85721, USA

²⁰Department of Astronomy, University of Wisconsin-Madison, 475 N Charter St, Madison, WI 53706, USA

(Received February 23, 2021; Revised July 13, 2021)

ABSTRACT

Sulfur-bearing molecules play an important role in prebiotic chemistry and planet habitability. They are also proposed probes of chemical ages, elemental C/O ratio, and grain chemistry processing. Commonly detected in diverse astrophysical objects, including the Solar System, their distribution and chemistry remain, however, largely unknown in planet-forming disks. We present CS (2 – 1) observations at ~0.′3 resolution performed within the ALMA-MAPS Large Program toward the five disks around IM Lup, GM Aur, AS 209, HD 163296, and MWC 480. CS is detected in all five disks, displaying a variety of radial intensity profiles and spatial distributions across the sample, including intriguing apparent azimuthal asymmetries. Transitions of C₂S and SO were also serendipitously covered but only upper limits are found. For MWC 480, we present complementary ALMA observations at ~0.′5, of CS, ¹³CS, C³⁴S, H₂CS, OCS, and SO₂. We find a column density ratio N(H₂CS)/N(CS) ~ 2/3, suggesting that a substantial part of the sulfur reservoir in disks is in organic form (i.e., CₓHᵧS_z). Using astrochemical disk modeling tuned to MWC 480, we demonstrate that N(CS)/N(SO) is a promising
probe for the elemental C/O ratio. The comparison with the observations provides a super-solar C/O. We also find a depleted gas-phase S/H ratio, suggesting either that part of the sulfur reservoir is locked in solid phase or that it remains in an unidentified gas-phase reservoir. This paper is part of the MAPS special issue of the Astrophysical Journal Supplement.

# 1. INTRODUCTION

Protoplanetary disks are a pivotal stage in the evolution from interstellar molecular clouds to planetary systems. Their chemical structures encode information both on the chemical evolution during star and planet formation, and on the future composition of planets. It is thus of fundamental importance to constrain and understand the chemistry of the principal chemical elements constituting these disks. During the past decade, a myriad of studies focused on oxygen, carbon, and nitrogen chemistry in protoplanetary disks (e.g., Öberg et al. 2011; Guilloteau et al. 2016; Kastner et al. 2018; Cleeves et al. 2018; Bergner et al. 2018, 2019; Pontoppidan et al. 2019), while very little is known about sulfur chemistry in disks. This is probably because, more generally, the chemistry of sulfur in the Universe has remained a long-standing mystery for the past two decades (Ruffle et al. 1999; Kama et al. 2019; Navarro-Almaida et al. 2020).

Sulfur plays an important role in prebiotic chemistry (Chen et al. 2015) and planet habitability (Ranjan et al. 2018; Ruf et al. 2019). It is also one of the most abundant elements in the Universe with a solar value of S/H ~ 1.5 × 10⁻⁵ (Asplund et al. 2009). In the diffuse interstellar medium (ISM) and photon-dissociation regions (PDR) the total amount of sulfur is close to the solar value (Goicoechea et al. 2006; Howk et al. 2006), while in dense molecular gas it is strongly depleted: less than ~1% of the sulfur solar abundance is observed in the gas phase (Tieftrunk et al. 1994; Wakelam et al. 2004; Vastel et al. 2018). Therefore, a question yet to be answered is: what causes the observed sulfur depletion from diffuse to dense gas? While most of the sulfur is suspected to be locked into icy grain mantles (e.g., Millar & Herbst 1990; Ruffle et al. 1999; Vidal et al. 2017; Laas & Caselli 2019), only ~4% of the solar S-abundance has been detected in interstellar ices so far (Palumbo et al. 1997; Boogert et al. 2015). Therefore, the identity of the sulfur reservoir(s) in the ISM remains an open question.

In the Solar System, sulfur-bearing species are routinely detected, in the remnants of our own planet-forming disk such as comets, meteorites, and on planets and their satellites (e.g., Calmonte et al. 2016;

Hirschmann 2016; Lellouch et al. 2007; Franz et al. 2019). In particular, in comets, a dozen of S-bearing species have now been detected (Meier & A'Hearn 1997; Bockelée-Morvan et al. 2004; Biver et al. 2016; Calmonte et al. 2016), including both complex S-molecules (CH₃SH and C₂H₆S<fnref n="1" />) and multi-sulfuretted molecules, such as S₂, CS₂, S₃, and S₄ which have not been detected yet, nor in the ISM, neither in protoplanetary disks. Studying the S-chemistry in disks is therefore crucial to understand the chemical origins of our own Solar System, and more generally, the role of sulfur in astrochemistry.

Among the approximately thirty different molecules detected in disks so far only five<fnref n="2" /> contain sulfur. These include CS, SO, H₂S, H₂CS and SO₂ with the former two detected during the past two decades, and the latter three detected within the past few years due to the significant sensitivity improvements made in radio-interferometry instruments. CS is the most readily detected S-bearing species in proto-planetary disks (Dutrey et al. 1997; Fuente et al. 2010; Guilloteau et al. 2016; Teague et al. 2018b; Le Gal et al. 2019b). SO was the sole oxygenated, sulfur-bearing species detected in disks until the recent detection of SO₂ in one disk (Booth et al. 2021a), which is probably indicative of a general highly reduced or O-poor gas chemistry in most disks. Another interesting point is that, so far, SO has only been detected toward a few young disks with signs of active accretion (Fuente et al. 2010; Guilloteau et al. 2013, 2016; Pacheco-Vázquez et al. 2016; Booth et al. 2018; Rivière-Marichalar et al. 2020). H₂S has long been thought to be a main sulfur reservoir and is a major sulfur carrier in comet 67P/C-G (Calmonte et al. 2016). However, it was only detected recently in the massive disk (~0.15 M☉) GG Tau with a H₂S/CS gas-phase column density ratio of ~1/20 (Phuong et al. 2018), after unsuccessful searches in a handful of other disks (GO Tau, MWC 480, DM Tau, and LkCa 15, Dutrey et al. 2011). While additional H₂S observations in disks are

<footnote>* CNES Fellowship Program Fellow</footnote>

<footnote>† NASA Hubble Fellowship Program Sagan Fellow</footnote>

<footnote>‡ NASA Hubble Fellow</footnote>

<footnote n="1">¹ Note that C₂H₆S has two isomers, ethanethiol (CH₃CH₂SH, also known as ethyl mercaptan and only detected in Orion KL so far, Kolesniková et al. 2014) and dimethyl sulphide ((CH₃)₂S, that is, to our knowledge, not yet detected elsewhere in Space). However, the 67P/C-G measurements did not allow to distinguish the ratio of these two isomers.</footnote>

<footnote n="2">² For this inventory we did not include isotopologues, but note that the isotopologues ¹³CS and C³⁴S are also detected in disks (Le Gal et al. 2019b; Loomis et al. 2020).</footnote>
required to draw firm conclusions, this result casts doubt
on the importance of H₂S in disk gas-phase S-chemistry
and has revived interest in the quest to identify the sul-
fur reservoir in disks (e.g., Kama et al. 2019). The re-
cent detection of H₂CS in the MWC 480 disk, with a
H₂CS/CS gas-phase column density ratio of ~ 1/3, is
in tension with recent models (Le Gal et al. 2019b) and
suggests an incomplete theoretical understanding of disk
S-chemistry. Thus, a better understanding of the S-
chemistry is needed to inform astrochemical models and
constrain the unseen reservoirs of S-bearing species, such
as those locked onto icy dust grains, in disks.

Disks are vertically stratified into atmospheres, warm
molecular layers, and cold midplanes, which are analogs
to PDR, lukewarm molecular clouds, and cold dense
cores, respectively (Aikawa et al. 2002; Bergin et al.
2007; Dutrey et al. 2014). Recent sulfur-bearing species
observations in each of these three types of astrophys-
ical environments – i.e., in a PDR (Fuente et al. 2017;
Rivière-Marichalar et al. 2019), in a protostellar enve-
lope (Drozdovskaya et al. 2018) and in dense cores (Vas-
tel et al. 2018; Navarro-Almaida et al. 2020) – have re-
vived interest in the global quest for understanding the
cycle of sulfur chemistry from molecular cloud to exo-
planetary systems, and are timely for disk S-chemistry
exploration. In this context, Le Gal et al. (2019b) scru-
tinized the reactions pertinent to the sulfur chemistry
within current gas-grain astrochemical models to con-
strain those molecules expected to be particularly abun-
dant in disks, and predicted their radial and vertical
distributions.

Here we present new observations of sulfur-bearing
species in disks taken with the Atacama Large
Millimeter/submillimeter Array (ALMA), as part of
the Molecules with ALMA at Planet-forming Scales
(MAPS) Large Program (Oberg et al. 2021). The ¹²CS
J = 2 – 1 rotational transition was observed toward the
five disks targeted within MAPS, i.e. the disks orbiting
IM Lup, GM Aur, AS 209, HD 163296, and MWC 480
(for which stellar and disk properties are described in
Table 1). The SO JN = 23 – 12 and JN = 54 – 4 and
C2S JN = 87–76 and JN = 1514–1415 rotational transi-
tions were also serendipitously covered. In addition, we
also present new complementary Cycle 6 ALMA obser-
vations (program 2018.1.01631.S, PI: R. Le Gal) toward
the MWC 480 disk, where ¹²CS J = 5 – 4 and its ¹³CS
and C³⁴S isotopologues were observed as well as several
rotational transitions of H₂CS, OCS, and SO₂.

The outline of the paper is as follows: we describe
the observations in Section 2, and we present the re-
sults, including the derivation of column densities and
excitation temperatures, in Section 3. In Section 4,

we present grids of disk chemistry models tuned to the
MWC 480 disk where we obtained the most observa-
tional constraints. In Section 5, we discuss the observa-
tional and modeling results, and summarize our conclu-
sions in Section 6.

## 2. OBSERVATIONS

We used three sets of ALMA observational data.
First, new observations obtained with MAPS (program
number: 2018.1.01055.L, PI: K. Öberg) of the ¹²CS
2 – 1 rotational transition, and of two rotational tran-
sitions of SO and C₂S that were serendipitously cov-
ered in the same program. Second, new observations
obtained with another ALMA program (program num-
ber: 2018.1.01631.S, PI: R. Le Gal) of CS, ¹³CS, and
C³⁴S 5 – 4, four H₂CS rotational transitions, two OCS
rotational transitions, and three SO₂ rotational tran-
sitions. Finally, to get better estimates of the column
densities and excitation temperatures of H₂CS, ¹²CS, and
C³⁴S, we also used already published complementary
ALMA data of additional detected rotational transitions
of these molecules (Le Gal et al. 2019b). The new obser-
vations are further described below and their molecular
transitions, their frequencies, and spectroscopic param-
eters are listed in Table 2.

### 2.1. MAPS observations

The CS 2 − 1 rotational transition was observed in
the five MAPS disks in Band 3 with an angular resolu-
tion of ~ 0.′3 (see Table 2) and a spectral resolution of
71 kHz, corresponding to ~ 0.22 km/s at 97 GHz. More
details about the observations can be found in Oberg
et al. (2021). For the descriptions of the reduction and
imaging procedures applied to the CS 2 − 1 MAPS ob-
servations, we refer the reader to Oberg et al. (2021) and
Czekala et al. (2021), respectively. Here we used the CS
2 − 1 images created with a robustness parameter of 0.5
for the Briggs weighting which results in slightly higher
resolution images than the fiducial images presented in
Oberg et al. (2021) and Law et al. (2021) which used
circularized 0.′3 beams.

While SO and C₂S were not part of the main targeted
molecules within the MAPS program, two of their rota-
tional transitions – namely the 2₃–1₂ (at 99.29987 GHz)
and 5₄ – 4₄ (at 100.2964 GHz) transitions for SO, and
the 8₇ – 7₆ (at 99.86652 GHz) and 15₁₄ – 14₁₅ (at
234.81596 GHz) lines for C₂S – were covered in Band
3 and 6, at lower spectral resolution (1.129 MHz, i.e.
~ 3.4 km/s at 100 GHz and ~ 1.4 km/s at 235 GHz).
After continuum subtraction with the CASA uvcontsub
function, we CLEANed (Högbom 1974) the C₂S and SO
data using the same procedure as outlined in Czekala
<table><caption>**Table 1.** Stellar and Disk Properties as presented in Oberg et al. (2021)</caption><thead><tr><th>Source</th><th>Spectral Type</th><th>dist. [pc]</th><th>incl [°]</th><th>PA [°]</th><th>T<sub>eff</sub> [K]</th><th>L<sub>*</sub> [L<sub>☉</sub>]</th><th>Age<sup>a</sup> [Myr]</th><th>M<sub>*</sub><sup>b</sup> [M<sub>☉</sub>]</th><th>log<sub>10</sub>(<i>dot</i>) [M<sub>☉</sub> yr<sup>-1</sup>]</th><th><i>v</i><sub>sys</sub> [km s<sup>-1</sup>]</th><th>References</th></tr></thead><tbody><tr><td>IM Lup</td><td>K5</td><td>158</td><td>47.5</td><td>144.5</td><td>4266</td><td>2.57</td><td>0.2 – 1.3</td><td>1.1</td><td>-7.9</td><td>4.5</td><td>1,2,3,4,5,6</td></tr><tr><td>GM Aur</td><td>K6</td><td>159</td><td>53.2</td><td>57.2</td><td>4350</td><td>1.2</td><td>~ 2.5</td><td>1.1</td><td>-8.1</td><td>5.6</td><td>1,7,8,9,10,11,12</td></tr><tr><td>AS 209</td><td>K5</td><td>121</td><td>35.0</td><td>85.8</td><td>4266</td><td>1.41</td><td>~ 1</td><td>1.2</td><td>-7.3</td><td>4.6</td><td>1,2,6,13,14</td></tr><tr><td>HD 162396</td><td>A1</td><td>101</td><td>46.7</td><td>133.3</td><td>9332</td><td>17.0</td><td>> 5</td><td>2.0</td><td>-7.4</td><td>5.8</td><td>1,2,6,15,16</td></tr><tr><td>MWC 480</td><td>A5</td><td>162</td><td>37</td><td>148</td><td>8250</td><td>21.9</td><td>~ 7</td><td>2.1</td><td>-6.9</td><td>5.1</td><td>1,17,18,19,20,21</td></tr></tbody></table>

<sup>a</sup>The stellar ages are uncertain by at least a factor of two and should only be considered as preliminary estimates.

<sup>b</sup>All stellar masses have been dynamically determined as described in Teague et al. (2021).

NOTE—References are 1. Gaia Collaboration et al. (2018); 2. Huang et al. (2018); 3. Alcalá et al. (2017); 4. Pinte et al. (2018); 5. Mawet et al. (2012); 6. Andrews et al. (2018); 7. Huang et al. (2020); 8. Macías et al. (2018); 9. Espaillat et al. (2010); 10. Kraus & Hillenbrand (2009); 11. Beck & Bary (2019); 12. Ingleby et al. (2015); 13. Salyk et al. (2013); 14. Huang et al. (2017); 15. Fairlamb et al. (2015); 16. Teague et al. (2019); 17. Liu et al. (2019); 18. Montesinos et al. (2009); 19. Simon et al. (2019); 20. Piétu et al. (2007); 21. Mendigutía et al. (2013)

et al. (2021). As these lines are expected to be weak, we applied a robustness parameter of 1 and 1" uv-taper to improve their imaging and signal-to-noise ratio (SNR). All MAPS images used here are available for download through the ALMA Archive via https://almascience.nrao.edu/alma-data/lp/maps. An interactive browser for this repository is also available on the MAPS project homepage at http://www.alma-maps.info.

## 2.2. Complementary ALMA observations of MWC 480

Independently from the MAPS program, Cycle 6 ALMA observations toward the MWC 480 disk were conducted on 2019 April 30 in three execution blocks (EB) with an angular resolution of ~ 0.55, as part of program 2018.1.01631.S (PI: R. Le Gal). We are presenting and using these data here to complement the data-set of sulfur-bearing molecules observed for this disk. The measurements used ALMA Band 6 receivers, with correlated data divided into thirteen spectral windows (SPWs). SPWs were centered on twelve different rotational transitions of sulfur-bearing molecules, including the CS, $^{13}$CS, and C$^{34}$S 5 – 4 lines, four H$_2$CS lines, two OCS lines, and three SO$_2$ lines. Each SPW contains 480 channels with a total bandwidth of 58.59 MHz, with a 0.141 MHz resolution per channel, corresponding to a velocity resolution of ~ 0.18 km/s. One SPW was reserved for high sensitivity continuum observations to aid in the self-calibration of the data. The total on-source integration time was ~ 43 minutes. A total of 42 and 43 antennas were included for the first EB and the remaining two EBs, respectively, and covered baselines from 15 to 740 m. All EBs used

the source J0510+1800 as their bandpass and flux calibrators and the source J0438+3004 as phase calibrator. Only one third of the proposed observations were performed, and both the RMS and beam size failed to meet our requested performance parameters. Therefore, the observations were deemed to QA2 SEMI-PASS and the data released by the observatory. However, the data quality already allows us to derive constraints on the S-chemistry as described below.

Data calibration was initially pipeline-processed. We then use the Common Astronomy Software Application package (CASA) version CASA 5.6.1-8 (McMullin et al. 2007) to reduce the data. Self-calibration was performed using the SPW reserved for continuum. We performed three iterations of phase self-calibrations, and then one amplitude self-calibration. After continuum subtraction with the CASA uvcontsub function, the data were CLEANed (Högbom 1974) using 3σ noise threshold and Briggs weighting with a robustness parameter of 0.5 for the main CS isotopologue and of 1 with a taper of 1" for the other lines to improve their imaging and SNR. The RMS per channel of all the observations presented in this study are listed in Table 2.

# 3. OBSERVATIONAL RESULTS

## 3.1. CS 2-1 fluxes and spatial distributions

Figure 1 displays the integrated intensity (zeroth moment) maps of the spatially resolved MAPS observations of the CS 2 – 1 rotational transition toward each of the five targeted disks. To build these maps we used the Python package bettermoments (Teague & Foreman-Mackey 2018) applied to the image cube available for
<table><caption>Table 2. List of Observations (molecular data from CDMS$^{(a)}$)</caption><thead><tr><th>Species</th><th>Transition</th><th>Frequency<br>(GHz)</th><th>E<sub>u</sub><br>(K)</th><th>Log<sub>10</sub>(A<sub>ij</sub>)<br>(s<sup>-1</sup>)</th><th>Source</th><th>RMS<sub>chan</sub><br>(mJy/beam)</th><th>Restored Beam<br>(" × ")</th><th>R<sub>max</sub><sup>(b)</sup><br>(")</th><th>S<sub>v</sub>Δ<sub>v</sub>(R<sub>max</sub><sup>(c)</sup>)<br>(mJy km/s)</th></tr></thead><tbody><tr><td colspan="10" style="text-align:center; font-weight:bold;">MAPS data (Project ID: 2018.1.01055.L)</td></tr><tr><td rowspan="5"><sup>12</sup>CS</td><td rowspan="5">2-1</td><td rowspan="5">97.98095</td><td rowspan="5">7.1</td><td rowspan="5">-4.7763</td><td>IM Lup</td><td>0.51</td><td>0.30 × 0.23</td><td>-80.2</td><td>3.0 ± 0.1</td><td>268 ± 11</td></tr><tr><td>GM Aur</td><td>0.46</td><td>0.39 × 0.27</td><td>5.1</td><td>2.1 ± 0.1</td><td>133 ± 3</td></tr><tr><td>AS 209</td><td>0.48</td><td>0.33 × 0.26</td><td>-66.7</td><td>0.9 ± 0.1</td><td>166 ± 4</td></tr><tr><td>HD 163296</td><td>0.42</td><td>0.31 × 0.24</td><td>-88.2</td><td>1.5 ± 0.1</td><td>302 ± 15</td></tr><tr><td>MWC 480</td><td>0.46</td><td>0.39 × 0.28</td><td>7.2</td><td>2.0 ± 0.1</td><td>48 ± 4</td></tr><tr><td rowspan="10">SO</td><td rowspan="5">2<sub>3</sub> - 1<sub>2</sub></td><td rowspan="5">99.29987</td><td rowspan="5">9.2</td><td rowspan="5">-4.9488</td><td>IM Lup</td><td>0.44</td><td>1.23 × 1.12</td><td>79.3</td><td>3.0 ± 0.1</td><td>~ 30</td></tr><tr><td>GM Aur</td><td>0.42</td><td>1.44 × 1.29</td><td>-5.9</td><td>2.1 ± 0.1</td><td>&lt; 10</td></tr><tr><td>AS 209</td><td>0.40</td><td>1.33 × 1.11</td><td>82.7</td><td>0.9 ± 0.1</td><td>&lt; 15</td></tr><tr><td>HD 163296</td><td>0.39</td><td>1.22 × 1.05</td><td>85.3</td><td>1.5 ± 0.1</td><td>&lt; 45</td></tr><tr><td>MWC 480</td><td>0.43</td><td>1.44 × 1.29</td><td>-3.6</td><td>1.5 ± 0.2</td><td>~ 14</td></tr><tr><td rowspan="5">5<sub>4</sub> - 4<sub>4</sub></td><td rowspan="5">100.02964</td><td rowspan="5">38.6</td><td rowspan="5">-5.9656</td><td>IM Lup</td><td>0.50</td><td>1.22 × 1.09</td><td>82.1</td><td>3.0 ± 0.1</td><td>&lt; 51</td></tr><tr><td>GM Aur</td><td>0.45</td><td>2.10 × 1.37</td><td>-29.0</td><td>2.1 ± 0.1</td><td>~ 4</td></tr><tr><td>AS 209</td><td>0.41</td><td>1.28 × 1.10</td><td>-80.9</td><td>0.9 ± 0.1</td><td>&lt; 18</td></tr><tr><td>HD 163296</td><td>0.41</td><td>1.26 × 1.07</td><td>85.0</td><td>1.5 ± 0.1</td><td>&lt; 53</td></tr><tr><td>MWC 480</td><td>0.41</td><td>2.08 × 1.38</td><td>-28.8</td><td>1.5 ± 0.2</td><td>&lt; 6</td></tr><tr><td rowspan="10">C<sub>2</sub>S</td><td rowspan="5">8<sub>7</sub> - 7<sub>6</sub></td><td rowspan="5">99.86652</td><td rowspan="5">28.1</td><td rowspan="5">-4.3562</td><td>IM Lup</td><td>0.49</td><td>1.22 × 1.09</td><td>82.1</td><td>3.0 ± 0.1</td><td>&lt; 51</td></tr><tr><td>GM Aur</td><td>0.45</td><td>2.10 × 1.37</td><td>-29.0</td><td>2.1 ± 0.1</td><td>&lt; 9</td></tr><tr><td>AS 209</td><td>0.41</td><td>1.29 × 1.10</td><td>-80.9</td><td>0.9 ± 0.1</td><td>~ 9</td></tr><tr><td>HD 163296</td><td>0.39</td><td>1.26 × 1.07</td><td>85.0</td><td>1.5 ± 0.1</td><td>&lt; 53</td></tr><tr><td>MWC 480</td><td>0.48</td><td>2.08 × 1.38</td><td>-28.8</td><td>1.5 ± 0.2</td><td>&lt; 5</td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">Complementary Cycle 6 unpublished ALMA data (project ID: 2018.1.01631.S)</td></tr><tr><td><sup>12</sup>CS</td><td>5-4</td><td>244.93564</td><td>35.3</td><td>-3.5271</td><td>MWC 480</td><td>3.9</td><td>0.71 × 0.45</td><td>-12.9</td><td>1.5 ± 0.2</td><td>98 ± 5</td></tr><tr><td><sup>13</sup>CS</td><td>5-4</td><td>231.220996</td><td>33.3</td><td>-3.6008</td><td>MWC 480</td><td>4.4</td><td>0.83 × 0.58</td><td>0.7</td><td>1.5 ± 0.2</td><td>&lt; 7</td></tr><tr><td>C<sup>34</sup>S</td><td>5-4</td><td>241.016088</td><td>27.8</td><td>-3.5568</td><td>MWC 480</td><td>3.4</td><td>0.86 × 0.59</td><td>6.8</td><td>1.5 ± 0.2</td><td>20 ± 5</td></tr><tr><td rowspan="4">H<sub>2</sub>CS</td><td>7<sub>16</sub> - 6<sub>15</sub></td><td>244.0485044</td><td>60.0</td><td>-3.6771</td><td>MWC 480</td><td>3.2</td><td>0.79 × 0.55</td><td>1.0</td><td>1.5 ± 0.2</td><td>29 ± 5</td></tr><tr><td>7<sub>26</sub> - 6<sub>25</sub></td><td>240.3820512</td><td>98.8</td><td>-3.7248</td><td></td><td>3.9</td><td>0.87 × 0.59</td><td>6.8</td><td></td><td>&lt; 13</td></tr><tr><td>7<sub>35</sub> - 6<sub>34</sub></td><td>240.3930370</td><td>164.6</td><td>-3.7760</td><td></td><td>4.0</td><td>0.86 × 0.59</td><td>6.3</td><td></td><td>&lt; 10</td></tr><tr><td>7<sub>34</sub> - 6<sub>33</sub></td><td>240.3937618</td><td>164.6</td><td>-3.7760</td><td></td><td>4.0</td><td>0.86 × 0.59</td><td>6.3</td><td></td><td>&lt; 10</td></tr><tr><td rowspan="2">OCS</td><td>19-18</td><td>231.06099</td><td>110.9</td><td>-4.4463</td><td>MWC 480</td><td>3.3</td><td>0.78 × 0.54</td><td>-0.02</td><td>1.5 ± 0.2</td><td>&lt; 5</td></tr><tr><td>20-19</td><td>243.21804</td><td>122.6</td><td>-4.3790</td><td></td><td>3.0</td><td>0.79 × 0.55</td><td>0.76</td><td></td><td>&lt; 13</td></tr><tr><td rowspan="3">SO<sub>2</sub></td><td>11<sub>57</sub> - 12<sub>48</sub></td><td>229.3476299</td><td>122</td><td>-4.7194</td><td>MWC 480</td><td>2.7</td><td>0.84 × 0.58</td><td>1.4</td><td>1.5 ± 0.2</td><td>18 ± 6</td></tr><tr><td>5<sub>24</sub> - 4<sub>13</sub></td><td>241.6157967</td><td>23.6</td><td>-4.0728</td><td></td><td>3.2</td><td>0.86 × 0.59</td><td>5.9</td><td></td><td>&lt; 6</td></tr><tr><td>5<sub>42</sub> - 6<sub>33</sub></td><td>243.0876473</td><td>53.1</td><td>-4.9886</td><td></td><td>2.9</td><td>0.79 × 0.55</td><td>1.2</td><td></td><td>&lt; 6</td></tr></tbody></table>

<footnote><sup>a</sup>https://cdms.astro.uni-koeln.de/cdms/portal/, Müller et al. (2001, 2005)</footnote>

<footnote><sup>b</sup> R<sub>max</sub> stands for the outer radius of the molecular line emission where 90% of the cumulative flux from the radial profiles is contained. The uncertainty is 1σ error.</footnote>

<footnote><sup>c</sup> S<sub>v</sub>Δ<sub>v</sub>(R<sub>max</sub>) corresponds to the flux density integrated out to the outer radius R<sub>max</sub> of the molecular line emission.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Zeroth moment maps and radial intensity profiles for the MAPS disk sample, ordered by increasing stellar mass, see Table 1, from left to right. *First row:* Zeroth moment maps of the dust continuum at 3 mm produced using an arcsinh color stretch for the AS 209 disk, and a power-law color stretch for all four other disks, to enhance the faint and extended emission. *Second row:* Zeroth moment maps of CS 2 − 1. Synthesized beams are shown in the lower left corner of each panel. *Third and fourth rows:* Radially de-projected and azimuthally averaged intensity profiles of the continuum and the CS 2 − 1 emission in y-linear and y-log scales, respectively. The vertical dashed gray lines indicate the outer radius $R_{max}$ of the molecular line emission where 90% of the cumulative flux from the radial profiles is contained within $1\sigma$ error. *Fifth row:* Integrated intensity spectra of CS 2 − 1. The uncertainties on radial and intensity profiles are calculated as the standard deviation on the mean in the radial annulus over which the emission was averaged, following the MAPS collaboration convention, described in detail in Law et al. (2021). So, these error bars do not include the absolute calibration uncertainty of 10%</figcaption></figure>
download in the MAPS data repository. We used an
hybrid mask combining a Keplerian mask (also avail-
able for download in the MAPS data repository) and
3σ clip to mask any pixels below this threshold. For
comparison, we also show the 3 mm continuum emis-
sion maps made in Sierra et al. (2021). Radially de-
projected and azimuthally averaged intensity profiles of
the continuum and CS 2 − 1 line are also shown for
each of the five MAPS disks. These were produced us-
ing the radial profile function of the Python package
GoFish (Teague 2019), considering the disk physical pa-
rameters (i.e., disk inclination, disk position angle, mass
of the central star, and distance) listed in Table 1. The
angular resolution is ~ 0.3'', i.e., ranging from 30 au
(HD 163296) to 49 au (MWC 480) depending on dis-
tance across the sample of MAPS disks. Finally, the CS
2-1 spectra are also depicted in Fig. 1 for each targeted
disks, showing a typical double-peaked profile indicative
of the Keplerian rotation of the disk.

Based on the radial intensity profiles of the CS 2 − 1 emission across the disk sample, central holes appear for IM Lup and HD 163296, with the largest radial hole extent found toward HD 163296. For the other three disks, the SNR and spatial resolution are not sufficient to infer the morphology of the inner disk emission. Beyond the inner holes we see a wide diversity in the morphology and extent of the CS 2 − 1 radial intensity profiles compared to the dust continuum. For instance, emission plateaus appear for IM Lup, GM Aur, and HD 163296, leading to outer CS 2 − 1 emission radii a factor of ≈ 2 larger than the dust continuum. The GM Aur disk – the only transition disk in the MAPS sample – is the only disk in which we see a tentative outer emission ring at ∼ 2''5 (i.e. ∼ 400 au).

Interestingly, the zeroth moment maps of the CS 2 − 1 emission show some asymmetries that do not appear in the dust emission. In particular, toward four of the sources (GM Aur, AS 209, HD 163296, and MWC 480), we find up to a factor 2 or 5σ of intensity difference between the near and far sides of the disks. Since the five MAPS disks have a non-zero inclination (see Table 1), the closest and farthest half disk sides, with respect to the disk semi-major axis, are defined as near and far disk sides, respectively. This is illustrated by the schematic views of the geometry of each disk that we show as insets in Fig. 2. For three of the five MAPS sources – namely AS 209, HD 163296, and GM Aur – the brightest CS emission sides coincide with the near side of the disks. Intriguingly, the reverse is observed for MWC 480, where the brightest CS emission side coincides with the far side of the disk. However, given the relatively low SNR, the robustness of these asymmetries is hard to assess. As for

IM Lup, which is the disk with the smallest CS 2-1 integrated intensity, we do not observe such asymmetries. These asymmetries are further discussed in Sect. 5.4.

3.2. *CS isotopologues and $H_2CS$ in MWC 480*

The ¹³CS and C³⁴S 5 – 4 rotational transitions were observed as part of our complementary ALMA program toward MWC 480 (see Sect. 2.2). We did not detect the ¹³CS 5 – 4 line neither with matched filtering method (VISIBLE, Loomis et al. 2018b) nor with line velocity shift and stacking techniques (GoFish, Teague 2019). The latter exploits the known geometry and velocity structure of the disk to de-project the rotation profile and combine Doppler shifted emission to a common centroid velocity reference frame. This results in a single disk-integrated spectrum for each transition. However, imaging the C³⁴S 5 – 4 line reveals a ~ 3 – 4σ detection that is shown in Fig. 3 and reported in Table 2. This detection is confirmed when we build the integrated spectrum of the line using the velocity shift and stacking methods of GoFish (see bottom panel in Fig. 3). These results are consistent with the ¹³CS and C³⁴S 6 – 5 observations reported in Le Gal et al. (2019b), where the C³⁴S line was tentatively detected toward the MWC 480 disk while ¹³CS was not.

Four H₂CS transitions (see Table 2) were also ob-
served as part of our complementary ALMA program
toward the MWC 480 disk. Among these four transi-
tions, only the H₂CS line with the lowest upper energy
level (i.e., H₂CS 7₁₆ – 6₁₅) is detected, with a ~ 5 – 6σ
detection. The zeroth moment map, the radially de-
projected and azimuthally averaged intensity profile,
and the shifted and stacked disk-integrated spectrum
of the H₂CS 7₁₆ – 6₁₅ detection are shown in Fig. 3. As
for the remaining three H₂CS lines, their non-detections
are not surprising regarding their upper energy levels
and line strengths (see Table 2, where we also report
upper limits in Table 2).

3.3. *Multi-line analysis*

To constrain the ¹²CS, C³⁴S, and H₂CS column densi-
ties and excitation temperatures toward the MWC 480
disk, we combined the new observations presented here
with complementary ALMA observations of ¹²CS (5 – 4
and 6 – 5), C³⁴S (6 – 5), and H₂CS (8₁₇ – 7₁₆, 9₁₉ – 8₁₈,
and 9₁₈–8₁₇), already published in Le Gal et al. (2019b).
Assuming optically thin lines and local thermal equi-
librium (LTE), we used a rotational diagram analysis
(Goldsmith & Langer 1999) to derive the disk-integrated
column densities and excitation temperatures of these
molecules. These quantities are derived from the disk-
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Upper and third rows: Zeroth moment maps of the CS 2-1 emission in the far and near sides of each of the five MAPS disks. Synthesized beams are shown in the lower left corner of the panels presenting the far side of each disk. Second and lower row: Radially de-projected and averaged intensity profiles of the furthest (blue) and near (orange) CS 2-1 emission sides of each disks compared to the azimuthally averaged intensity profiles (gray dotted line). The uncertainties on the radial intensity profiles are calculated as the standard deviation on the mean in the radial annulus over which the emission was averaged. Synthesized beams are shown by the gray error bar below each radial intensity profiles. The insets represent a schematic view of each disk inclination.</figcaption></figure>

integrated flux densities, as described in Le Gal et al. (2019b) and summarized below.

To illustrate that the LTE assumption is justified we show in Figure 4 the profiles of the main physical parameters used to build the MWC 480 disk physical structure along with the number density computed for CS and H₂CS with our corresponding published model Le Gal et al. (2019b). This allows us to show that the gas density in the main emitting molecular layers is larger than 10⁷ cm⁻<fnref n="3" />, i.e. well above the critical densities of the observed CS and H₂CS transitions, which justifies well the LTE assumption. For temperatures in the range 20–50 K, the critical densities for CS are in the ranges ~7 × 10⁴ – 3 × 10⁶ cm⁻<fnref n="3" /> (Shirley 2015), and in the range 2–4 × 10⁶ cm⁻<fnref n="3" /> for H₂CS, using scaled H₂CO collisional

rates from Wiesenfeld & Faure (2013) (see the Leiden Atomic and Molecular Database (LAMDA)<fnref n="3" />, van der Tak et al. 2020).

Assuming optically thin transitions, the disk-integrated flux densities S_νΔv, can be related to the column density of their respective upper energy state, N_u, as follows:

$$N_u = \frac{4\pi S_\nu \Delta v}{A_{ul} \Omega hc}, \qquad (1)$$

where S_ν is the flux density, Δv the line width, A_{ul} the Einstein coefficient, c the speed of light, and Ω the solid angle subtended by the source (e.g., Bisschop et al. 2008;

<footnote n="3">³ https://home.strw.leidenuniv.nl/ moldata/</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. Zeroth moment maps (top panels), radially de-projected and azimuthally averaged intensity profiles within 1σ - as in Fig. 1 - (middle panels), and shifted and stacked disk-integrated line spectra of the C$^{34}$S 5–4 (left) and H$_{2}$CS 7$_{16}$ – 6$_{15}$ (right) rotational transitions observed toward the MWC 480 disk with 1σ (bottom panels). These last uncer-tainties are calculated on a per channel basis, taking into account de-correlation along the spectral axis (see also Yen et al. 2016; Ilee et al. 2021).</figcaption></figure>

Loomis et al. 2018a). For this analysis, we use the disk flux densities S$_{ν}$Δv integrated out to the outer radius of the molecular line emission, referred to as R$_{max}$ in Table 2.

The disk-integrated column density, N$_{tot}$ , and excita- tion temperature, T , can then be derived from the upper level population, N$_{u}$ , which follows the Boltzmann dis- tribution:

$$
N_{\text{tot}} = \frac{N_u}{g_u} Q_{\text{rot}}(T) e^{E_u / k_B T}, \quad (2)
$$

with g$_{u}$ and E$_{u}$ the degeneracy and energy of the upper energy level u , respectively; k$_{B}$ the Boltzmann constant; and Q$_{rot}$ the partition function of the molecule, which for a diatomic molecule such as CS can be approximated by:

$$
Q_{\text{rot}}(T) \approx \frac{k_B T}{h B_0} + \frac{1}{3}. \qquad (3)
$$

In this expression h is the Planck constant and B$_{0}$ is the rotational constant of the molecule. For CS we used

B$_{0}$ = 24495.562 × 10 6 Hz (see CDMS). For H$_{2}$CS we interpolated the { T , Q$_{rot}$ ( T ) } values provided by CDMS.

Using Eq. 2 and appendix B of Le Gal et al. (2019b), the optical depth of a given transition at temperature T can be expressed as:

$$
\tau_{\nu} = \sqrt{\frac{4 \ln 2}{\pi} \frac{N_u A_{ul} c^3}{\Delta v_{\text{FWHM}} 8 \pi \nu^3} (e^{h\nu/k_B T_{\text{ex}}} - 1)}, \quad (4)
$$

where Δ v$_{FWHM}$ = √ 8 ln 2 σ$_{v}$ is the full width at half maximum of the observed transition. σ$_{v}$ is the width of the Gaussian fit, since for optically thin lines, the line profiles remain Gaussian.

As described in appendix B of Le Gal et al. (2019b), we can substitute Eq. 4 in C$_{τ}$ = τ 1 − e − τ , which corresponds to the ”optical depth correction factor” for a square line profile in case τ ≪ 1 (Goldsmith & Langer 1999). This allows us to build a likelihood function L (data, N$_{tot}$ , T$_{ex}$ ) that we used with the Python implementation emcee (Foreman-Mackey et al. 2013) of the affine-invariant en- semble sampler for Markov Chain Monte Carlo (MCMC) (Goodman & Weare 2010) to compute posterior proba- bility distributions for T$_{ex}$ and N$_{tot}$ . The following uni- form and permissive priors were assumed:

$$
T_{\text{ex}}(K) = \mathcal{U}(3, 300) \qquad (5)
$$

$$
N_{\text{tot}}(\text{cm}^{-2}) = \mathcal{U}(10^7, 10^{20}). \quad (6)
$$

3.4. Disk-integrated column densities in MWC 480

Using the method described in Section 3.3, we de- rived the disk-integrated column densities of CS, C$^{34}$S, and H$_{2}$CS in MWC 480. The random draws from the posterior distributions for each molecule are depicted in gray in Fig. 5. The uncertainties are derived from the median and 16th–84th percentiles of the posterior dis- tributions, respectively. The 16th and 84th percentiles are chosen as equivalent to ± 1 σ uncertainties on the fit. The results converged toward:

• T$_{ex}$ ≃ 22 . 7 +2 . 2 − 1 . 8 K and N$_{tot}$ ≃ 6 . 2 +0 . 5 − 0 . 5 × 10 12 cm − 2 for CS,

• T$_{ex}$ ≃ 14 . 2 +16 . 0 − 6 . 8 K and N$_{tot}$ ≃ 6 . 9 +15 . 6 − 2 . 9 × 10 11 cm − 2 for C$^{34}$S,

• T$_{ex}$ ≃ 29 . 6 +14 . 4 − 8 . 6 K and N$_{tot}$ ≃ 4 . 0 +2 . 1 − 0 . 9 × 10 12 cm − 2 for H$_{2}$CS.

This leads to N (CS)/N (C$^{34}$S) ≃ 9 +20 − 3 and N (CS)/N (H$_{2}$CS) ≃ 1 . 6 +0 . 8 − 0 . 4 . While the uncertainties on the former do not allow us to draw any firm conclusion, the latter is about a factor of two lower than previously found using fewer rotational transitions with a smaller
<figure><img src="image_5.png" /></figure>

**Figure 4.** First row: Visual extinction and UV flux profiles fed in our MWC 480 protoplanetary disk astrochemical model (Le Gal et al. 2019b). Second row: First two panels show the temperature and density 2D profiles fed in our MWC 480 protoplanetary disk astrochemical model. The third and fourth panels show the modeled number densities (i.e., absolute abundances) of CS, H₂CS. All panels are represented as functions of disk radius vs. height.

dynamic range in upper energy for H₂CS (i.e. here we have $E_u = 55.9 - 88.5$ K versus $E_u = 73.4 - 88.5$ K in Le Gal et al. 2019b). This illustrates the need for multiple line observations for a given molecule, to better constrain its excitation temperature and column density with rotational diagram methods.

### 3.5. Radially resolved column density of CS in MWC 480

Applying the same rotational diagram analysis to the radially de-projected and azimuthally averaged intensities, we compute the excitation temperature and column density of CS as a function of the distance from the star. All CS transitions were re-imaged to have matching beam sizes (i.e., ~ 0.5"). The results are presented in Fig. 6, along with the three CS lines radial intensity profiles. They are in good agreement with the disk-integrated results which are also depicted in Fig. 6 to facilitate the comparison. As expected from the CS radial intensity profiles, the CS column density decreases with increasing radius. One can note that the disk average values appear biased toward small distances from the central star because the bulk of the emission is coming from these small distances. So this is why the radially resolved column densities are preferred when derivable. While the temperature gradient is consistent with typical earlier derived radial temperature profiles, it is interesting to notice that the typical model temperatures

are higher than the ones derived from the Boltzmann analysis of the observations.

### 3.6. disk-integrated column density of CS in MAPS and literature

Next, we estimate the disk-integrated CS column densities for the remaining four MAPS disks. As each of these disks only has a single CS transition observed with MAPS (i.e., 2–1) we fix the excitation temperature to a minimum of 10 K and maximum of 30 K. This temperature range is based on the constraints derived for the MWC 480 disk (see Sect. 3.3), assuming that CS resides in similar temperature layers in each disk. We calculate the column densities associated with this temperature range using Equations 1 and 2. To enlarge our sample, we extended this calculation to another CS ALMA survey we performed in a sample of five additional disks (Le Gal et al. 2019b). The resulting CS disk-integrated column densities are shown in Fig. 7, sorted by stellar mass. The CS disk-integrated column density varies by ≈ 1.5 order of magnitude across the sample of disks, ranging from ≈ (0.2 – 6.0) × 10¹³ cm⁻². There are no obvious trends with stellar mass or spectral type. The two Herbig Ae stars, MWC 480 and HD 163296, are close to the sample average.

### 3.7. Upper limits and tentative detections of SO, C₂S, OCS, and SO₂
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Rotational diagrams of (i) the CS 2–1, 5–4, and 6–5 rotational transitions (top panel), (ii) the C$^{34}$S 5–4 and 6–5 rotational transitions (middle panel), (iii) H$_{2}$CS 7$_{16}$–6$_{15}$, 8$_{17}$–7$_{16}$, 9$_{19}$–8$_{18}$, and 9$_{18}$–8$_{17}$ rotational transitions (bottom panel), integrated over the outer radius of the molecular line emission, $R_{\text{max}}$, toward MWC 480. A 10% calibration uncertainty on the flux values is also included.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Top panel: Radially de-projected and azimuthally averaged intensity profiles of the three rotational transitions CS 2–1, 5–4, and 6–5 observed toward the MWC 480 disk. Middle and bottom panels: Radially de-projected and azimuthally averaged excitation temperature and column density profiles of the MCMC rotational diagram results applied to the aforementioned three CS lines. Median values and uncertainties based on the 16th, 50th, and 84th percentiles of the samples are depicted. For comparison, the disk-integrated CS column density and excitation temperature are over-plotted in green on the middle and bottom panels.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

**Figure 7.** Estimated CS, SO, and C₂S column densities disk-integrated up to R<sub>max</sub> for each MAPS disk (see Table 2) and computed for T<sub>ex</sub> = 10–30 K. For CS, the disk sample is extended to the additional five disks surveyed in Le Gal et al. (2019b). The disks are sorted by increasing stellar mass. The averaged column density of CS is represented by the dashed orange line and its standard deviation by the orange rectangle. Upper limits are indicated by the downward triangles.

While the only sulfur-bearing molecular transition targeted in a dedicated SPW within MAPS was ¹²CS 2–1, transitions of C₂S and SO were covered within the MAPS program (see Table 2). Assuming Keplerian emission and using matched filtering (Loomis et al. 2018b) these lines were not detected. To check for non-Keplerian emission, we also imaged these lines using the Briggs weighting with a robust parameter of 1 and a taper of 1" to improve the SNR and image quality. The corresponding zeroth moment maps and radially de-projected and azimuthally averaged intensity profiles are shown in Appendix A, in Fig. 12. The integrated intensities and upper limits for the non-detections are reported in Table 2 for the corresponding ¹²CS 2–1 emitting area. To estimate the upper limits of the SO and C₂S column densities, we use Eqs. 1 and 2 with the constraints on the excitation temperature of CS derived in Sect. 3.3 (as done in Sect. 3.6). The results are over-plotted in red and blue in Fig. 7.

Finally, two other oxygenated sulfur-bearing molecules, SO₂ and OCS, were observed toward the MWC 480 disk, as part of our complementary ALMA program (see Table 2). Figure 13, in Appendix A, shows the zeroth moment maps, radially de-projected and azimuthally averaged intensity profiles, and spectra of SO₂ and OCS respectively. Integrating the intensity over the CS 2–1 emitting area and FWHM, we find a 3σ tentative detection of SO₂ (see Table 2) that is also

distinguishable from the zeroth moment map where we see a subtle flux peak toward the disk. However, we do not reproduce a detection when using matched filtering (implemented in VISIBLE), nor with velocity shifting and stacking (implemented in GoFish). As for OCS, it is not detected and although its radial intensity profile shows a tentative peak toward the central star, the signal shown on the zeroth moment map is shifted from the disk location. Assuming LTE and using the CS excitation temperature derived toward the MWC 480 disk (see Sect. 3.3), we derived upper limits on the column densities that we compare with results from disk chemistry modeling in Fig. 8, presented in Sect. 4.4.

## 4. ASTROCHEMICAL MODELING

To further investigate the S-chemistry in protoplanetary disks, we computed a grid of astrochemical models tuned to the MWC 480 disk, which is the disk in which we observed the most S-bearing molecules (see Sect. 2.2).

### 4.1. Protoplanetary disk physical structure

Our fiducial protoplanetary disk astrochemical model is based on the MWC 480 disk model developed in Le Gal et al. (2019b). It consists of a 2D parametric physical structure in which the chemistry is post-processed (see Sect. 4.2). We consider here a simplistic physical structure in the sense that the disk is assumed to be symmetric azimuthally and with respect to the midplane.
Such a disk physical structure can thus be described in cylindrical coordinates centered on the star along two perpendicular axes characterizing the radius and height in the disk. Figure 4 shows the profiles of the gas temperature and density throughout the disk, for which the physical parameters used to compute the physical structure of MWC 480 are summarized in Table 3 and the parameterization is briefly summarized below, following Le Gal et al. (2019b).

For a given radius $r$ from the central star, the vertical temperature profile is computed following the formalism developed by Dartois et al. (2003):

$$T(z) = \begin{cases} T_{\text{mid}} + (T_{\text{atm}} - T_{\text{mid}}) \left[ \sin\left(\frac{\pi z}{2z_q}\right) \right]^{2\delta} & \text{if } z < z_q \\ T_{\text{atm}} & \text{if } z \ge z_q, \end{cases} \quad (7)$$

where $T_{\text{mid}}$ and $T_{\text{atm}}$ are respectively the midplane and atmosphere temperatures that vary as power law of the radii (Beckwith et al. 1990; Piétu et al. 2007; Le Gal et al. 2019b). $z_q = 4H$ with $H$ the pressure scale height that, assuming vertical hydrostatic equilibrium, can be expressed as follows:

$$H = \sqrt{\frac{k_B T_{\text{mid}} r^3}{\mu m_H G M_{\star}}}, \quad (8)$$

with $k_B$ the Boltzmann constant, $\mu = 2.4$ the reduced mass of the gas, $m_H$ the proton mass, $G$ the gravitational constant, and $M_{\star}$ the mass of the central star. The mid-plane temperature $T_{\text{mid}}$ is estimated following a simple irradiated passive flared disk approximation (e.g., Chiang & Goldreich 1997; Dullemond et al. 2001):

$$T_{\text{mid}}(r) \approx \left( \frac{\varphi L_{\star}}{8\pi r^2 \sigma_{\text{SB}}} \right)^{1/4}, \quad (9)$$

with $L_{\star} = 24 L_{\odot}$ the stellar luminosity (Andrews et al. 2013), $\sigma_{\text{SB}}$ the Stefan-Boltzman constant and $\varphi = 0.05$, a typical flaring angle value (e.g., Brauer et al. 2008; Baillié & Charnoz 2014). The atmosphere temperature, $T_{\text{atm}}$, is based on observational constraints. So here we consider $T_{\text{atm}} = T_{\text{atm,100 au}}(\frac{r}{100 \text{ au}})$, with $T_{\text{atm,100 au}}=48$ K from Guilloteau et al. (2011).

The disk is assumed to be in hydrostatic equilibrium. Thus, for a given vertical temperature profile, the vertical density structure is determined by solving the equation of hydrostatic equilibrium, as described from Eq. (17) to (20) in Le Gal et al. (2019b).

The surface density of the disk is assumed to follow a simple power law varying as $r^{-3/2}$ (Shakura & Sunyaev 1973; Hersant et al. 2009):

$$\Sigma(r) = \Sigma_{R_c} \left(\frac{r}{R_c}\right)^{-3/2}, \quad (10)$$

where $\Sigma_{R_c}$ is the surface density at the characteristic radius that can be expressed as function of the mass of the disk, $M_{\text{disk}}$, and its outer radius, $R_{\text{out}}$:

$$\Sigma_{R_c} = \frac{M_{\text{disk}} R_c^{-3/2}}{4\pi\sqrt{R_{\text{out}}}}, \quad (11)$$

with here $M_{\text{disk}} = 0.18 M_{\odot}$ (Guilloteau et al. 2011).

The visual extinction profile is derived from the hydrostatic density profile using the gas-to-extinction ratio of $N_{\text{H}}/A_{\text{V}} = 1.6 \times 10^{21}$ (Wagenblast & Hartquist 1989), with $N_{\text{H}} = N(\text{H})+2N(\text{H}_2)$ the vertical hydrogen column density of hydrogen nuclei. This gas-to-extinction ratio assumes a typical mean grain radius size of 0.1 µm and dust-to-mass ratio of 0.01. While the use of a grain size distribution including both large and small grains would be more realistic, its impact on the chemistry remains poorly constrained and would require a dedicated study such as, e.g., the one recently performed in Gavino et al. (2021). We therefore opt for the simpler approximation, which should be sufficient to provide an interpretative framework for the presented observations.

Finally, to compute the UV flux profile we multiplied the UV flux factor impinging on the disk with $e^{-x}$, where $x$ contains the visual extinction profile. The unattenuated UV flux factor, $f_{\text{UV}}$, at a given radius $r$ depends on both the photons coming directly from the central embedded star and on the photons that are downward-scattered by small grains in the upper atmosphere of the disk. Thus, following Wakelam et al. (2016), we consider:

$$f_{\text{UV}} = \frac{f_{\text{UV,R}_c}/2}{\left(\frac{r}{R_c}\right)^2 + \left(\frac{4H}{R_c}\right)^2}. \quad (12)$$

## 4.2. Protoplanetary disk chemical model

The disk chemistry is computed time-dependently in 1+1D based on the gas-grain astrochemical code Nautilus, which includes gas-phase, grain-surface, and grain-bulk chemistry (Wakelam et al. 2017; Le Gal et al. 2019b,a). This rate-equation gas-grain chemical code follows the formalism described in Hasegawa et al. (1992) and Hasegawa & Herbst (1993). We used the same chemical network as Le Gal et al. (2019b), which is based on the Kinetic Database for Astrochemistry (KIDA)$^4$, and includes recent updates (Vidal et al. 2017; Le Gal et al. 2017; Fuente et al. 2017; Le Gal et al. 2019b). It contains 589 gas-phase species and 540 solid-state species interacting together through a total of 13402 reactions. Chemical exchanges in between the

<footnote>$^4$ (http://kidaobs.u-bordeaux1.fr/)</footnote>
<table><caption>**Table 3.** Physical parameters used for our disk models</caption><thead><tr><th>Parameters</th><th>MWC 480<sup>a</sup></th></tr></thead><tbody><tr><td>Stellar mass: M<sub>*</sub> (M<sub>☉</sub>)</td><td>1.8</td></tr><tr><td>Disk mass: M<sub>d</sub> (M<sub>☉</sub>)</td><td>0.18</td></tr><tr><td>Characteristic radius: R<sub>c</sub> (au)</td><td>100</td></tr><tr><td>Outer cut-off radius: R<sub>out</sub> (au)</td><td>500</td></tr><tr><td>Density power-law index: γ</td><td>1.5</td></tr><tr><td>Midplane temperature at R<sub>c</sub><sup>b</sup>: T<sub>100au</sub>(K)</td><td>30</td></tr><tr><td>Atmosphere temperature at R<sub>c</sub>: T<sub>100au</sub>(K)</td><td>48</td></tr><tr><td>Surface density at R<sub>c</sub>:</td><td>5.7</td></tr><tr><td>Temperature power-law index: q</td><td>0.5</td></tr><tr><td>Vertical temperature gradient index: β</td><td>2</td></tr><tr><td>UV Flux: f<sub>UV,R<sub>c</sub></sub> (in Draine (1978)'s units)</td><td>8500<sup>c</sup></td></tr></tbody></table>

<sup>a</sup> These are the values used for the model developed in Le Gal et al. (2019b) and that we are using here to interpret the observations presented in the present work.

<sup>b</sup> The midplane temperature is estimated from Eq. (9), the luminosity and a typical flaring angle $\phi = 0.05$.

<sup>c</sup> from Dutrey et al. (2011), originally computed from the Kurucz (1993) ATLAS9 of stellar spectra.

gas-phase, grain-surface, and grain-bulk phases are included, with adsorption and desorption processes linking the gas and surface phases, and swapping processes linking the mantle and surface of grains. Several desorption mechanisms are taken into account: thermal desorption (Hasegawa et al. 1992), cosmic-ray induced desorption (Hasegawa & Herbst 1993), photodesorption and chemical desorption (e.g., Le Gal et al. 2017, and references therein). In the gas phase typical bi-molecular ion-neutral and neutral-neutral reactions are considered, as well as cosmic-ray induced processes, photoionizations and photodissociations caused by both stellar and interstellar UV photons.

First, we model the chemical evolution of a representative starless dense molecular cloud, with a characteristic age of 1 Myr (e.g., Elmegreen 2000; Hartmann et al. 2001) and typical constant physical conditions: grain and gas temperatures of 10 K, a total gas density of $2 \times 10^4 \text{ cm}^{-3}$, $\zeta = 1 \times 10^{-17} \text{ s}^{-1}$ per H$_2$, and a visual extinction of 30 mag. For this first stage model, we consider the initial chemical conditions to be close to diffuse gas conditions, i.e. all the elements are initially in atomic form (see Table 4) except hydrogen which is assumed to be already fully molecular. The elements taken into account in our simulation with an ionization potential lower than that of hydrogen (13.6 eV) are thus assumed to be initially singly ionized, see Table 4. The

<table><caption>**Table 4.** Initial Elemental Abundances</caption><thead><tr><th>Species</th><th>n<sub>i</sub>/n<sub>H</sub></th><th>Reference</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>0.5</td><td></td></tr><tr><td>He</td><td>9.0×10<sup>-2</sup></td><td>1</td></tr><tr><td>C<sup>+</sup></td><td>1.7×10<sup>-4</sup></td><td>2</td></tr><tr><td>N</td><td>6.2×10<sup>-5</sup></td><td>2</td></tr><tr><td>O</td><td>3.4 × 10<sup>-4</sup> – 1.1×10<sup>-4</sup></td><td>3</td></tr><tr><td>S<sup>+</sup></td><td>8.0 × 10<sup>-8</sup> – 1.5 × 10<sup>-5</sup></td><td>4</td></tr><tr><td>Si<sup>+</sup></td><td>8.0×10<sup>-9</sup></td><td>5</td></tr><tr><td>Fe<sup>+</sup></td><td>3.0×10<sup>-9</sup></td><td>5</td></tr><tr><td>Na<sup>+</sup></td><td>2.0×10<sup>-9</sup></td><td>5</td></tr><tr><td>Mg<sup>+</sup></td><td>7.0×10<sup>-9</sup></td><td>5</td></tr><tr><td>P<sup>+</sup></td><td>2.0×10<sup>-10</sup></td><td>5</td></tr><tr><td>Cl<sup>+</sup></td><td>1.0×10<sup>-9</sup></td><td>5</td></tr><tr><td>F<sup>+</sup></td><td>6.7×10<sup>-9</sup></td><td>6</td></tr></tbody></table>

NOTE— (1) Wakelam & Herbst (2008); (2) Jenkins (2009); (3) We varied the oxygen elemental abundance in this range to test the impact of the C/O ratio (see § 4.4). (see Sect. 4.4); (4) We varied the sulfur elemental abundance in this range to test the impact of the S/H ratio (see § 4.3); (5) Graedel et al. (1982); (6) Neufeld et al. (2015);

outcoming chemical gas and ice compositions of this representative parent molecular cloud serve as the initial chemistry for our 1+1D disk model, for which the physical parameters are described in Le Gal et al. (2019b). Second, we run the chemistry of our 1+1D disk model up to 1 Myr, the typical chemical age of a disk when dust evolution is not included (e.g., Cleeves et al. 2015), which is the case for our disk model. While the disk chemistry has not reached steady state at that time, its evolution is slow enough that the results presented here hold for a disk twice younger or older.

### 4.3. Impact of the S/H ratio

In the context of S-bearing molecules, a crucial parameter to study is the S/H elemental ratio, i.e. the total amount of S not locked into refractory compounds and thus available for the volatile S-chemistry. Figure 8 shows the modeled column densities of CS, H$_2$CS, SO, SO$_2$, C$_2$S, and OCS as function of distance from the central star for a range of C/O ratios (further described in Sect. 4.4) and for three different elemental S/H ratios: the usual highly depleted S-abundance value of $8.0 \times 10^{-8}$, corresponding to the "low metal" abundances from Graedel et al. (1982), an intermediate S-abundance value of $3.5 \times 10^{-6}$, corresponding to the value derived in PDR regions (Goicoechea et al. 2006; Le Gal et al. 2019a), and the solar abundance ($1.5 \times 10^{-5}$, Asplund et al. 2009). For comparison, estimated and upper limits of the column densities of the six S-bearing species we observed toward the MWC 480 disk are also displayed in Fig. 8.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** CS, H₂CS, C₂S, SO, SO₂, and OCS modeled column densities tuned to the MWC 480 disk, vertically integrated from the disk upper layer to the midplane and convolved to a resolution of 0.′5 to facilitate the comparison with the observations. The modeled column densities are shown by the solid lines investigating the impact of the C/O and S/H ratios. Observational error bars and upper limits derived toward the MWC 480 disk are indicated in gray. Note that the scales are replicated in all panels.</figcaption></figure>

While for a low S-elemental abundance, the column density of CS can be reproduced for C/O $\gtrsim$ 0.9, similarly to what has been found to reproduce the column densities of CH₃CN and HC₃N in the same disk (Le Gal et al. 2019a), the H₂CS column densities is underpredicted. On the contrary, no S-depletion, i.e. considering that all the solar S-abundance is available for S-chemistry in disks, allows the reproduction of H₂CS but cannot reproduce the column density of CS, which is then over-predicted. Models without S-depletion also require different C/O ratios to reproduce the column densities of C₂S and SO.

Since H₂CS is a more complex molecule as compared to CS, we suspect that its underproduction in our model is more likely to be due to missing formation pathways than differences in elemental abundances between models and observations. Experimental and theoretical chemical studies are needed to better constrain the formation pathways of H₂CS.

## 4.4. Impact of the C/O ratio

The relative gas-phase abundances of the chemical elements are known to strongly influence the chemistry of star-forming regions (van Dishoeck & Blake 1998). At the onset of star formation, substantial amount of the total budget of the main chemical elements, such as oxygen (O) and carbon (C), are locked in refractory materials. Furthermore, for some of them, huge uncertainties remain on the nature and the form of substantial part of

their reservoir. This is in particular the case for oxygen where ~40% of the O budget remains unaccounted for (Whittet 2010; Jones & Ysard 2019; Öberg & Bergin 2021) which results in a non-negligible uncertainty on the C/O ratio in the gas-phase.

In order to mimic the differential depletion of volatiles, we varied the C/O ratio from 0.5 to 1.5 (see Table 4). The impact of the gas-phase C/O ratio on the column densities of CS, H₂CS, SO, SO₂, C₂S, and OCS is shown in Fig. 8 and summarized in Table 5. As can be expected, for the carbonated sulfur molecules, i.e. the S-bearing species containing C-S bond, an O-poor chemistry (i.e., a high C/O ratio) results in higher column densities, while the reverse is seen for the oxygenated sulfur molecules (i.e. the molecule containing an O-S bond). This behavior is most prominent in the inner 1′′0 (i.e., ~160 au) of the disk for most molecules. Interestingly, the best model to fit the CS data is the most depleted S/H model.

The MAPS observations provide upper limits on SO, which allows us to calculate lower limits on the CS/SO ratio to which we can compare our model. Figure 9 shows how the modeled radial profile of the CS/SO column density ratio varies as a function of the elemental C/O ratio and total amount of sulfur. We only consider the two depleted S-abundance models, since we ruled out models with solar S in Section 4.3, based on comparisons between observed and modeled CS radial profiles. We find that the CS/SO ratio is highly sensitive to the
<table><caption>**Table 5.** Observed versus modeled disk-integrated column densities (in cm⁻²) and CS/SO ratio in MWC 480 out to 1.5".</caption><thead><tr><th></th><th>CS</th><th>H<sub>2</sub>CS</th><th>C<sub>2</sub>S</th><th>SO</th><th>SO<sub>2</sub></th><th>OCS</th><th>CS/SO<sup>†</sup></th><th>CS/SO<sup>‡</sup></th></tr></thead><tbody><tr><td>Observed value</td><td>6.6<sup>+0.5</sup><sub>-0.5</sub> × 10<sup>12</sup></td><td>4.0<sup>+2.1</sup><sub>-0.9</sub> × 10<sup>12</sup></td><td>~ 2.2 × 10<sup>12</sup></td><td>~ 1.2 × 10<sup>12</sup></td><td>~ 1.4 × 10<sup>15</sup></td><td>~ 3.7 × 10<sup>13</sup></td><td>~ 5.5</td><td></td></tr><tr><td>S/H=8 × 10<sup>-8</sup>; C/O=0.5</td><td>1.1 × 10<sup>12</sup></td><td>8.1 × 10<sup>99</sup></td><td>1.1 × 10<sup>11</sup></td><td>5.3 × 10<sup>12</sup></td><td>9.6 × 10<sup>12</sup></td><td>1.4 × 10<sup>11</sup></td><td>0.21</td><td>0.71</td></tr><tr><td>C/O=0.7</td><td>1.5 × 10<sup>12</sup></td><td>4.3 × 10<sup>10</sup></td><td>1.7 × 10<sup>11</sup></td><td>4.8 × 10<sup>12</sup></td><td>1.1 × 10<sup>13</sup></td><td>1.9 × 10<sup>11</sup></td><td>0.31</td><td>1.5</td></tr><tr><td>C/O=0.9</td><td>1.8 × 10<sup>12</sup></td><td>4.9 × 10<sup>10</sup></td><td>2.1 × 10<sup>11</sup></td><td>3.1 × 10<sup>12</sup></td><td>8.2 × 10<sup>12</sup></td><td>1.5 × 10<sup>11</sup></td><td>0.58</td><td>3.2</td></tr><tr><td>C/O=1.1</td><td>2.0 × 10<sup>13</sup></td><td>1.1 × 10<sup>11</sup></td><td>8.6 × 10<sup>11</sup></td><td>1.4 × 10<sup>12</sup></td><td>4.5 × 10<sup>12</sup></td><td>1.3 × 10<sup>11</sup></td><td>14</td><td>15</td></tr><tr><td>C/O=1.3</td><td>3.3 × 10<sup>13</sup></td><td>1.1 × 10<sup>11</sup></td><td>6.1 × 10<sup>11</sup></td><td>1.3 × 10<sup>12</sup></td><td>3.9 × 10<sup>12</sup></td><td>1.7 × 10<sup>11</sup></td><td>25</td><td></td></tr><tr><td>C/O=1.5</td><td>3.5 × 10<sup>13</sup></td><td>1.2 × 10<sup>11</sup></td><td>5.4 × 10<sup>11</sup></td><td>1.1 × 10<sup>12</sup></td><td>3.4 × 10<sup>12</sup></td><td>1.9 × 10<sup>11</sup></td><td>31</td><td></td></tr><tr><td>S/H=3.5 × 10<sup>-6</sup>; C/O=0.5</td><td>2.6 × 10<sup>13</sup></td><td>9.1 × 10<sup>10</sup></td><td>9.0 × 10<sup>11</sup></td><td>7.9 × 10<sup>13</sup></td><td>1.1 × 10<sup>14</sup></td><td>3.1 × 10<sup>13</sup></td><td>0.33</td><td>0.92</td></tr><tr><td>C/O=0.7</td><td>3.9 × 10<sup>13</sup></td><td>2.0 × 10<sup>12</sup></td><td>1.5 × 10<sup>12</sup></td><td>3.9 × 10<sup>13</sup></td><td>6.4 × 10<sup>13</sup></td><td>3.0 × 10<sup>13</sup></td><td>1.0</td><td>1.1</td></tr><tr><td>C/O=0.9</td><td>7.7 × 10<sup>13</sup></td><td>2.7 × 10<sup>12</sup></td><td>3.0 × 10<sup>12</sup></td><td>2.1 × 10<sup>13</sup></td><td>3.5 × 10<sup>13</sup></td><td>1.9 × 10<sup>13</sup></td><td>3.7</td><td>1.4</td></tr><tr><td>C/O=1.1</td><td>3.1 × 10<sup>14</sup></td><td>5.3 × 10<sup>12</sup></td><td>1.0 × 10<sup>13</sup></td><td>1.0 × 10<sup>13</sup></td><td>2.5 × 10<sup>13</sup></td><td>2.2 × 10<sup>13</sup></td><td>31</td><td>11</td></tr><tr><td>C/O=1.3</td><td>7.0 × 10<sup>14</sup></td><td>6.0 × 10<sup>12</sup></td><td>9.6 × 10<sup>12</sup></td><td>9.0 × 10<sup>12</sup></td><td>2.0 × 10<sup>13</sup></td><td>2.0 × 10<sup>13</sup></td><td>78</td><td></td></tr><tr><td>C/O=1.5</td><td>8.2 × 10<sup>14</sup></td><td>6.8 × 10<sup>12</sup></td><td>1.1 × 10<sup>13</sup></td><td>8.1 × 10<sup>12</sup></td><td>1.6 × 10<sup>13</sup></td><td>1.8 × 10<sup>13</sup></td><td>101</td><td></td></tr><tr><td>S/H=1.5 × 10<sup>-5</sup>; C/O=0.5</td><td>5.8 × 10<sup>13</sup></td><td>4.2 × 10<sup>11</sup></td><td>9.2 × 10<sup>11</sup></td><td>1.7 × 10<sup>14</sup></td><td>2.0 × 10<sup>14</sup></td><td>1.8 × 10<sup>14</sup></td><td>0.34</td><td></td></tr><tr><td>C/O=0.7</td><td>1.7 × 10<sup>14</sup></td><td>5.3 × 10<sup>12</sup></td><td>3.3 × 10<sup>12</sup></td><td>8.3 × 10<sup>13</sup></td><td>1.0 × 10<sup>14</sup></td><td>1.5 × 10<sup>14</sup></td><td>2.0</td><td></td></tr><tr><td>C/O=0.9</td><td>3.0 × 10<sup>14</sup></td><td>7.8 × 10<sup>12</sup></td><td>5.3 × 10<sup>12</sup></td><td>4.9 × 10<sup>13</sup></td><td>5.5 × 10<sup>13</sup></td><td>1.0 × 10<sup>14</sup></td><td>6.1</td><td></td></tr><tr><td>C/O=1.1</td><td>6.2 × 10<sup>14</sup></td><td>1.3 × 10<sup>13</sup></td><td>1.7 × 10<sup>13</sup></td><td>1.8 × 10<sup>13</sup></td><td>3.2 × 10<sup>13</sup></td><td>9.7 × 10<sup>13</sup></td><td>34</td><td></td></tr><tr><td>C/O=1.3</td><td>1.1 × 10<sup>15</sup></td><td>1.5 × 10<sup>13</sup></td><td>1.9 × 10<sup>13</sup></td><td>1.6 × 10<sup>13</sup></td><td>2.3 × 10<sup>13</sup></td><td>7.7 × 10<sup>13</sup></td><td>69</td><td></td></tr><tr><td>C/O=1.5</td><td>1.5 × 10<sup>15</sup></td><td>1.6 × 10<sup>13</sup></td><td>2.0 × 10<sup>13</sup></td><td>1.4 × 10<sup>13</sup></td><td>1.7 × 10<sup>13</sup></td><td>6.1 × 10<sup>13</sup></td><td>107</td><td></td></tr></tbody></table>

† Models w/o dust settling.

‡ Models with dust settling (see Sect. 4.5).

C/O ratio; a change in C/O from 0.5 to 1.5 increases the CS/SO ratio by up to 4 orders of magnitude. This is consistent with previous disk modeling results from (Semenov et al. 2018), as well as with cloud chemistry predictions (e.g. Bergin et al. 1997; Nilsson et al. 2000). We can compare these model results with our observational lower limit of > 5.5 (see Table 5). Based on the visual comparison between models and data in Fig. 9, the C/O ratio needs to be ≥ 0.9 in order to reproduce the CS/SO ratio observation in the MWC 480 disk. We also provide disk-integrated CS/SO ratios for the relevant disk models in Table 5, and these confirm that only models with C/O > 0.9 are consistent with observations.

### 4.5. Impact of dust evolution

Dust evolution and in particular dust settling can have a profound effect on disk chemistry, including on the ratios of molecules that have been proposed as tracers of C/O. Wakelam et al. (2019) recently explored the impact of several disk parameters on the vertically-integrated column densities of a set of molecules observed in the DM Tau disk. They found that dust settling can have a strong impact on the disk chemistry and can, in particular, enhance the chemical abundances of several carbon-bearing molecules such as CH$_3$CN and HC$_3$N. To test if dust settling could change our conclusions

on the C/O ratio in the MWC 480 disk, we ran an additional set of models, for a smaller grid of C/O values (from 0.5 to 1.1) and including similar dust settling as the fiducial one proposed in Wakelam et al. (2019) (i.e. their E2 model, where the settling occurs at $z/h = 1$). The results are depicted in Fig. 10.

Comparing Figs. 9 and 10, we see that dust settling indeed influences the variation of CS/SO ratio as function of C/O ratio. With dust settling, the models with varying C/O produce column density ratios within 1 order of magnitude for radius ≥ 0.3". While without dust settling, at least for the inner disk (i.e. radius < 0.5"), the models show a spread of 4 to 8 orders of magnitude. Therefore, it seems that if dust settling is present in a disk, one can only derive whether C/O is smaller or larger than 1. Furthermore, one can note that for the outer disk (> 0.7"), the models with varying C/O are almost indistinguishable when dust settling is present (in particular for the most depleted S/H model).

In summary, despite the apparent dust settling impact on the CS/SO ratio found with our modeling, we find that the results can be consistent with a high C/O ratio under the specific conditions assumed in the modeling. While out of the scope of the present study, a larger and dedicated deeper study that simultaneously explores the impact of dust settling and different C/O
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Caculated N(CS)/N(SO) column density ratios for a grid of models tuned to the MWC 480 disk investigating the impact of C/O and S/H ratios. Observations toward the MWC 480 disk</figcaption></figure>

are indicated by the gray horizontal lower limits. Because SO is not detected, the spatial distribution of SO in the disk is unknown. So, we extracted the upper limit on the SO emission for the exact same region of the disk area in which the CS emission is detected (see Section 3.7).

ratios on global disk chemistry would be very interesting to pursue in the future.

## 5. DISCUSSION

### 5.1. *Is the CS/SO column density ratio a good C/O ratio proxy?*

As shown in Fig. 9, our modeling results suggest that the CS/SO column density ratio is a promising probe of the C/O ratio in disks (see Section 4.4). We find that an elevated C/O ratio (i.e. a super-solar C/O) is required for the MWC 480 disk in order to reproduce the observed CS/SO ratio. A C/O ratio $\ge$ 0.9 seems reasonable for the MWC 480 chemistry as it results in both a detectable column density of the S-organic compounds H$_2$CS and a good match to the observed column densities of ni-

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Same as Fig. 9 but with models considering dust settling as described in Section 4.5.</figcaption></figure>

triles (Le Gal et al. 2019a). Moreover, while dust settling seems to impact the CS/SO ratio with varying C/O, we are still finding results converging toward a super-solar C/O. It is worth noticing that this is in very good agreement with the results found from other molecules probed within the MAPS program, i.e., CO isotopologues and C$_2$H. Using an independent disk model, we found that a super-solar C/O is also required to reproduce the CO isotopologues and C$_2$H observations in the same disk (Bosman et al. 2021; Alarcón et al. 2021). Furthermore, we also checked the predicted water vapor abundances from our models, and in the elevated C/O case they are consistent with the upper limits provided by the WISH project (van Dishoeck et al. 2021), while water is over-predicted in the low C/O models. Thus, 1) there seems to be a robust convergence toward an elevated C/O ratio and 2) the CS/SO ratio appears to be an additional and independent good probe of the C/O ratio.

In Figure 11, we show the lower limits found for the N(CS)/N(SO) column density ratio derived in each of the protoplanetary disks observed with MAPS. Among this sample, the ratio varies by a factor ~ 10, which
leads to similar C/O ratio constraints for each of the disks, i.e. a super-solar C/O if we consider that our disk model results holds for the other four MAPS disks. However, these preliminary results would need to be corroborated by deeper upper limits on SO and further modeling investigations which would be address in forthcoming studies. Additional CS/SO measurements toward a larger sample of protoplanetary disks would also be a good way to measure how common is the C/O ratio expected to be elevated in disks.

<figure><img src="image_12.png" /><figcaption>Figure 11. N(CS)/N(SO) column density ratio derived from the MAPS observations.</figcaption></figure>

## 5.2. Interpretation of disk S/H ratio

In protoplanetary disks, the S/H elemental ratio has been much less studied, and therefore less well constrained, than the C/O and C/H ratios. As of today, we still do not know what the major S reservoir(s) in disks is (are) and in which form it resides (solid or gaseous). However, this is an important parameter to constrain as well, not only to solve the current disk modeling tension found to interpret the high H₂CS/CS ratio in MWC 480, but also because, more generally, many S-bearing species are observed in comets and do play an important role in the building-up of pre-biotic molecules and on planet habitability.

Recently, based on the abundances of B star photospheres, Kama et al. (2019) found that ~ 81 – 97% of the S-budget should be locked in disk refractory material. Following the findings of Keller et al. (2002), the former authors proposed that most of the sulfur should be locked in the form of solid FeS in disks, rather than in polymeric S<sub>n</sub> (n = 2 – 8) molecules, where the lat-

ter has been proposed for decades as a potential hidden S-reservoir in the ISM (e.g., Wakelam et al. 2004). However, the observations of solid FeS has not yet been confirmed in disks (Keller et al. 2002) whereas hints of solid S<sub>n</sub> were recently reported on the comet 67P/C-G with the detections of S<sub>2</sub>, S<sub>3</sub>, and S<sub>4</sub> from the Rosetta mission (Calmonte et al. 2016), although these detections constitute only ~ 0.2% of the total detected S-content of 67 P/C-G.

The comparison with comets’ sulfur-bearing molecules content is, however, instructive. In particular, it is interesting to note that in the inventory of the molecular abundances detected in comets compiled by Bockelée-Morvan & Biver (2017), a H₂CS/CS ratio of ≈ 0.45 is reported, i.e. close to the value we measured in the MWC 480 protoplanetary disk (see Section 3.4 and Le Gal et al. 2019b). Another relevant point stressed in the Bockelée-Morvan & Biver (2017) review is that one of the most abundant S-bearing molecules detected in comets is OCS, which could be another potential S-reservoir. OCS is indeed the only S-bearing molecule unambiguously detected in ice mantles so far (Geballe et al. 1985; Palumbo et al. 1995). Furthermore, astrochemical shock modeling benchmarked to protostellar shock observations predicted that ≥ 50% of the sulfur ice reservoir resides in OCS (Podio et al. 2014; Holdship et al. 2016). The latter is therefore a promising S-reservoir to interpret the high H₂CS/CS ratio we observed toward MWC 480 that is not reproducible with our model (see Section 4).

## 5.3. Could the sulfur organic chemistry be underappreciated in models?

In our model, and as previously described in Le Gal et al. (2019b), H₂CS is mainly formed from the following neutral-neutral and electronic dissociative recombination gas-phase reactions:

$$S + CH_3 \longrightarrow H_2CS + H, \qquad (13)$$

$$H_3CS^+ + e^- \longrightarrow H_2CS + H, \qquad (14)$$

with H₃CS⁺ originating from the S⁺ + CH₄ reaction.

It is also formed, for a smaller contribution, from gas-grain chemistry where H₂CS is produced from successive hydrogenation on icy dust mantles and released for ~ 1% in the gas phase by chemical reactive desorption. However, in our current model, there is no consideration of OCS grain surface processing that could lead to the formation of S-organics such as H₂CS and thus maybe help in better reproducing the observations. Because H₂CS is a more complex species compared to CS, we
suspect that its underproduction is more likely due to missing formation pathways than an un-depleted S gas-phase reservoir. The latter is inconsistent with CS observations (see Section 4.3). Laboratory experiments and theoretical chemical calculations for such mechanisms are required to further test our hypothesis. While the S-reservoir could be changing from the inner to the outer disk regions, dedicated disk resolved S-observations are also needed to further investigate the nature and identity of the S-reservoir in disks, and how and if there is any chemical inheritance from molecular cloud stage to the planet-forming environment.

### 5.4. CS disk-emission asymmetries

In Sect. 3.1, we highlight asymmetries in the CS 2–1 emission spatial distribution toward four of the five targeted disks. Intriguingly, these asymmetries are not observed in the other molecular lines targeted within MAPS nor in higher CS transitions (i.e., the 5–4 and 6–5 transitions published in Le Gal et al. 2019b). We also do not see any such asymmetries in the dust emission of any of the targeted disks. So these asymmetries appear to be different than, for instance, the one detected in CS J = 7–6 toward the HD 142527 disk (van der Plas et al. 2014). However, we should note that in the present study we are reporting CS 2–1 asymmetries which are likely emitting from disk layers closer to the midplane than the CS 7–6 would. Toward the GM Aur, AS 209, and HD 163296 disks, the CS 2–1 emission asymmetries seem to correlate with their respective disk inclination, i.e. the CS 2–1 emission is brightest in the near side of the disk (see Fig. 2). However, toward MWC 480, we observe the reverse: the CS 2–1 emission is brightest in the far side of the disk. Here we investigate what could cause such asymmetries.

The HD 163296 disk is known to harbor both a jet (Grady et al. 2000) and a disk wind. The latter was discovered through ¹²CO observations by Klaassen et al. (2013) and is also further characterized using CO isotopologue observations as a part of MAPS (Booth et al. 2021b). According to the geometry of jet, wind and viewing angle proposed in Ellerbroek et al. (2014), the far side of the disk is supposed to be viewed through the disk’s wind and jet, which is also the disk side where we found the CS 2–1 emission to be the weakest. Thus, a speculative interpretation could be that the disk’s wind and jet impact the CS content, or simply the emission of the 2–1 line, and could therefore explain the decline in CS 2–1 flux in the side of the disk affected by the wind and jet. For instance, if the wind impacts the local C/O ratio thus it could impact the local disk chemistry and maybe the total amount of the CS bulk lying closer

to the midplane layer; or/and if the wind is also made of dust, the line emission below the wind with respect to the observer could be hampered due to dust wind opacity. However, follow-up observations of CS toward both the HD 163296 disk and wind would be required to test this hypothesis.

Similarly, MWC 480 is known to be driving a bipolar jet aligned with the disk semi-minor axis (Grady et al. 2010). Notably, the jet flow appears denser in the SW direction which could explain the decrease in CS 2–1 emission we observe in the same direction. A better characterization of this jet is required to assess how it could impact the CS 2–1 disk emission.

Hints for a disk wind are also found toward the AS 209 disk (e.g., Banzatti et al. 2019; Fang et al. 2018, and references therein), but, to our knowledge, the orientation and spatial distribution of the latter remain to be determined. Another point to mention about the AS 209 disk is that its west half side is known to be cloud-contaminated (Öberg et al. 2011; Huang et al. 2016; Teague et al. 2018a). While it strongly impacts the ¹²CO 2–1 and HCO⁺ 1–0 line emission in this disk (see Fig. 2 and 4 in Law et al. 2021), this cloud-contamination does not match with the CS 2–1 asymmetries we are finding in the present work.

As for GM Aur, Macías et al. (2018) discuss the possibility of a radio jet. Furthermore, the GM Aur disk is also the only transitional disk of our sample – i.e., the only one with a central dust cavity – and, as characterized by its complex gas structures, it is known to be affected by much more prominent gas dynamics than the other disks of our sample (e.g., Huang et al. 2021, and reference therein). Therefore, the CS 2–1 asymmetry of this disk is probably the least difficult to justify but would require further observations to be linked with the other gas kinematics features observed in this disk.

While all these hypothesis seem appealing, further investigations are required to truly determine the nature of these CS 2–1 asymmetries and in particular to identify if they are tracing one specific characteristic of disk evolution or if they could be explained by multiple phenomena.

## 6. CONCLUSION

We presented ALMA observations of S-bearing molecules observed toward the five protoplanetary disks targeted by the MAPS ALMA Large Program, orbiting the IM Lup, GM Aur, and AS 209 T Tauri stars, and the two Herbig Ae stars HD 163296 and MWC 480.

Our main findings are summarized below:

1. The CS 2–1 line was observed within MAPS and detected toward all five disks displaying a variety
of radial intensity profiles and spatial distributions across the sample, including intriguing apparent azimuthal asymmetries.

2. Using complementary ALMA observations of CS 5–4 and 6–5 in one of the disks, i.e., the MWC 480 disk, allows us to assess the CS column density across the full sample, assuming a temperature in the range 10-30 K, which results in $N_{\text{tot}}(\text{CS}) \approx (0.2–5) \times 10^{13} \text{ cm}^{-2}$.

3. $C_2S$ and SO lines were also covered within MAPS. While no detection can be robustly claimed from these observations, we provide upper limits on their column densities, with upper limits in the range $10^{12} – 10^{13} \text{ cm}^{-3}$ for $C_2S$ and $[1–5] \times 10^{12} \text{ cm}^{-3}$ for SO. In particular, we used the upper limit on SO to derive lower limits on the CS/SO ratio across the MAPS sample, which is found to range from ~4 to 14.

4. Using complementary ALMA programs, we find $N(\text{H}_2\text{CS})/N(\text{CS}) \approx 2/3$ in MWC 480. This high ratio suggests that substantial S-reservoirs in disks may be in the form of S-organics (i.e., $C_xH_yS_z$).

5. Using astrochemical disk models, we find that the CS/SO ratio is a promising probe for the elemental C/O ratio. CS/SO varies by more than 4 orders of magnitude when C/O varies from 0.5 to 1.5.

6. For MWC 480, without considering dust settling, we find C/O $\ge 0.9$, consistent with constraints from nitriles observations (Le Gal et al. 2019a). When considering dust settling, one can only derive whether C/O is smaller or larger than 1, but this remains consistent with a high C/O ratio under the specific conditions assumed in the modeling. More interestingly, this is confirmed with independent disk chemistry models predicting super-solar C/O based on the CO and $C_2H$ MAPS data (Bosman et al. 2021).

7. We find a depleted gas-phase S/H ratio, suggesting either that part of the sulfur reservoir is locked in solid phase or that it remains in an unidentified gas-phase reservoir. More sulfur observations are required to confirm this and, to a larger extent, to identify the nature of the S-reservoir(s).

Together these results illustrate the importance of sulfur chemistry in protoplanetary disks, demonstrating that, not only, sulfur-bearing molecules observations in such disks can serve to constrain the sulfur chemistry itself and its reservoir(s), but

also that sulfur-bearing molecules are powerful tools to constrain other fundamental parameters, such as the elemental C/O ratio. Furthermore, sulfur-bearing molecules seem to uniquely probe disk gas substructures, but this require deeper observations to be further investigated and confirmed. Therefore, to fully comprehend the role of sulfur in disks, further theoretical and observational studies on the sulfur chemistry in disks are still required.

## ACKNOWLEDGMENTS

We thank the referees for their valuable comments which helped us to improve the quality of this manuscript. R.L.G. also thanks Clément Baruteau for helpful discussions. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01055.L and ADS/JAO.ALMA#2018.1.01631.S. 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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

R.L.G. acknowledges support from CNES fellowship grant. K.I.Ö. and R.L.G. acknowledge support from a Simons Foundation award (SCOL # 321183, KO) and an NSF AAG Grant (#1907653). R.T. acknowledges support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow. C.J.L. acknowledges funding from the National Science Foundation Graduate Research Fellowship under Grant DGE1745303. C.W. acknowledges financial support from the University of Leeds, STFC and UKRI (grant numbers ST/R000549/1, ST/T000287/1, MR/T040726/1). E.A.B. and A.D.B. acknowledge support from NSF AAG Grant #1907653. F.M. acknowledges support from ANR of France under contract ANR-16-CE31-0013 (Planet-Forming-Disks) and ANR-15-IDEX-02 (through CDP “Origins of Life”). S. M. A. and J. H. acknowledge funding support from the National Aeronautics and Space Administration under Grant No. 17-XRP17 2-0012 issued through the Exoplanets Research Program. J.H. acknowledges support for this work provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51460.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Inc., for NASA, under contract
NAS5-26555. Y.A. acknowledges support by NAOJ
ALMA Scientific Research Grant Numbers 2019-13B,
and Grant-in-Aid for Scientific Research 18H05222 and
20H05847. A.S.B. acknowledges the studentship funded
by the Science and Technology Facilities Council of
the United Kingdom (STFC). G.C. is supported by
NAOJ ALMA Scientific Research Grant Code 2019-13B.
J.B.B. acknowledges support from NASA through the
NASA Hubble Fellowship grant \#HST-HF2-51429.001-
A, awarded by the Space Telescope Science Institute,
which is operated by the Association of Universities for
Research in Astronomy, Inc., for NASA, under contract
NAS5-26555. L.I.C. gratefully acknowledges support
from the David and Lucille Packard Foundation and
Johnson & Johnson's WiSTEM2D Program. IC was
supported by NASA through the NASA Hubble Fel-
lowship grant HST-HF2-51405.001-A awarded by the
Space Telescope Science Institute, which is operated
by the Association of Universities for Research in As-
tronomy, Inc., for NASA, under contract NAS5-26555.
V.V.G. acknowledges support from FONDECYT Ini-
ciación 11180904 and ANID project Basal AFB-170002.
J.D.I. acknowledges support from the Science and Tech-
nology Facilities Council of the United Kingdom (STFC)
under ST/T000287/1. H.N. acknowledges support by
NAOJ ALMA Scientific Research Grant Code 2018-
10B and Grant-in-Aid for Scientific Research 18H05441.
K.R.S. acknowledges the support of NASA through
Hubble Fellowship Program grant HST-HF2-51419.001,
awarded by the Space Telescope Science Institute, which
is operated by the Association of Universities forRe-
search in Astronomy, Inc., for NASA, under contract
NAS5-26555. T.T. is supported by JSPS KAKENHI
Grant Numbers JP17K14244 and JP20K04017. Y.Y. is
supported by IGPEES, WINGS Program, the Univer-
sity of Tokyo. K.Z. acknowledges the support of the
Office of the Vice Chancellor for Research and Gradu-
ate Education at the University of Wisconsin – Madi-
son with funding from the Wisconsin Alumni Research
Foundation, and the support of NASA through Hubble
Fellowship grant \#HST-HF2-51401.001. awarded by the
Space Telescope Science Institute, which is operated by
the Association of Universities for Research in Astron-
omy, Inc., for NASA, under contract NAS5-26555.

Facilities: ALMA

Software: **CASA** (McMullin et al. 2007), **Astropy** (Astropy Collaboration et al. 2013), **Matplotlib** (Hunter 2007), **NumPy** (van der Walt et al. 2011), **emcee** (Foreman-Mackey et al. 2013), **SciPy** (Jones et al. 2001–), **scikit-image** (van der Walt et al. 2014), **Gofish** (Teague 2019), **VISIBLE** (Loomis et al. 2018b)
<figure><img src="image_13.png" /><figcaption>**Figure 12.** *First and third rows:* Zeroth moment maps of C<sub>2</sub>S 8<sub>7</sub> − 7<sub>6</sub> (upper row) and SO 2<sub>3</sub> − 1<sub>2</sub> (third row). Synthesized beams are shown in the lower left corner of each panel. *Second and fourth rows:* Radially de-projected and azimuthally averaged intensity profiles of the C<sub>2</sub>S 8<sub>7</sub> − 7<sub>6</sub> (second row) and SO 2<sub>3</sub> − 1<sub>2</sub> (fourth row) emission lines.</figcaption></figure>

# APPENDIX

## A. IMAGING OF SO, C₂S, OCS, AND SO₂

In this section, we present the zeroth moment map, radially de-projected and azimuthally averaged intensity profile, and shifted and stacked disk-integrated spectrum zeroth moment maps, radial intensity profiles for the strongest C₂S and SO lines serendipitously covered by MAPS (Fig. 12) and for the SO₂ and OCS lines covered in our complementary Cycle 6 ALMA data (Fig. 13). All of them are derived for the same emitting area as the detected CS 2 − 1 transition (see Table 2).

## REFERENCES
<figure><img src="image_14.png" /></figure>

**Figure 13.** Zeroth moment map (top panels), radially de-projected and azimuthally averaged intensity profile within 1σ - similarly too Fig. 1 and 3 (middle panels), and shifted and stacked disk-integrated line spectra within 1σ (bottom panels). These last uncertainties are calculated on a per channel basis, taking into account de-correlation along the spectral axis (see also Yen et al. 2016; Ilee et al. 2021).

Aikawa, Y., van Zadelhoff, G. J., van Dishoeck, E. F., &
Herbst, E. 2002, A&A, 386, 622,
doi: 10.1051/0004-6361:20020037

Alarcón, F., Bosman, A., Bergin, E., et al. 2021, arXiv e-prints, arXiv:2109.06263.
https://arxiv.org/abs/2109.06263

Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20, doi: 10.1051/0004-6361/201629929

Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner, D. J. 2013, ApJ, 771, 129, doi: 10.1088/0004-637X/771/2/129

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, The Astrophysical Journal, 869, L41, doi: 10.3847/2041-8213/aaf741

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068

Baillié, K., & Charnoz, S. 2014, ApJ, 786, 35, doi: 10.1088/0004-637X/786/1/35

Banzatti, A., Pascucci, I., Edwards, S., et al. 2019, ApJ, 870, 76, doi: 10.3847/1538-4357/aaf1aa

Beck, T. L., & Bary, J. S. 2019, The Astrophysical Journal, 884, 159, doi: 10.3847/1538-4357/ab4259

Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924, doi: 10.1086/115385
Bergin, E. A., Aikawa, Y., Blake, G. A., & van Dishoeck, E. F. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 751.
https://arxiv.org/abs/astro-ph/0603358

Bergin, E. A., Goldsmith, P. F., Snell, R. L., & Langer, W. D. 1997, ApJ, 482, 285, doi: 10.1086/304108

Bergner, J. B., Guzmán, V. G., Öberg, K. I., Loomis, R. A., & Pegues, J. 2018, ApJ, 857, 69, doi: 10.3847/1538-4357/aab664

Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, ApJ, 876, 25, doi: 10.3847/1538-4357/ab141e

Bisschop, S. E., Jørgensen, J. K., Bourke, T. L., Bottinelli, S., & van Dishoeck, E. F. 2008, A&A, 488, 959, doi: 10.1051/0004-6361:200809673

Biver, N., Moreno, R., Bockelée-Morvan, D., et al. 2016, A&A, 589, A78, doi: 10.1051/0004-6361/201528041

Bockelée-Morvan, D., & Biver, N. 2017, Philosophical Transactions of the Royal Society of London Series A, 375, 20160252, doi: 10.1098/rsta.2016.0252

Bockelée-Morvan, D., Crovisier, J., Mumma, M. J., & Weaver, H. A. 2004, The composition of cometary volatiles, ed. M. C. Festou, H. U. Keller, & H. A. Weaver, 391

Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, 53, 541, doi: 10.1146/annurev-astro-082214-122348

Booth, A. S., van der Marel, N., Leemker, M., van Dishoeck, E. F., & Ohashi, S. 2021a, arXiv e-prints, arXiv:2104.08908. https://arxiv.org/abs/2104.08908

Booth, A. S., Walsh, C., Kama, M., et al. 2018, A&A, 611, A16, doi: 10.1051/0004-6361/201731347

Booth, A. S., Tabone, B., Ilee, J. D., et al. 2021b, arXiv e-prints, arXiv:2109.06586.
https://arxiv.org/abs/2109.06586

Bosman, A. D., Alarcón, F., Bergin, E. A., et al. 2021, arXiv e-prints, arXiv:2109.06221.
https://arxiv.org/abs/2109.06221

Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859, doi: 10.1051/0004-6361:20077759

Calmonte, U., Altwegg, K., Balsiger, H., et al. 2016, MNRAS, 462, S253, doi: 10.1093/mnras/stw2601

Chen, Y. J., Juang, K. J., Nuevo, M., et al. 2015, ApJ, 798, 80, doi: 10.1088/0004-637X/798/2/80

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368, doi: 10.1086/304869

Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204, doi: 10.1088/0004-637X/799/2/204

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2018, ApJ, 865, 155, doi: 10.3847/1538-4357/aade96

Czekala, I., Loomis, R. A., Teague, R., et al. 2021, arXiv e-prints, arXiv:2109.06188.
https://arxiv.org/abs/2109.06188

Dartois, E., Dutrey, A., & Guilloteau, S. 2003, A&A, 399, 773, doi: 10.1051/0004-6361:20021638

Draine, B. T. 1978, ApJS, 36, 595, doi: 10.1086/190513

Drozdovskaya, M. N., van Dishoeck, E. F., Jørgensen, J. K., et al. 2018, MNRAS, 476, 4949, doi: 10.1093/mnras/sty462

Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957, doi: 10.1086/323057

Dutrey, A., Guilloteau, S., & Guelin, M. 1997, A&A, 317, L55

Dutrey, A., Wakelam, V., Boehler, Y., et al. 2011, A&A, 535, A104, doi: 10.1051/0004-6361/201116931

Dutrey, A., Semenov, D., Chapillon, E., et al. 2014, Protostars and Planets VI, 317, doi: 10.2458/azu_uapress_9780816531240-ch014

Ellerbroek, L. E., Podio, L., Dougados, C., et al. 2014, A&A, 563, A87, doi: 10.1051/0004-6361/201323092

Elmegreen, B. G. 2000, ApJ, 530, 277, doi: 10.1086/308361

Espaillat, C., D'Alessio, P., Hernández, J., et al. 2010, The Astrophysical Journal, 717, 441, doi: 10.1088/0004-637x/717/1/441

Fairlamb, J. R., Oudmaijer, R. D., Mendigutía, I., Ilee, J. D., & van den Ancker, M. E. 2015, Monthly Notices of the Royal Astronomical Society, 453, 976, doi: 10.1093/mnras/stv1576

Fang, M., Pascucci, I., Edwards, S., et al. 2018, ApJ, 868, 28, doi: 10.3847/1538-4357/aae780

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067

Franz, H. B., King, P. L., & Gaillard, F. 2019, in Volatiles in the Martian Crust, ed. J. Filiberto & S. P. Schwenzer (Elsevier), 119 – 183, doi: https://doi.org/10.1016/B978-0-12-804191-8.00006-4

Fuente, A., Cernicharo, J., Agúndez, M., et al. 2010, A&A, 524, A19, doi: 10.1051/0004-6361/201014905

Fuente, A., Goicoechea, J. R., Pety, J., et al. 2017, ApJL, 851, L49, doi: 10.3847/2041-8213/aaa01b

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051

Gavino, S., Dutrey, A., Wakelam, V., et al. 2021, arXiv e-prints, arXiv:2106.05888.
https://arxiv.org/abs/2106.05888

Geballe, T. R., Baas, F., Greenberg, J. M., & Schutte, W. 1985, A&A, 146, L6

Goicoechea, J. R., Pety, J., Gerin, M., et al. 2006, A&A, 456, 565, doi: 10.1051/0004-6361:20065260
Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209,
doi: 10.1086/307195

Goodman, J., & Weare, J. 2010, Communications in
Applied Mathematics and Computational Science, Vol. 5,
No. 1, p. 65-80, 2010, 5, 65,
doi: 10.2140/camcos.2010.5.65

Grady, C. A., Devine, D., Woodgate, B., et al. 2000, ApJ,
544, 895, doi: 10.1086/317222

Grady, C. A., Hamaguchi, K., Schneider, G., et al. 2010,
ApJ, 719, 1565, doi: 10.1088/0004-637X/719/2/1565

Graedel, T. E., Langer, W. D., & Frerking, M. A. 1982,
ApJS, 48, 321, doi: 10.1086/190780

Guilloteau, S., Di Folco, E., Dutrey, A., et al. 2013, A&A,
549, A92, doi: 10.1051/0004-6361/201220298

Guilloteau, S., Dutrey, A., Piétu, V., & Boehler, Y. 2011,
A&A, 529, A105, doi: 10.1051/0004-6361/201015209

Guilloteau, S., Reboussin, L., Dutrey, A., et al. 2016, A&A,
592, A124, doi: 10.1051/0004-6361/201527088

Hartmann, L., Ballesteros-Paredes, J., & Bergin, E. A.
2001, ApJ, 562, 852, doi: 10.1086/323863

Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 263, 589,
doi: 10.1093/mnras/263.3.589

Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS,
82, 167, doi: 10.1086/191713

Hersant, F., Wakelam, V., Dutrey, A., Guilloteau, S., &
Herbst, E. 2009, A&A, 493, L49,
doi: 10.1051/0004-6361:200811082

Hirschmann, M. M. 2016, American Mineralogist, 101, 540,
doi: 10.2138/am-2016-5452

Högbom, J. A. 1974, A&AS, 15, 417

Holdship, J., Viti, S., Jimenez-Serra, I., et al. 2016,
MNRAS, 463, 802, doi: 10.1093/mnras/stw1977

Howk, J. C., Sembach, K. R., & Savage, B. D. 2006, ApJ,
637, 333, doi: 10.1086/497352

Huang, J., Öberg, K. I., & Andrews, S. M. 2016, ApJL,
823, L18, doi: 10.3847/2041-8205/823/1/L18

Huang, J., Öberg, K. I., Qi, C., et al. 2017, The
Astrophysical Journal, 835, 231,
doi: 10.3847/1538-4357/835/2/231

Huang, J., Andrews, S. M., Cleeves, L. I., et al. 2018, ApJ,
852, 122, doi: 10.3847/1538-4357/aaa1e7

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2020,
The Astrophysical Journal, 891, 48,
doi: 10.3847/1538-4357/ab711e

Huang, J., Bergin, E. A., Öberg, K. I., et al. 2021, arXiv
e-prints, arXiv:2109.06224.
https://arxiv.org/abs/2109.06224

Hunter, J. D. 2007, Computing in Science Engineering, 9,
90, doi: 10.1109/MCSE.2007.55

Ilee, J. D., Walsh, C., Booth, A. S., et al. 2021, arXiv
e-prints, arXiv:2109.06319.
https://arxiv.org/abs/2109.06319

Ingleby, L., Espaillat, C., Calvet, N., et al. 2015, The
Astrophysical Journal, 805, 149,
doi: 10.1088/0004-637x/805/2/149

Jenkins, E. B. 2009, ApJ, 700, 1299,
doi: 10.1088/0004-637X/700/2/1299

Jones, A. P., & Ysard, N. 2019, A&A, 627, A38,
doi: 10.1051/0004-6361/201935532

Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy:
Open source scientific tools for Python.
http://www.scipy.org/

Kama, M., Shorttle, O., Jermyn, A. S., et al. 2019, ApJ,
885, 114, doi: 10.3847/1538-4357/ab45f8

Kastner, J. H., Qi, C., Dickson-Vandervelde, D. A., et al.
2018, ApJ, 863, 106, doi: 10.3847/1538-4357/aacff7

Keller, L. P., Hony, S., Bradley, J. P., et al. 2002, Nature,
417, 148, doi: 10.1038/417148a

Klaassen, P. D., Juhasz, A., Mathews, G. S., et al. 2013,
A&A, 555, A73, doi: 10.1051/0004-6361/201321129

Kolesniková, L., Tercero, B., Cernicharo, J., et al. 2014,
ApJL, 784, L7, doi: 10.1088/2041-8205/784/1/L7

Kraus, A. L., & Hillenbrand, L. A. 2009, The Astrophysical
Journal, 704, 531, doi: 10.1088/0004-637x/704/1/531

Kurucz, R. L. 1993, SYNTHE spectrum synthesis programs
and line data

Laas, J. C., & Caselli, P. 2019, A&A, 624, A108,
doi: 10.1051/0004-6361/201834446

Law, C. J., Loomis, R. A., Teague, R., et al. 2021, arXiv
e-prints, arXiv:2109.06210.
https://arxiv.org/abs/2109.06210

Le Gal, R., Brady, M. T., Öberg, K. I., Roueff, E., & Le
Petit, F. 2019a, ApJ, 886, 86,
doi: 10.3847/1538-4357/ab4ad9

Le Gal, R., Herbst, E., Dufour, G., et al. 2017, A&A, 605,
A88, doi: 10.1051/0004-6361/201730980

Le Gal, R., Öberg, K. I., Loomis, R. A., Pegues, J., &
Bergner, J. B. 2019b, ApJ, 876, 72,
doi: 10.3847/1538-4357/ab1416

Lellouch, E., McGrath, M. A., & Jessup, K. L. 2007, Io's
atmosphere, ed. R. M. C. Lopes & J. R. Spencer, 231,
doi: 10.1007/978-3-540-48841-5_10

Liu, Y., Dipierro, Giovanni, Ragusa, Enrico, et al. 2019,
A&A, 622, A75, doi: 10.1051/0004-6361/201834157

Loomis, R. A., Cleeves, L. I., Öberg, K. I., et al. 2018a,
ApJ, 859, 131, doi: 10.3847/1538-4357/aac169

Loomis, R. A., Öberg, K. I., Andrews, S. M., et al. 2018b,
AJ, 155, 182, doi: 10.3847/1538-3881/aab604
—. 2020, ApJ, 893, 101, doi: 10.3847/1538-4357/ab7cc8
Macías, E., Espaillat, C. C., Ribas, Á., et al. 2018, ApJ, 865, 37, doi: 10.3847/1538-4357/aad811

Mawet, D., Absil, O., Montagnier, G., et al. 2012, A&A, 544, A131, doi: 10.1051/0004-6361/201219662

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127

Meier, R., & A'Hearn, M. F. 1997, Icarus, 125, 164, doi: 10.1006/icar.1996.5600

Mendigutía, I., Brittain, S., Eiroa, C., et al. 2013, The Astrophysical Journal, 776, 44, doi: 10.1088/0004-637x/776/1/44

Millar, T. J., & Herbst, E. 1990, A&A, 231, 466

Montesinos, B., Eiroa, C., Mora, A., & Merín, B. 2009, A&A, 495, 901, doi: 10.1051/0004-6361:200810623

Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, Journal of Molecular Structure, 742, 215, doi: 10.1016/j.molstruc.2005.01.027

Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49, doi: 10.1051/0004-6361:20010367

Navarro-Almaida, D., Le Gal, R., Fuente, A., et al. 2020, A&A, 637, A39, doi: 10.1051/0004-6361/201937180

Neufeld, D. A., Godard, B., Gerin, M., et al. 2015, A&A, 577, A49, doi: 10.1051/0004-6361/201425391

Nilsson, A., Hjalmarson, Å., Bergman, P., & Millar, T. J. 2000, A&A, 358, 257

Öberg, K. I., & Bergin, E. A. 2021, PhR, 893, 1, doi: 10.1016/j.physrep.2020.09.004

Öberg, K. I., Qi, C., Fogel, J. K. J., et al. 2011, ApJ, 734, 98, doi: 10.1088/0004-637X/734/2/98

Oberg, K. I., Guzman, V. V., Walsh, C., et al. 2021, arXiv e-prints, arXiv:2109.06268. https://arxiv.org/abs/2109.06268

Pacheco-Vázquez, S., Fuente, A., Baruteau, C., et al. 2016, A&A, 589, A60, doi: 10.1051/0004-6361/201527089

Palumbo, M. E., Geballe, T. R., & Tielens, A. G. G. M. 1997, ApJ, 479, 839, doi: 10.1086/303905

Palumbo, M. E., Tielens, A. G. G. M., & Tokunaga, A. T. 1995, ApJ, 449, 674, doi: 10.1086/176088

Phuong, N. T., Chapillon, E., Majumdar, L., et al. 2018, A&A, 616, L5, doi: 10.1051/0004-6361/201833766

Piétu, V., Dutrey, A., & Guilloteau, S. 2007, A&A, 467, 163, doi: 10.1051/0004-6361:20066537

Pinte, C., Ménard, F., Duchêne, G., et al. 2018, A&A, 609, A47, doi: 10.1051/0004-6361/201731377

Podio, L., Lefloch, B., Ceccarelli, C., Codella, C., & Bachiller, R. 2014, A&A, 565, A64, doi: 10.1051/0004-6361/201322928

Pontoppidan, K. M., Salyk, C., Banzatti, A., et al. 2019, ApJ, 874, 92, doi: 10.3847/1538-4357/ab05d8

Ranjan, S., Todd, Z. R., Sutherland, J. D., & Sasselov, D. D. 2018, Astrobiology, 18, 1023, doi: 10.1089/ast.2017.1770

Rivière-Marichalar, P., Fuente, A., Goicoechea, J. R., et al. 2019, A&A, 628, A16, doi: 10.1051/0004-6361/201935354

Rivière-Marichalar, P., Fuente, A., Le Gal, R., et al. 2020, A&A, 642, A32, doi: 10.1051/0004-6361/202038549

Ruf, A., Bouquet, A., Boduch, P., et al. 2019, ApJL, 885, L40, doi: 10.3847/2041-8213/ab4e9f

Ruffle, D. P., Hartquist, T. W., Caselli, P., & Williams, D. A. 1999, MNRAS, 306, 691, doi: 10.1046/j.1365-8711.1999.02562.x

Salyk, C., Herczeg, G. J., Brown, J. M., et al. 2013, The Astrophysical Journal, 769, 21, doi: 10.1088/0004-637x/769/1/21

Semenov, D., Favre, C., Fedele, D., et al. 2018, A&A, 617, A28, doi: 10.1051/0004-6361/201832980

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

Shirley, Y. L. 2015, PASP, 127, 299, doi: 10.1086/680342

Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, arXiv e-prints, arXiv:2109.06433. https://arxiv.org/abs/2109.06433

Simon, M., Guilloteau, S., Beck, T. L., et al. 2019, The Astrophysical Journal, 884, 42, doi: 10.3847/1538-4357/ab3e3b

Teague, R. 2019, The Journal of Open Source Software, 4, 1632, doi: 10.21105/joss.01632

Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378, doi: 10.1038/s41586-019-1642-0

Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018a, ApJ, 868, 113, doi: 10.3847/1538-4357/aae836

Teague, R., & Foreman-Mackey, D. 2018, Research Notes of the American Astronomical Society, 2, 173, doi: 10.3847/2515-5172/aae265

Teague, R., Henning, T., Guilloteau, S., et al. 2018b, ApJ, 864, 133, doi: 10.3847/1538-4357/aad80e

Teague, R., Bae, J., Aikawa, Y., et al. 2021, arXiv e-prints, arXiv:2109.06218. https://arxiv.org/abs/2109.06218

Tieftrunk, A., Pineau des Forets, G., Schilke, P., & Walmsley, C. M. 1994, A&A, 289, 579

van der Plas, G., Casassus, S., Ménard, F., et al. 2014, ApJL, 792, L25, doi: 10.1088/2041-8205/792/2/L25

van der Tak, F. F. S., Lique, F., Faure, A., Black, J. H., & van Dishoeck, E. F. 2020, Atoms, 8, 15, doi: 10.3390/atoms8020015
van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37

van der Walt, S., Schönberger, J. L., Nunez-Iglesias, J., et al. 2014, ArXiv e-prints. https://arxiv.org/abs/1407.6245

van Dishoeck, E. F., & Blake, G. A. 1998, ARA&A, 36, 317, doi: 10.1146/annurev.astro.36.1.317

van Dishoeck, E. F., Kristensen, L. E., Mottram, J. C., et al. 2021, A&A, 648, A24, doi: 10.1051/0004-6361/202039084

Vastel, C., Quénard, D., Le Gal, R., et al. 2018, MNRAS, 478, 5514, doi: 10.1093/mnras/sty1336

Vidal, T. H. G., Loison, J.-C., Jaziri, A. Y., et al. 2017, MNRAS, 469, 435, doi: 10.1093/mnras/stx828

Wagenblast, R., & Hartquist, T. W. 1989, MNRAS, 237, 1019, doi: 10.1093/mnras/237.4.1019

Wakelam, V., Caselli, P., Ceccarelli, C., Herbst, E., & Castets, A. 2004, A&A, 422, 159, doi: 10.1051/0004-6361:20047186

Wakelam, V., Chapillon, E., Dutrey, A., et al. 2019, MNRAS, 484, 1563, doi: 10.1093/mnras/stz068

Wakelam, V., & Herbst, E. 2008, ApJ, 680, 371, doi: 10.1086/587734

Wakelam, V., Loison, J. C., Mereau, R., & Ruaud, M. 2017, Molecular Astrophysics, 6, 22, doi: 10.1016/j.molap.2017.01.002

Wakelam, V., Ruaud, M., Hersant, F., et al. 2016, A&A, 594, A35, doi: 10.1051/0004-6361/201628748

Whittet, D. C. B. 2010, ApJ, 710, 1009, doi: 10.1088/0004-637X/710/2/1009

Wiesenfeld, L., & Faure, A. 2013, MNRAS, 432, 2573, doi: 10.1093/mnras/stt616

Yen, H.-W., Koch, P. M., Liu, H. B., et al. 2016, ApJ, 832, 204, doi: 10.3847/0004-637X/832/2/204