Sulfur chemistry in protoplanetary disks: CS and H2CS

TL;DR

This study investigates sulfur chemistry in protoplanetary disks using ALMA observations of sulfur molecules like CS and H2CS, revealing diverse distributions and challenging existing chemical models, with implications for planetary habitability.

Contribution

It provides new observational data of sulfur molecules in PPDs and tests chemical models, highlighting discrepancies and the need for further high-resolution studies.

Findings

CS shows diverse radial distributions including gaps and rings.

First detections of $^{13}$CS, C$^{34}$S, and H$_2$CS in PPDs.

Model predictions match observations within a factor of a few for CS, but underpredict H$_2$CS by 1-2 orders of magnitude.

Abstract

The nature and abundance of sulfur chemistry in protoplanetary disks (PPDs) may impact the sulfur inventory on young planets and therefore their habitability. PPDs also present an interesting test bed for sulfur chemistry models, since each disk present a diverse set of environments. In this context, we present new sulfur molecule observations in PPDs, and new S-disk chemistry models. With ALMA we observed the CS 5-4 rotational transition toward five PPDs (DM Tau, DO Tau, CI Tau, LkCa 15, MWC 480), and the CS 6-5 transition toward three PPDs (LkCa 15, MWC 480 and V4046 Sgr). Across this sample, CS displays a range of radial distributions, from centrally peaked, to gaps and rings. We also present the first detection in PPDs of $^{13}$CS 6-5 (LkCa 15 and MWC 480), C$^{34}$S 6-5 (LkCa 15), and H$_2$CS $8_{17}-7_{16}$, $9_{19}-8_{18}$ and $9_{18}-8_{17}$ (MWC 480) transitions. Using LTE models to constrain column densities and excitation temperatures, we find that either $^{13}$C and $^{34}$S are enhanced in CS, or CS is optically thick despite its relatively low brightness temperature. Additional lines and higher spatial resolution observations are needed to distinguish between these scenarios. Assuming CS is optically thin, CS column density model predictions reproduce the observations within a factor of a few for both MWC 480 and LkCa 15. However, the model underpredicts H$_2$CS by 1-2 orders of magnitude. Finally, comparing the H$_2$CS/CS ratio observed toward the MWC~480 disk and toward different ISM sources, we find the closest match with prestellar cores.

Figures

Tables

Table 1. Table 1: Disk characteristics.

Source	Stellar Type	R.A.(a)	Dec.(a)	Dist.(a)	$L_{\star }$	$M_{\star }$	Disk inc.b	Disk P.A.b	Age	VLSR
		(J2000)	(J2000)	(pc)	($L_{\odot }$)	($M_{\odot }$)	(∘)	(∘)	(Myr)	(km/s)
DM Tau	M1-M2[1,2]	04:33:48.7	18:10:10.0	145	0.2[3]	0.5-0.6[3,4]	54.4	-25.8	3.5-8[3,5]	6.05[6]
DO Tau	M0-M1[1,2]	04:38:28.6	26:10:49.5	139	1.4[3]	0.5[3]	30.3	5.2	0.8[3]	6.5[6]
CI Tau	K4-K7[1,2]	04:33:52.0	22:50:30.1	159	0.9[3]	0.8[3,5]	50.3	194.5	1.5-2.5[3,5]	5.77[6]
MWC 480	A1-A3/4 [7,2]	04:58:46.3	29:50:37.0	162	19-24[8,3]	1.7-2.3[3,4,8]	35.3	-37.0	6-7.1[3,4,8]	5.1[9]
LkCa 15	K3-K5[1,2]	04:39:17.8	22:21:03.4	159	0.8[3]	1.0[3,4]	49.7	61.4	3-5[3,4,5]	6.3[9]
V4046 Sgr(c)	K5, K7[10]	18:14:10.5	-32:47:35.3	72	0.5, 0.3[10]	0.9, 0.9[11,12]	40.3	255.5	4-13[10,13]	2.9[9]

Table 2. Table 2: List of Observations (data from CDMS (a) )

Species	Transition	Frequency	$E_{up}$	S${}_{ij}{}^{}\mu _{}^2$	Source	RMSchan	Restored Beam		$R_{\sigma }^{\left(b\right)},R_{\max }$	$S_{\nu }{\Delta }_v$($R_{\sigma }$)	$S_{\nu }{\Delta }_v$($R_{\max }$)
		(GHz)	(K)	(D2)		(mJy/beam)	(${}_{}{}^{\prime \prime }\times _{}^{\prime \prime }$)	(∘)	(AU)	(mJy km/s)	(mJy km/s)
12CS	$5-4$	244.93564	35.3	19.3	MWC 480	4.0	$0.45\times 0.71$	-12.9	150, 400	$153\pm 6$	$536\pm 19$
							$1.05\times 1.26$	-15.9		$119\pm 6$	$509\pm 25$
					LkCa 15	2.8	$0.47\times 0.57$	20.8	200, 600	$301\pm 5$	$1050\pm 14$
							$1.01\times 1.15$	37.6		$233\pm 9$	$1065\pm 28$
					CI Tau	2.8	$0.46\times 0.57$	16.0	200, 600	$182\pm 6$	$756\pm 17$
					DO Tau	3.1	$0.45\times 0.59$	16.7	100, 200	$60\pm 3$	$159\pm 6$
					DM Tau	3.1	$0.47\times 0.54$	20.4	200, 500	$159\pm 5$	$388\pm 9$
	$6-5$	293.91224	49.4	23.1	MWC 480	8.8	$0.36\times 0.78$	-26.5	150, 400	$185\pm 6$	$688\pm 13$
							$0.98\times 1.28$	-36.6		$143\pm 3$	$628\pm 12$
					LkCa 15	8.3	$0.38\times 0.70$	-32.0	200, 600	$415\pm 10$	$1282\pm 25$
							$0.96\times 1.15$	-47.1		$319\pm 5$	$1339\pm 17$
					V4046 Sgr				–, 200	–	$886\pm 11$
13CS	$6-5$	277.45548	49.6	23.1	MWC 480	11.5	$1.05\times 1.28$	-10.7	150, 400	$\lesssim 4$	$\lesssim 6$
					LkCa 15	13.2	$1.02\times 1.15$	-43.3	200, 600	$17\pm 4$	$71\pm 17$
C34S	$6-5$	289.20923	38.2	22.3	MWC 480	13.5	$1.01\times 1.22$	-11.1	150, 400	$\lesssim 6$	$21\pm 12$
					LkCa 15	15.2	$0.98\times 1.11$	-45.7	200, 600	$24\pm 4$	$87\pm 17$
H2CS	$8_{17}-7_{16}$	278.88766	73.4	64.1	MWC 480	13.8	$1.04\times 1.28$	-10.7	–, 300	–	$52\pm 9$
					LkCa 15	16.5	$1.02\times 1.15$	-43.3	–, 320	–	$\lesssim 16$
	$9_{19}-8_{18}$	304.30770	86.2	72.3	MWC 480	8.8	$0.62\times 0.94$	-1.0	–, 300	–	$49\pm 10$
					LkCa 15	8.8	$0.63\times 0.82$	-8.4	–, 320	–	$\lesssim 22$
	$9_{18}-8_{17}$	313.71665	88.5	72.3	MWC 480	15.8	$0.53\times 0.95$	-24.8	–, 300	–	$51\pm 20$
					LkCa 15	18.3	$0.53\times 0.86$	-33.2	–, 320	–	$\lesssim 50$

Table 3. Table 3: Column densities and excitation temperatures derived from the methods presented in Sect.  3 when applied

Species	Source	Radius	$N_{\mathrm{tot}}$	$T_{\mathrm{ex}}$	Method
		(AU)	(${\mathrm{cm}}^{-2}$)	(K)
${}_{}{}^{12}\text{ce}CS$	MWC 480	150	$\left(7.3\pm 0.8\right)\times {10}^{12}$	$22.8\pm 2.5$	Rotational diagram
	LkCa 15	200	$\left(1.2\pm 0.1\right)\times {10}^{13}$	$26.8\pm 2.3$	Rotational diagram
			$\left[1.5\times {10}^{13}-5\times {10}^{14}\right]$	$\sim 7$	Isotopologue ratio analysis
\ceH2CS	MWC 480	300	${3.0}_{-0.4}^{+2.7}\times {10}^{12}$	${40.5}_{-20.2}^{+38.4}$	Rotational diagram
	LkCa 15	320	$\lesssim 2\times {10}^{12}$	$40$	Using Eqs. A1 and A2 and a fixed $T_{\mathrm{ex}}$

Table 4. Table 4: Physical parameters used for our disk models

Parameters	MWC 480a	LkCa 15b
Stellar mass: $M_{\star }$ ($M_{\odot }$)	1.8	1.0
Disk mass: $M_d$ ($M_{\odot }$)	0.18	0.03
Characteristic radius: $R_c$ (AU)	100	100
Outer cut-off radius: $R_{\mathrm{out}}$ (AU)	500	550
Density power-law index: $\gamma$	1.5	1.5
Midplane temperature at $R_c$c: $T_{100AU}$(K)	30	15
Atmosphere temperature at $R_c$: $T_{100AU}$(K)	48	20
Surface density at $R_c$:	5.7	0.9
Temperature power-law index: $q$	0.5	0.5
Vertical temperature gradient index: $\beta$	2	2
UV Flux of reference: $f_{\mathrm{UV},R_c}$ (in Draine (1978)’s units)	8500d	2550e

Table 5. Table 5: Initial Elemental Abundances

Species	$n_i/n_{\text{H}}$	Reference
H2	0.5
He	9.0$\times$10-2	1
C+	1.7$\times$10-4	2
N	6.2$\times$10-5	2
O	2.4$\times$10-4	3
S+	8.0$\times$10-8	4
Si+	8.0$\times$10-9	4
Fe+	3.0$\times$10-9	4
Na+	2.0$\times$10-9	4
Mg+	7.0$\times$10-9	4
P+	2.0$\times$10-10	4
Cl+	1.0$\times$10-9	4
F+	6.7$\times$10-9	5

Equations

$$T_{\mathrm{ex}}(\mathrm{K})=\mathcal{U}(3,300)$$

$$N_{\mathrm{tot}}(\mathrm{c}{\mathrm{m}}^{-2})=\mathcal{U}(10^7,10^{20}).$$

$$\frac{T_B({}^{13}\mathrm{CS})}{T_B({}^{12}\mathrm{CS})}=\frac{1-e^{-{\tau }_{{}^{13}\mathrm{CS}}}}{1-e^{-{\tau }_{{}^{12}\mathrm{CS}}}},$$

$$\frac{T_B({}^{13}\mathrm{CS})}{T_B({}^{12}\mathrm{CS})}=\frac{{\tau }_{{}^{13}\mathrm{CS}}}{{\tau }_{{}^{12}\mathrm{CS}}}\approx R.$$

$$T_{\mathrm{B}}\Delta v=\frac{c^2{S}_{\nu }\Delta v}{2k_{\mathrm{B}}{\nu }^2\Omega }.$$

$$\frac{T_B({}^{12}\mathrm{CS})}{T_B({}^{13}\mathrm{CS})}\simeq 17\ll {\left(\frac{{}^{12}\text{ce}C}{{}^{13}\text{ce}C}\right)}_{\mathrm{LISM}}.$$

$$T_{ex}({}^{12}\mathrm{CS})=\frac{\mathrm{h}\nu /{\mathrm{k}}_{\mathrm{B}}}{\ln \left(1+\frac{\mathrm{h}\nu /{\mathrm{k}}_{\mathrm{B}}}{{\mathrm{T}}_{\mathrm{R}}({}^{12}\mathrm{CS})/f+{\mathrm{J}}_{\nu }({\mathrm{T}}_{\mathrm{bg}})}\right)}.$$

$$\frac{T_B({}^{13}\mathrm{CS})}{T_B({}^{12}\mathrm{CS})}\approx {\tau }_{{}^{13}\mathrm{CS}}.$$

$$\begin{array}{r}{N}_{\mathrm{tot}}=\frac{3h}{8{\pi }^{3}S{\mu }^2}\frac{Q_{rot}(T_{\mathrm{ex}})}{g_u}\frac{e^{E_u/k_B{T}_{ex}}}{e^{h\nu /k_B{T}_{ex}}-1}\\ \times \frac{1}{J_{\nu }(T_{ex})-J_{\nu }(T_{bg})}\int T_{R}dv{C}_{\tau }.\end{array}$$

$$T_{\mathrm{ex}}(\mathrm{K})=\mathcal{U}(10,100),$$

$$N_{\mathrm{tot}}(\mathrm{c}{\mathrm{m}}^{-2})=\mathcal{U}(10^7,10^{20}).$$

$$\small T(z)=\left\{\begin{array}[]{ll}T_{\rm{mid}}+(T_{\rm{atm}}-T_{\rm{mid}})\left[\sin\left(\frac{\pi z}{2z_{q}}\right)\right]^{2\delta}&\mbox{if}\,z<z_{q}\\ T_{\rm{atm}}&\mbox{if}\,z\geq z_{q},\end{array}\right.$$

$$T_{\mathrm{mid}}=T_{\mathrm{mid},R_c}{\left(\frac{r}{R_{\mathrm{c}}}\right)}^{-q},$$

$$T_{\mathrm{atm}}=T_{\mathrm{atm},R_c}{\left(\frac{r}{R_{\mathrm{c}}}\right)}^{-q},$$

$$H=\sqrt{\frac{k_{\mathrm{B}}{T}_{\mathrm{mid}}{r}^3}{\mu m_{\mathrm{H}}G{M}_{\star }}},$$

$$T_{\mathrm{mid}}(r)\approx {\left(\frac{\phi L_{\star }}{8\pi r^2{\sigma }_{\mathrm{SB}}}\right)}^{1/4},$$

$$\frac{\delta \ln \rho (z)}{\delta z}=-\left[\frac{G{M}_{\star }z}{r^3}\frac{\mu m_{\mathrm{H}}}{k_{B}T}+\frac{\delta \ln T}{\delta z}\right],$$

$$\frac{\delta \rho (z)}{\rho (z)}=-\left[\frac{G{M}_{\star }z}{r^3}\frac{\mu m_{\mathrm{H}}}{k_{B}T{}_{\mathrm{mid}}}\right]\delta z,$$

$$\rho (z)={\rho }_o\exp \left(-\frac{z^2}{2H^2}\right),$$

$$\Sigma ={\rho }_0{\int }_{-\infty }^{+\infty }\exp \left(-\frac{z^2}{2H^2}\right)=\sqrt{2\pi }{\rho }_{0}H.$$

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

$$2\pi r\frac{\delta \Sigma }{\delta t}=-\frac{\dot{M}}{\delta r},$$

$$M_{\mathrm{disk}}={\int }_{R_{\mathrm{int}}}^{R_{\mathrm{out}}}2\pi r\Sigma (r)dr.$$

$${\Sigma }_R=\frac{M_{\mathrm{disk}}{R}^{-3/2}}{4\pi \sqrt{R_{\mathrm{out}}}},$$

$$f_{\mathrm{UV}}=\frac{f_{\mathrm{UV},{\mathrm{R}}_{\mathrm{c}}}/2}{{\left(\frac{r}{R_c}\right)}^2+{\left(\frac{4\mathrm{H}}{R_c}\right)}^2}.$$

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

$$\frac{N_u}{g_u}=\frac{N_{\mathrm{tot}}}{Q_{\mathrm{rot}}(T)}{e}^{-E_u/k_{B}T},$$

$$Q_{\mathrm{rot}}(T)\approx \frac{k_{B}T}{h{B}_0}+\frac{1}{3},$$

$$\ln \frac{N_u}{g_u}=\ln N_{\mathrm{tot}}-\ln Q_{\mathrm{rot}}(T)-\ln C_{\tau }-\frac{E_u}{k_{B}T}.$$

$$\int {\tau }_{\nu }\mathrm{d}v=\frac{N_u{A}_{ul}{c}^3}{8\pi {\nu }^3}(e^{h\nu /k_B{T}_{\mathrm{ex}}}-1).$$

$$\int {\tau }_{\nu }\mathrm{d}v={\tau }_{\nu }\Delta v_{\mathrm{FWHM}}\sqrt{\frac{\pi }{4\ln 2}},$$

$${\tau }_{\nu }=\sqrt{\frac{4\ln 2}{\pi }}\frac{N_u{A}_{ul}{c}^3}{\Delta v_{\mathrm{FWHM}}8\pi {\nu }^3}(e^{h\nu /k_B{T}_{\mathrm{ex}}}-1).$$

Full text

astrochemistry — ISM: molecules — methods: observational — methods: numerical — protoplanetary disks — techniques: interferometric

††journal: ApJ††software: CASA

1 Introduction

Planets form in dust and gas-rich disks, protoplanetary disks (PPDs), around young stars. Disk chemical structures thus link interstellar chemistry to the chemistry of nascent planets. The chemical composition of these disks is set by a combination of inheritance from their parent molecular cloud and in situ processes (e.g. Aikawa & Herbst, 1999; Willacy, 2007; Willacy & Woods, 2009; Dutrey et al., 2014; Sakai et al., 2014; Cleeves et al., 2014; Öberg et al., 2015b). Much of this chemical composition is still poorly constrained, however, and in the case of sulfur, its chemistry and main reservoirs are yet largely unknown (e.g. Semenov et al., 2018; Phuong et al., 2018).

Sulfur chemistry is not just poorly understood in disks, but also in many other circumstellar and interstellar environments. In the diffuse interstellar medium (ISM) and in photon-dominated regions (PDRs) the observed sulfur abundance is close to the cosmic value (e.g. Howk et al., 2006; Neufeld et al., 2015; Goicoechea et al., 2006), whereas in dense molecular gas it is found highly depleted: only 0.1% of the sulfur cosmic abundance is observed in the gas phase (Tieftrunk et al., 1994), implying a depletion factor of three orders of magnitude (e.g. Wakelam et al., 2004; Vastel et al., 2018). This level of depletion suggests that most sulfur is locked into icy mantles coating interstellar dust grains (Millar & Herbst, 1990; Ruffle et al., 1999; Vidal et al., 2017; Laas & Caselli, 2019). However, OCS and tentatively \ceSO2 are yet the sole molecules detected in icy grain mantles toward high-mass protostars (e.g. Geballe et al., 1985; Palumbo et al., 1997; Zasowski et al., 2009) and their abundances ($\sim{\,\rm 10^{-7}}$) make up less than 4% of the sulfur cosmic abundance. Icy \ceH2S, the most obvious sulfur ice chemistry product from successive hydrogenations of atomic sulfur in solid phase, has not been detected yet (Smith, 1991; Boogert et al., 2015).

While most of the sulfur is hidden from us in star forming regions, recent spectral line surveys have increased the number of known interstellar sulfur molecules to a dozen of species in prestellar cores (Vastel et al., 2018), protostellar envelopes (Drozdovskaya et al., 2018) and PDRs (Riviere-Marichalar et al. submitted). Together with new models, these new detections have advanced our understanding of the ISM gas-phase S-chemistry. In particular, we now know that a rich interstellar sulfur chemistry is present in various ISM environments. The main S-reservoirs have not yet been identified, however, and there is still much theoretical work left to do to identify the chemical pathways that produce the observed distributions of sulfur species.

It is also unclear how much of the interstellar sulfur molecules survive during star formation and become incorporated into planet-forming disks, or whether the sulfur chemistry in such disks is largely reset. Sulfur species are present in the remnants of our own PPD, comets (e.g. Boissier et al., 2007; Le Roy et al., 2015) which, could actually even be relics from its parent ISM molecular cloud according to the Rosetta analysis of the 67P/Churyumov-Gerasimenko cometary ices (e.g. Altwegg et al., 2016). In theory, the cometary sulfur molecule abundance patterns and isotopic ratios could be used to constrain whether they are of interstellar or disk origin, but this requires a much better understanding of disk sulfur chemistry than is currently available. A first motivation for pursuing disk sulfur studies is therefore to better constrain the chemical origins of the Solar System. A second motivation is to enable predictions of what kind of sulfur molecules become incorporated into young planets, since this may affect their hospitality to origins of life (Ranjan et al., 2018).

A third motivation is that disks provide a diverse set of chemical environments and could therefore be used to test and develop a holistic sulfur chemistry model. In particular, protoplanetary disks are vertically stratified into disk atmospheres, warm molecular layers and cold midplanes that are analogous to PDRs, molecular clouds and dense cloud cores, respectively. Dependent on which region that dominates the observed emission, disks may express PDR, cloud or core-like sulfur chemistries. Which one is observed may vary between disks around T Tauri and Herbig Ae stars, and perhaps even across individual disks.

Sulfur bearing molecules have been detected in PPDs, but the detections are scarce compared to the earlier stages of star formation described above. So far, only three S-bearing species detections have been reported in the literature. CS is detected in dozens of disks (e.g. Dutrey et al., 1997; Fuente et al., 2010; Guilloteau et al., 2016; Dutrey et al., 2017; Teague et al., 2018; Phuong et al., 2018), while SO is only found toward a few young disks presenting active accretion signs, including the Herbig A0 star AB Aur and Class I sources (e.g. Fuente et al., 2010; Guilloteau et al., 2013, 2016; Pacheco-Vázquez et al., 2016; Sakai et al., 2016; Booth et al., 2018). Finally, very recently, the first detection of \ceH2S was reported toward the T Tauri star GG Tau A (Phuong et al., 2018).

Here we present a holistic analysis of new observations of S-species in disks using the Atacama Large Millimeter/submillimeter Array (ALMA). The survey includes the CS $5-4$ rotational transition observed toward five Taurus PPDs (DM Tau, DO Tau, CI Tau, LkCa 15, MWC 480), the CS $6-5$ transition observed toward two of these disks (LkCa 15 and MWC 480) and V4046 Sgr, and the first detection in PPDs of the species 13CS, C34S, and \ceH2CS in their $6-5$, and $8_{17}-7_{16}$, $9_{19}-8_{18}$ and $9_{18}-8_{17}$ rotational transitions, respectively. The latter three are observed toward LkCa 15 and/or MWC 480 (see also Loomis et al. in prep.). These observations are compared to two dimensional disk chemistry models.

The outline of the paper is the following. The observational details are presented in §2 and the observational results in §3, including derived column densities and excitation temperatures for CS and \ceH2CS toward LkCa 15 and MWC 480. To better understand the observed PPD S-chemistry, we ran two disk chemistry models tuned to MWC 480 and LkCa 15, which we present in §4. The observational and modeling results are compared and discussed in §5 and the conclusions presented in §6.

2 Observations

2.1 Sample

Our disk sample is composed of six PPDs. Five of them – surrounding the DM Tau, DO Tau, CI Tau, MWC 480, and LkCa 15 stars – are situated in the nearby Taurus star forming region (from 139 to 169 pc, (Gaia Collaboration, 2018), see Table 1). Their ages range from $\sim 1$ to $\sim 10$ Myr (see Table 1), and host large ($>200$ au) gas-rich disks (e.g. Chiang et al., 2001; Piétu et al., 2007; Öberg et al., 2010; Isella et al., 2012; Guilloteau et al., 2016; Huang et al., 2017).

The sixth PPD is orbiting V4046 Sgr, an isolated binary system located in the $~{}23\pm 3$Myr old $\beta$-Pictoris Moving Group (Mamajek & Bell, 2014). Among the six disks, three are known to be transitional disks, namely DM Tau, LkCa 15 and V4046 Sgr, meaning that they present large gaps and/or inner cavities in their submillimeter/millimeter dust continuum emission (Andrews et al., 2009; Piétu et al., 2006; Rosenfeld et al., 2013).

2.2 Observational details

The observed species transitions, their frequencies and spectroscopic parameters are listed in Table 2.

The CS $5-4$ transition was observed in the five Taurus disks (ALMA project code 2016.1.00627.S, PI: K. Öberg), in Band 6 during ALMA cycle 4, with a spectral resolution of 61 kHz. These observations were performed in two execution blocks on December 1, 2016 with an angular resolution of $\sim 0.5^{\prime\prime}$. The total on-source integration times were 22-24 min. 44 antennas were included and covered baselines from 15 to 704 m. The first block used the source J0510+1800 as its band-pass calibrator, while the second block used J0237+2848. Both blocks used the J0423-0120 and J0426+2327 source as a flux calibrator and phase calibrator, respectively.

The three \ceH2CS $8_{17}-7_{16}$, $9_{19}-8_{18}$ and $9_{18}-8_{17}$ rotational transitions, CS $6-5$ and corresponding 13CS and C34S isotopologues were observed in Band 7 during ALMA cycles 3 & 4 (project code 2013.1.01070.S, PI: K. Öberg), toward two of the Taurus disks: the Herbig Ae disk MWC 480, and the T Tauri disk LkCa 15. The spectral resolution for these lines was $\sim 975$ kHz, i.e. corresponding to a velocity resolution of $\sim 1$ km/s.

The 12CS $6-5$ transition was observed in three execution blocks on January 17, 2016, April 23, 2016 and December 12, 2016. For the first execution block 31 antennas were included and covered baselines from 15 to 331 m. 36 antennas were included for the two other execution blocks, covering baselines from 15 to 463 m and 15 to 650 m, respectively. The first and third blocks used the source J0510+1800 as band-pass and flux calibrators, and the source J0438+3004 as phase calibrator. The second block of execution used the source J0238+1636 as band-pass calibrator, the source J0433+2905 as phase calibrator, and the source J0510+1800 as flux calibrator. The total on-source integration time were 17.6, 12.6 and 20.7 min., respectively.

The 13CS and C34S $6-5$, and \ceH2CS $8_{17}-7_{16}$ rotational transitions were observed in a single execution block, on January 17, 2016, with the same calibrator sources as those used for the first and second execution blocks used for 12CS $6-5$ (see above) but with 36 antennas and a total on-source integration time of $\sim 19.2$ min.

The \ceH2CS $9_{19}-8_{18}$ rotational transition was observed in two execution blocks, on December 17 and 18, 2016, with the same calibrator sources as those used for the first and second execution blocks used for 12CS $6-5$ (see above). 40 and 42 antennas were used with baseline coverages from 15 to 460 m and 15 to 492 m, respectively. The total on-source integration times were $\sim 13.1$ min for each block.

The \ceH2CS $9_{18}-8_{17}$ rotational transition was observed in two execution blocks, on December 13 and 14, 2016, with the same calibrator sources as those used for the first and second execution blocks used for 12CS $6-5$ (see above). 36 and 39 antennas were used with baseline coverages of 15 to 650 m and 15 to 460 m, respectively. The total on-source integration times were $\sim 12.6$ min for each block.

More details on the ALMA Band 7 observations presented above can be found in Loomis et al. in prep.

Finally, the CS $6-5$ transition was observed toward V4046 Sgr as part of the ALMA project 2016.1.01046.S (PI: R. Loomis). These observations were performed in a single execution block on May 06, 2018 with an angular resolution of $0.54^{\prime\prime}$ and a spectral resolution of 0.1 km/s. The total on-source integration times was 29 min. 44 antennas were included and covered baselines from 15 to 500 m. The quasar J1924-2914 was used as a band-pass and flux calibrator, and J1802-3940 was used as a phase calibrator.

2.3 Data reduction

Data calibration was initially performed by the ALMA/NAASC staff. We then use the Common Astronomy Software Application package (CASA) version 4.7.2 (McMullin et al., 2007) to reduce the data. Self-calibrations were performed using the continuum emission from each spectral window containing the targeted line.

After continuum subtraction with the CASA uvconsub function, the data were CLEANed (Högbom, 1974) using a $3\sigma$ noise threshold and Briggs weighting with a robustness parameter of 0.5 for all the lines. The one exception was the high signal to noise ratio CS $6-5$ line observed toward V4046 Sgr, for which a value of 0 was adopted in order to improve the angular resolution. The RMS per channel are listed in Table 2. The masks used for the CLEANing procedure were drawn manually to follow the shape of the emission in each velocity channel.

We applied a Keplerian mask (e.g. Salinas et al., 2017, and see details in Appendix of Pegues et al. in prep for our methodology) to the CLEANed data cubes to improve the signal-to-noise of extracted intensity maps, spectra, and integrated flux for each line. Briefly, the Keplerian mask selects the disk area in which the line emits within each channel, given the beam size, velocity resolution, stellar mass, and disk geometry, and assuming that the disk is flat. In theory the flat disk assumption could result in emission that is coming from the flaring layers being cut out. However we pad the mask sufficiently to avoid this and have validated it against CO emission (Pegues et al. in prep.), which is originating in at least as high layers as our molecules due to high optical depths. The Keplerian masks calculated for the different lines studied here are depicted on the channel maps in Appendix Fig. 17-33. The disk and star parameters used to build the Keplerian masks are listed Table 1.

3 Results

3.1 CS fluxes and emission morphologies

Figure 1 shows the spatially and spectrally resolved CS $5-4$ line detections toward the five Taurus disks of our sample, together with their corresponding $\sim 1.3$ mm continuum emission maps. The angular resolution is $\sim 0.5^{\prime\prime}$, which corresponds to $\sim 80$ au. The spectra of the CS $5-4$ are also depicted in Fig. 1. For four of our targets, a typical double-peaked profile tracing the Keplerian rotation of the disk, while the double-peaked signature is not obvious toward the faintest, DO Tau, PPD. This could be due to spectral resolution issue, but DO Tau is supposed to be very young so there may also be a non-Keplerian component. Further DO Tau studies would help in better constrain its physical structure. Figure 2 shows CS $6-5$ line and the $\sim 1$ mm continuum emission maps toward V4046 Sgr at an angular resolution of $\sim 0.5^{\prime\prime}$, which corresponds to $\sim 32$ au (see source distance in Table 1). As shown in the last panel of the figure, the spectra presents typical double-peaked profile corresponding to Keplerian rotation.

Across the sample the CS radial emission profile is centrally peaked in 5/6 disks; the only exception is the disk around Herbig Ae star MWC 480, where a small emission dip is seen at the disk center. In the outer disk the CS emission morphology appears related to the dust emission outer edge. In the radial profiles, depicted on the second and third row panels in Fig. 1, there is a break in the CS slope close to the continuum edge in all disks. In the V4046 Sgr disk, the relationship between the continuum disk and the CS morphology is even more pronounced; as there is a ring of CS gas at $\sim 100$ au, at the edge of the dust continuum (see Fig. 2).

The spectral profiles also encode information about the CS distributions in the observed disks. All disk spectra present substantial wings. This is again most pronounced toward V4046 Sgr, where the spectral line looks like a superposition of two spectra, a narrow Keplerian profile that carries most of the flux and traces emission from the outer disk ring, and a broad fainter component that traces emission from the inner disk.

The disk-integrated flux densities, $S_{\nu}\Delta v$, derived from this work are listed Table 2. We report integrated fluxes out to two different radii for each source, first out to the radius, $R_{\sigma}$, where the CS flux is $\sim 33$% of the maximum CS emission peak (thus focusing on the part of the disk where most of the CS emission originates when deriving a characteristic CS column density, see Figs. 1 and 2), and second out to the furthest radius, $R_{max}$, where CS emission is detected.

3.2 CS 6-5 and isotopologues in MWC 480 and LkCa 15 disks

Toward two disks in the sample, MWC 480 and LkCa 15, three additional CS lines were observed (Fig. 3). The 12CS $6-5$ line is strongly detected in both disks. The peak fluxes are comparable (within 20%), while the integrated fluxes are somewhat higher compared to the $5-4$ lines (Table 2). The 12CS $6-5$ emission radial morphologies follow those found for 12CS $5-4$ in the respective disk (Fig. 4), but the breaks observed at the edge of the continuum are more pronounced.

Figures 3 and 4 show that the $6-5$ line of the 13CS and C34S isotopologues are detected toward LkCa 15, and tentatively toward MWC 480. The stronger isotopologue emission toward LkCa 15 compared to MWC 480 is consistent with the fact that the former also presents stronger main isotopologue line emission.

The disk integrated emission of CS, 13CS and C34S $6-5$ in the two disks are reported in Table 2. First, whether integrating out to 150/200 au, the radius out to which isotopologue emission is clearly detected, or across the full CS disk, LkCa 15 consistently presents about a factor of two higher flux densities. Second, toward LkCa 15, the flux ratio between 13CS and the main isotopologue is quite high, which is discussed in detail in §3.3.2.

3.3 CS column densities and excitation temperatures

Emission line fluxes encode information about the molecular column density and excitation conditions. For the two disks we have observed in two 12CS lines, we can simultaneously determine the excitation temperature and column density using a rotational diagram analysis (e.g. Linke et al., 1979; Blake et al., 1987; Goldsmith & Langer, 1999). This approach implies the simplifying assumptions that the lines are optically thin and at local thermal equilibrium (LTE), and are pursued in § 3.3.1. For one disk, LkCa 15, we have additional information in the form of CS isotopologue detections, and we use these lines to evaluate the CS column density and excitation temperature in the other limiting case where CS main isotopologue lines are optically thick, in § 3.3.2. Finally, in § 3.3.3, we use the constraints on the CS rotational temperature from § 3.3.1 to estimate the CS column density in all of our disk sample.

3.3.1 Rotational diagram analysis

The first assumption for rotational diagram analysis is that the lines are in LTE. Lines can be assumed to be in LTE when their critical densities are surpassed by the gas densities in their emitting regions. The critical densities of the CS $5-4$ and $6-5$ transitions for temperatures of 20–50 K are $\sim[1.7-1.3]{\,\times 10^{6}}$ and $\sim[2.9-2.2]{\,\times 10^{6}}{\,\rm cm^{-3}}$ Shirley (2015), respectively. These are lower than expected in molecular disk layers (see e.g. Fig. 11) and LTE thus seems like a reasonable approximation. Under LTE, the population of an energy level $u$ can be approximated as $T_{{\rm ex}}=T_{{\rm kin}}$. The second assumption for rotational diagram analysis is that the lines are optically thin. This assumption can be checked by inspecting the line brightness temperature, which will be substantially below the ambient temperature for optically thin lines. Fig. 1 and 4 show that CS brightness temperatures are below 7 and 9 K, for the MWC 480 and LkCa 15 disks respectively. Thermalized lines usually have brightness temperatures which are $\sim 0.6-0.8*T_{{\rm kin}}$ (unless optical depths are extreme), meaning that a 9 K brightness temperature line could partly be optically thick if $T_{{\rm kin}}\sim 11-15$ K. But, in particular for MWC 480, the lowest expected gas temperature in CS emitting regions is $\sim 20-30$ K (e.g. Teague et al. (2018); Semenov et al. (2018), and § 4), which suggests that the lines are probably not optically thick.

Following the LTE and optically thin assumptions, the total column densities and rotational temperatures of CS can be derived from the disk-integrated flux densities, as described in appendix A. We then used the likelihood function $\mathcal{L}(N_{u},T_{{\rm ex}})$ described in appendix A with the Python implementation emcee (Foreman-Mackey et al., 2013) of the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) (Goodman & Weare, 2010) to fit the data and compute posterior probability distributions for $T_{{\rm ex}}$ and $N_{u}$ describing the range of total column densities, $N_{\rm{tot}}$, and excitation temperatures, $T_{{\rm ex}}$, consistent with the data. The following uninformative priors were assumed:

[TABLE]

Random draws and results from the posterior distributions are depicted in gray in Fig. 5, along with the least square fit results, shown in blue. For the MWC 480 disk, the results converged toward {$T_{{\rm ex}}\simeq 22.8\pm 2.5$ K, $N_{\rm{tot}}\simeq 7.3\pm 0.9{\,\times 10^{12}}{\,\rm cm^{-2}}$} for a flux integrated over a 150 au radius. For the LkCa 15 disk, the results converged toward {$T_{{\rm ex}}\simeq 26.8\pm 2.3$ K, $N_{\rm{tot}}\simeq 1.2\pm 0.1{\,\times 10^{13}}{\,\rm cm^{-2}}$} for a flux integrated over a 200 au radius. The calculated optical depths are $0.07-0.3$ and $0.1-0.4$, toward the MWC 480 and LkCa 15 disks, respectively. While not optically thick the upper boundaries are close enough to $\tau=1$ to be somewhat affected by opacity. Column densities derived using the optically thin approximation may therefore be underestimated by up to a bit more than an order of magnitude. Higher resolution observations and additional CS are needed to better determine the true CS opacity.

We next calculate the column density of CS as a function of the radius, using the same rotational diagram analysis to the azimuthally deprojected radial intensities. The results are presented in Fig. 6, showing, as expected, that the CS column density decreases with increasing radius. We note that in the case of MWC 480, the column density decreases toward the central star, in line with its previously noticed centrally dipped flux profiles.

3.3.2 Isotopologue ratio analysis

The analysis in § 3.3.1 assumed the main isotopologue CS lines are optically thin. This assumption can be explored for one disk, LkCa 15, where two CS isotpopologue lines, ${}^{13}\ce{CS}$ and C34S $6-5$, are well detected.

The ratio of the two isotopologues 12CS and 13CS can be expressed as:

[TABLE]

with $T_{B}$ the brightness temperature (see Appendix B), and provided that the transitions have the same uniform excitation temperature along the line of sight, and emit over the same volume of gas. The former is approximately true when comparing CS and 13CS $6-5$, and the latter is a reasonable approximation if the main isotologue line does not become optically thick high up in the disk atmosphere. Furthermore, if 12CS and 13CS are both optically thin then we have $(1-e^{-\tau})\rightarrow\tau$ which leads to:

[TABLE]

where $R$ is the abundance ratio of 12CS to 13CS. Assuming that their abundance ratio is equal to the 12C/13C ratio determined in the local interstellar medium (LISM), i.e. that their has been no isotopic fractionation, $R\simeq 68\pm 15$. (Milam et al., 2005; Asplund et al., 2009; Manfroid et al., 2009).

For optically thin lines we can calculate the disk-integrated brightness temperatures from the disk integrated fluxes $S_{\nu}\Delta v$ using the Rayleigh-Jeans law:

[TABLE]

To consistently calculate disk integrated fluxes for the main and rare CS isotopologues we applied a uv-taper to the extended baselines of the 12CS visibilities before integrating the fluxes, to match the synthesized beams. As can be seen in Table 2 this reduces the calculated flux by $\sim$25%. For the LkCa 15 disk we then obtain:

[TABLE]

If the 12CS/13CS abundance ratio is equal to the 12C/13C LISM ratio, Eq. 6 implies that the 12CS $6-5$ transition is optically thick. This implication can be used to calculate the 13CS column density and the ${}^{12}\ce{CS}$ excitation temperature. For an optically thick line $(1-e^{-\tau_{{}^{12}\ce{CS}}})\rightarrow 1$. Thus, Eq. B1 directly provides $T_{ex}(\rm{{}^{12}CS})$ from the peak of its brightness temperature:

[TABLE]

Eq. 7 gives an excitation temperature of $T_{ex}(\rm{{}^{12}CS})\simeq 7$ K which, inserted into B1, gives an opacity of $\tau_{\rm{{}^{13}CS}}\simeq 0.06$ (this last step assumes that CS and 13CS have the same excitation temperatures). Moreover, since we assume that 13CS is optically thin we have $(1-e^{-\tau_{{}^{13}\ce{CS}}})\rightarrow\tau_{{}^{13}\ce{CS}}$, which inserted into Eq. 3 gives:

[TABLE]

Eq. 8 gives an opacity consistent with the value obtained from Eq. 7 and B1.

Once the opacity and excitation temperature derived, the column density can be computed as follows:

[TABLE]

This gives a column density of $N_{\rm{tot}}$(13CS)$\simeq 7{\,\times 10^{12}}{\,\rm cm^{-2}}$. Assuming that 12CS is present at the LISM ratio, $N_{\rm{tot}}$(12CS)$\simeq 4.8{\,\times 10^{14}}{\,\rm cm^{-2}}$, which is significantly higher than the rotational diagram result, by a factor of forty (see Table 3). This discrepancy can be understood when considering that $\tau_{\rm{{}^{13}CS}}\simeq 0.06$ thus implies that that $\tau_{\rm{{}^{12}CS}}>1$ (assuming LISM $\ce{{}^{12}C}/\ce{{}^{13}C}$ ratio). We discuss further below whether the LISM assumption and other assumptions implicit or explicit in the analysis are defendable.

A similar computation can be done for the C32S/C34S ratio considering the elemental 32S/34S$=24.4\pm 5.0$ in the vicinity of the Sun (Chin et al., 1996). Both observed 12CS($6-5$)/C34S($6-5$) and 12CS($5-4$)/C34S($6-5$) give a ratio of $\sim$13, which is somewhat lower than the cosmic 32S/34S ratio, but not nearly as much below the cosmic ratio as 12CS/13CS. Taken at face value, the excitation temperature and opacity are $T_{{\rm ex}}$(12CS)$\simeq 7.4$ K and $\tau_{{}^{34}\ce{CS}}=0.07$, leading to $N_{\rm{tot}}$(C34S)$\simeq 6{\,\times 10^{11}}{\,\rm cm^{-2}}$. Using the LISM ratio we thus obtain $N_{\rm{tot}}$(12CS)$\simeq 1.5{\,\times 10^{13}}{\,\rm cm^{-2}}$, which is an order of magnitude lower than the value derived from the 12CS/13CS ratio, and similar to the rotational diagram result.

The results obtained from the different methods presented in §3.3.1 and §3.3.2 to derive CS column densities and excitation temperatures are summarized in Table 3. Both isotopologue analysis suggest that the 12CS line is optically thick, but the large difference in derived 12CS opacities and column densities, as well as the low CS excitation temperatures cast some doubt on the approach. This makes it important to revisit some of the assumptions that went into the calculations. First, these calculations assumed a beam filling factor of unity. Considering the realization that both dust and gas are highly structured in disks (Huang et al., 2017; Favre et al., 2018; Andrews et al., 2018; Huang et al., 2018) this is not obviously true. A beam filling factor of 20 to 40% would increase the excitation temperature to $20-30$ K, in better agreement with our model expectations (see § 4.3). Although Figure 3 shows apparent asymmetries for the 13CS and C34S emissions, these are likely due to noise in the image which propagates to the 0th moment map. So, as such, there is likely little effect on the extracted flux and thus on the isotopologue ratios. Therefore, it seems reasonable to not expect any clear morphological differences between CS and its isotopologues. However, differences in beam-filling factors of the different isotopologues could impact this analysis. Higher signal-to-noise ratio and higher resolution data are needed to improve this analysis. A second assumption is that LTE conditions hold. This should be true if CS is present in the warm molecular layer of disk midplane, but LTE may not hold if CS is instead emitting from the disk atmosphere.

Third, and perhaps most importantly, 13CS and C34S may not be present at the cosmic ratios. There are many examples of isotopic fractionation in disks and the ISM, including in carbon (e.g. Woods & Willacy, 2009; Sakai et al., 2010, 2013; Ossenkopf et al., 2013; Smith et al., 2015; Yoshida et al., 2015; Yu et al., 2016; Huang et al., 2017). For instance, the 13CS enhanced abundances in disk molecular layers could be explained by excessive 13C-carrier, such as 13CO (while 12CO remains self-shielded). Rapid ion-molecule gas-phase reactions could thus transfer 13C from 13CO to hydrocarbons, and then to 13CS. The degree of fractionation that would be needed to explain the 13CS emission is quite high, and detailed modeling on carbon fractionation in CS would be needed and warranted to support this scenario. In the meantime we adopt the results from the rotational diagram analysis, which also agree with the C34S data within error bars, for model comparisons and calculations of column densities in disks without constraints on the excitation temperature.

3.3.3 CS column density estimates for the disk sample

In 4/6 disks we observed a single CS line and therefore have no constraint on the CS excitation temperature. To estimate column densities in these sources, we use the constraints on the CS rotational temperature in the LkCa 15 and MWC 480 disks and apply them to the disk sample. In § 3.3.1 we found disk averaged rotational temperatures of $\sim 23\pm 3$, and $\sim 27\pm 2$ K toward MWC 480 and LkCa 15, respectively. We adopt the average of the two, 25 K, as the typical disk-averaged CS rotational temperature, and an uncertainty of $\pm 5$ K. We then use Eqs. A1 and A2 to calculate the column density. The results are presented in red in Fig. 7. The CS disk-integrated column density varies by less than an order of magnitude across the sample of disks. The average column density is $\approx 7{\,\times 10^{12}}{\,\rm cm^{-2}}$ and the standard deviation $\approx 2{\,\times 10^{11}}{\,\rm cm^{-2}}$. There are no obvious trends with stellar spectral type; the only Herbig Ae star, MWC 480, is close to the sample average. There is also no trend with disk type, i.e. we see no systematic difference in CS column density between full disks and transition disks. It is interesting that the youngest disk in the sample, DO Tau, presents the lowest value, but apart from this there is also no obvious trend with age as V4046 Sgr, the oldest disk in the sample, lies close to the sample average. In the literature some of those disks are found to be a bit colder (e.g. the DM Tau case, Semenov et al., 2018). For colder disk temperatures the derived CS column densities increase, as expected from Eq. A2, leading for instance to an average column density of $\approx 5{\,\times 10^{13}}{\,\rm cm^{-2}}$ with a standard deviation of $\approx\pm 5{\,\times 10^{12}}{\,\rm cm^{-2}}$ for a disk temperature of $10\pm 2$ K, i.e. a factor 5 higher column density.

3.4 \ceH2CS in MWC 480 and LkCa 15

The three \ceH2CS ortho lines ($8_{17}-7_{16}$, $9_{19}-8_{18}$ and $9_{18}-8_{17}$) are detected toward MWC 480 and tentatively detected toward LkCa 15 (2–3$\sigma$). The integrated flux density maps are shown in Fig. 8 and the radial profiles in Fig. 9. The emission toward MWC 480 shows no sub-structure, including no central depression. However, the resolution and signal-to-noise are too low to rule out a similar central dip as is seen in CS. As mentioned in § 3.3.2 for the 13CS and C34S emissions, the quality of the \ceH2CS data is also not high enough to enable an unambiguous distinction between centrally peaked or centrally depleted distributions for these weak lines, and despite appearances in the image plane the emission is consistent with azimuthally symmetric profiles. Higher signal-to-noise ratio and higher resolution data would help in disentangling potential differences in beam-filling factors for the different transitions that could impact the analysis.

\ce

H2CS is an asymmetric rotor with two identical hydrogen nuclei. The molecule is thus found in two different nuclear spin configurations, ortho (where the quantum number $K_{a}$ is odd) and para (where the quantum number $K_{a}$ is even). The three \ceH2CS lines presented here are ortho lines, which can be used to derive a column density of the $o$-H2CS, and a total column density assuming an $o/p$ ratio. This is done for MWC 480 in Fig. 10 using the rotational diagram method described Sect. 3.3.1. The following flat priors were assumed:

[TABLE]

Random draws and results from the posterior distributions are depicted in gray in Fig. 10, along with the least square fit results, shown in blue. The results converged toward {$T_{{\rm ex}}\simeq 41^{+38}_{-20}$ K, $N_{\rm{tot}}\simeq 3.0^{+2.7}_{-0.4}{\,\times 10^{12}}{\,\rm cm^{-2}}$}, assuming the statistical ortho-to-para ratio of 3, toward the MWC 480 disk. The large uncertainties are due to a combination of low signal to noise and closely spaced energy levels. Still, it is noteworthy that the best-fit rotational temperature is close to the rotational temperature of CS, assuming CS is optically thin.

We used Eqs. A1 and A2 to derived an upper limit of $N_{\rm{tot}}\lesssim 2{\,\times 10^{12}}{\,\rm cm^{-2}}$ when considering the three \ceH2CS line tentatively detected toward LkCa 15, a rotational temperature of $T\approx 40$ K and an ortho-to-para ratio of 3. This corresponds to less than a factor of 2 lower than the value observed toward MWC 480.

4 Astrochemical modeling

To explore the origin of the observed sulfur chemistry properties across our disk sample in general, and in MWC 480 and LkCa 15 in particular, we use a 1+1D protoplanetary disk model, based on the gas-grain chemical model Nautilus (v.1.1) (Wakelam et al., 2016), tuned to the physical structures of MWC 480 and LkCa 15.

4.1 Physical structure

Assuming azimuthally symmetric disks, PPD physical parameters are commonly described in cylindrical coordinates centered on the inner star with a radius $r$ and the vertical axis $z$ perpendicular to the disks. Thus, in the following, we describe our PPD modeling along these two radial and vertical axes. The disk models extend from an inner radius of 1 au to an outer radius of 400 au. The other physical parameters used to compute the physical structures of MWC 480 and LkCa 15 are listed in Table 4. We use the previously calculated temperature and density profiles from Piétu et al. (2007) and Guilloteau et al. (2011),assuming that the temperature and surface density, governing line emissions, vary as power laws of the radii (Beckwith et al., 1990; Piétu et al., 2007). Below we describe the parameterization of temperature and density, visual extinction and UV fields, we used for our 1+1D disk modeling.

4.1.1 Temperature profile

At a given radius $r$ from the central star, the vertical temperature gradient is computed following the modified prescription of Dartois et al. (2003) used in Rosenfeld et al. (2013) and Williams & Best (2014):

[TABLE]

where

[TABLE]

vary as power law of the radii (Beckwith et al., 1990; Piétu et al., 2007). $R_{c}$ is a characteristic radius and $z_{q}=4H$, where $H$ is the pressure scale height modeled as the height at which the optical depth to UV photons, $A_{{\rm V}}$, becomes $\lesssim 1$. $H$, defined by the midplane temperature and assuming vertical static equilibrium, can be expressed as follows:

[TABLE]

with $k_{\rm{B}}$ the Boltzmann constant, $\mu=2.4$ the reduced mass of the gas, $m_{\rm{H}}$ the proton mass, $G$ the gravitational constant, and $M_{\star}$ the mass of the central star. The midplane temperature $T_{\mathrm{mid},R_{c}}$ can be estimated from the following approximation for a simple irradiated passive flared disk:

[TABLE]

with $L_{\star}$ the stellar luminosity, $\sigma_{\rm{SB}}$ the Stefan-Boltzman constant and $\varphi$ the flaring angle (e.g., (e.g. Chiang & Goldreich, 1997; D’Alessio et al., 1998; Dullemond et al., 2001; Dullemond & Penzlin, 2018; Huang et al., 2018).

The dust temperature is assumed equal to that of the gas. This is not a valid assumption in the uppermost layers of the disk (Bergin et al., 2007), but is a reasonable simplification in this study focusing on S-chemistry in the molecular disk layers (e.g. Woitke et al., 2009; Akimkin et al., 2013; Dutrey et al., 2017; Semenov et al., 2018). For instance, in the Flying Saucer edge-on T Tauri disk, the CS emission is observed to be less vertically extended than the CO one (Dutrey et al., 2017). The resulting 2D temperature profiles for the MWC 480 and LkCa 15 disks are shown in the first column panels of Fig. 11.

4.1.2 Density profile

The disk is assumed to be in hydrostatic equilibrium. Thus, for a given vertical temperature profile (see sect. 4.1.1), the vertical density structure is determined by solving the equation of hydrostatic equilibrium:

[TABLE]

which gives for the isothermal case:

[TABLE]

with solution:

[TABLE]

where $\rho_{o}$ corresponds to the volume density of the gas in the midplane. The surface density of the disk, $\Sigma$, can thus be expressed as follows:

[TABLE]

Furthermore, 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) with sharp cutoff at specified inner and outer radii (i.e. $R_{\rm{int}}<r<R_{\rm{out}}$):

[TABLE]

where $\Sigma_{c}$ is the surface density at the characteristic radius.

During its evolution, the disk loses materials toward the central star which can be quantified by the accretion rate $\dot{M}=-2\pi r\Sigma(r)v_{r}$, with $v_{r}$ the radial velocity of the disk. For axis-symmetric disk, mass conservation thus reads:

[TABLE]

If the mass loss of the disk is constant with time, i.e. $\frac{\delta\Sigma}{\delta t}=0$, the disk mass can simply be expressed as follows:

[TABLE]

Thus, the surface density at a radius $R$ reads:

[TABLE]

which highlights the link between the mass of the disk, $M_{\rm{disk}}$, and its outer radius, $R_{\rm{out}}$.

The resulting 2D density profile is represented in the second column panels of Fig. 11, for each model.

4.1.3 Visual extinction profile

The vertical visual extinction gradient is computed from the hydrostatic density profile using the conversion factor $N_{\ce{H}}/A_{\rm{V}}=1.6{\,\times 10^{21}}$ (Wagenblast & Hartquist, 1989), with $N_{\ce{H}}=N(\ce{H})+2N(\ce{H2})$ the vertical hydrogen column density of hydrogen nuclei. This conversion factor assumes a typical mean grain radius size of 0.1 $\mu$m and dust-to-mass ratio of 0.01, consistent with model assumptions. The resulting 2D visual extinction profile is represented in the third column panels of Fig. 11, for each model.

4.1.4 UV flux distribution

In the disk model the simplifying assumption of parameterizing the UV flux from the star as a multiple of the standard ISM radiation field in Draine (1978)’s units is done. This is a reasonable assumption for UV excess spectra of Herbig Ae stars, which are found to have UV field shapes similar to the interstellar radiation field (ISRF) (Chapillon et al., 2008). For T Tauri stars UV excess spectrum can be dominated by accretion luminosity and therefore Lyman alpha, instead. Thus, for T Tauri stars, we are probably underestimating UV penetration into the disk, since Ly-alpha photons can scatter more efficiently into the disk (Bergin et al., 2003). However, for simplicity and ease of comparison we use this approximation in both disk models presented here.

The unattenuated UV flux factor, $f_{\rm{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. Since the physical structure is built from a parametric model approach (i.e. without radiative transfer treatment), the input UV flux of reference $f_{\rm{UV},R_{c}}$ (i.e. the UV flux factor at a radius $R_{c}$, see Table 1), is divided by two (e.g. Wakelam et al., 2016). This 0th order approximation assumes that half of photons irradiating from the embedded young star diffuse into the disk, and half in the perpendicular direction. The UV flux factor at a given radius $r$ thus reads:

[TABLE]

To calculate the UV flux at any point in the disk the UV flux factor is convolved with the visual extinction profile. The resulting 2D UV flux is represented in the last column panels of Fig. 11, for each model.

4.2 Chemical modeling

Once built, the 1+1D physical structure described above (and depicted in Fig. 11 for the MWC 480 and LkCa 15 cases) is fed as an input in the pseudo time-dependent astrochemical modeling code Nautilus v1.1 (Wakelam et al., 2016) used in three-phase mode (Ruaud et al., 2016). This rate-equation gas-grain code – based on Hasegawa et al. (1992) and Hasegawa & Herbst (1993) – simulates the time-dependent chemistry of $\sim 1100$ species (half in the gas phase and half in solid phase) linked together via more than $\sim 12000$ reactions, in the vertical direction at each radius in three different phases: gas, grain surface (top two ice layers on grains), and grain mantle (deeper ice layers on grains). Exchanges in between all the different phases are included: adsorption and desorption processes link the gas and surface phases, and swapping processes link the mantle and surface of grains. Several desorption mechanisms are considered, thermal desorption (Hasegawa et al., 1992), and non thermal ones, such as cosmic-ray induced desorption (Hasegawa & Herbst, 1993), photodesorption and chemical desorption (for further details see e.g. Garrod et al., 2007; Ruaud et al., 2016; Wakelam et al., 2016; Le Gal et al., 2017). For photodesorption we use the prescriptions discussed in Le Gal et al. (2017), where we followed the recommendations of Bertin et al. (2013) to use a simplistic approach. This consists in considering a single photo-desorption yield for all the molecules rather than individual ones, only experimentally determined for a few number of species. Thus, we set a constant generic photo-desorption yield of $1{\,\times 10^{-4}}$ molecule/photon (Andersson & van Dishoeck, 2008) for all the species contained in our chemical network. 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. A dust-to-mass ratio of 0.01 is assumed, with spherical grains of $0.1~{}\mu$m single radius size.

A key parameter of the new modeling study we present here, is the use of an up-to-date sulfur chemical network, based on the KInetic Database for Astrochemistry (KIDA) (http://kida.obs.u-bordeaux1.fr/), and including recent updates (Vidal et al., 2017; Le Gal et al., 2017; Fuente et al., 2017), that are relevant for the purpose of this work.

To take into account chemical inheritance from previous stages, we first simulate the chemical evolution of a starless dense molecular cloud up to a characteristic age of $1{\,\times 10^{6}}$ years (e.g. Elmegreen, 2000; Hartmann et al., 2001). For this 0D model typical constant physical conditions were used: grain and gas temperatures of 10 K, a gas density $n_{\rm{H}}=n(\ce{H})+n(\ce{H2})=2{\,\times 10^{4}}{\,\rm cm^{-3}}$ and a cosmic ionization rate of $1.3{\,\times 10^{-17}}{\,\rm s^{-1}}$; this parent molecular cloud is also considered to be shielded from external UV photons by a visual extinction of 30 mag. For this first simulation stage, we consider diffuse gas starting conditions, i.e. that initially all the elements are in atomic form (see Table LABEL:tab:elemental_ab) except for hydrogen assumed to be initially 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 LABEL:tab:elemental_ab. The chemical gas and ice compositions of this representative parent molecular cloud serve as the initial chemistry for our 1+1D disk model.

We ran the disk chemistry models for one million years. This is on the low side for our disk sample (see Table 1), but has been shown in previous studies to be a reasonable approximation of the chemical age of a disk when grain growth is not taken into account (e.g. Cleeves et al., 2015). Even though the chemistry in the disk has not reached steady state at such a time, its evolution is slow enough that the results presented here hold for a disk twice younger or older. We note that the purpose of these models is not to quantitatively fit observations, but rather to obtain a first intuition on whether existing sulfur chemistry networks can reproduce observed sulfur species.

With the disk model presented above, we thus ran two different disk chemistry simulations that differentiate from one another only by their physical structures: one simulating the MWC 480 disk and a second the LkCa 15 disk (see Fig. 11). Both models start from the same initial chemical composition, that results from the 0D model described above and considers a ”low” sulfur elemental abundance of $8{\,\times 10^{-8}}$ (see Table LABEL:tab:elemental_ab, Graedel et al., 1982; Wakelam & Herbst, 2008), i.e. $\sim 2$ orders of magnitude below the cosmic value $\sim 1.5{\,\times 10^{-5}}$ (Asplund et al., 2009). This ”low” sulfur elemental abundance corresponds to the typical sulfur depleted value used in models to reproduce the abundances observed in dense astrophysical molecular environments (e.g. Tieftrunk et al., 1994; Wakelam et al., 2004; Vidal et al., 2017; Vastel et al., 2018). For the present study, this approximation seems reasonable since CS and \ceH2CS are supposed to reside in disk molecular layers, which are similar to dense astrophysical molecular environments.

4.3 Model results

The resulting number densities of CS, \ceH2CS and their main precursors, \ceS+ and S, at 1 Myr are presented in Fig. 12 for the MWC 480 and LkCa 15 disk models (the corresponding fractional abundances are shown in Appendix D). Fig. 12 shows that in the warmer and denser Herbig Ae disk MWC 480, all gas-phase S carriers occur in higher disk layers than in the colder T Tauri disk LkCa 15 (see Table 4). Even if a common expectation is that the higher UV fluxes found around more massive stars will push the molecular layers down to deeper disk layers, since the disk is more massive, the visual extinction is much larger (as can be seen in Fig. 11), and so the overall outcome is that UV does not penetrate. MWC 480 has a factor of 6 larger disk mass than LkCa 15. In the more tenuous LkCa 15 disk the shielded layers thus occur deeper into the disk than would be expected if it was as massive as MWC 480. Furthermore, chemistry is not just set by UV radiation, but also by density, which, e.g., controls the recombination rate. When comparing Fig. 12 with Fig. 11, we see that the layers where CS and H2CS reside have similar densities of $\sim{\,\rm 10^{6}}-{\,\rm 10^{8}}{\,\rm cm^{-3}}$.

In both disks the S and \ceH2CS layers reside deeper toward the disk midplane than \ceS+ and \ceCS. The depth of the S/\ceS+ transition is set by a balance of, on the one hand, production of S+ from S through photo-ionization and charge transfer reaction with \ceC+ and, on the other hand, recombination of S+ with electrons to form S. The main formation pathway to form CS starts with S+, while H2CS mainly forms from atomic S, explaining their difference in vertical distribution.

In more detail, the CS-producing sulfur chemistry in the upper disk layers is driven by rapid ion-neutral reactions between \ceS+ and small hydrocarbons \ceCHx and \ceCyH (with $x=1-4$ and $y=2-3$). This produces carbonated S-ions including \ceHCS+, \ceCS+, \ceHC3S+, and \ceC2S+, that subsequently recombine with electrons to form neutral S-containing species. CS is thus a daughter molecule of \ceS+. CS also has slower neutral-neutral formation pathways starting with S and neutral small hydrocarbons, which produces a second CS reservoir at deeper disk layers. Its main destruction pathways are photodissociation processes mainly in the disk atmosphere, reactions with protons and protonated ions (i.e. with \ceH+, \ceH3+, \ceHCO+) and with \ceHe+, and freeze-out onto grain surfaces.

\ce

H2CS is mainly produced by the gas-phase neutral-neutral reaction of atomic S and \ceCH3 and the electronic dissociative recombination of \ceH3CS+, whose formation begins with the reaction between \ceS+ and \ceCH4. At high densities and low temperatures, i.e. in regions closer to the midplane, some grain-surface reactions can also contribute to its formation, such as the hydrogenation reaction \ces-H + s-HCS -¿ s-H2CS releasing some ($\sim 1\%$) of its \ceH2CS product in the gas-phase by chemical reactive desorption. Similarly to CS, \ceH2CS is mainly destroyed by photodissociation processes in the disk atmosphere, reactions with protons and protonated ions (i.e. with \ceH+, \ceH3+, \ceHCO+), and freeze-out onto grain surfaces.

The resulting radial distributions of the vertically integrated column densities along with the observationally constrained column densities are presented in Fig. 13. For both the MWC 480 and LkCa 15 models, the CS column density peaks at a radius of $\sim 100$ au, and decreases from there with increasing radius. The theoretical outer disk CS column density profiles appear similar for both disks, but in the inner disk, the LkCa 15 column density dips, while the MWC 480 column density peaks. Apart from this inner peak, the theoretical LkCa 15 column densities are systematically higher than those of MWC 480, by a factor of $\sim 2-10$, consistent with observations. Across the disks, the model predictions are generally within a factor of a few for both MWC 480 and LkCa 15, which must be considered as good agreement, considering that there is no tuning of the model involved. We also note, that the CS column density derived from the 13CS observations toward LkCa 15 being an order of magnitude higher are thus inconsistent with our model predictions. However, we note that our disk models were run with an initially depleted S-element abundance (see § 4.2). Running the same models with undepleted S-element abundance results in an overprediction of CS by 1-2 orders of magnitude, while \ceH2CS is in better agreement compared to the depleted models. This emphasizes the need for additional S-species disk observations to anchor S-disk chemistry modeling.

The predicted H2CS column densities are 1–2 orders of magnitude below the predicted CS column densities in each disks. Similarly to CS, the H2CS profiles are characterized by fall-off in the outer disk, a local maximum at $50-100$ au, and then (similar to CS in MWC 480) an increasing column density toward the central source. The dips in the inner disks are not seen in the CS and H2CS precursors, with the possible exception of S toward LkCa 15. Compared with observations toward MWC 480, the model underpredicts H2CS by several orders of magnitude. Furthermore the model predicts that more H2CS should be present toward LkCa 15 compared to MWC 480, which is opposite to what we observed. Finally the model predicts that H2CS should be centrally very peaked, while observations suggest more extended emission.

Figure 14 presents the calculated mean disk gas temperature weighted by the density of CS and H2CS, respectively in the two disks, which provides a measure of the typical temperature environment within which these molecules reside. In the outer disk (i.e. for radius $>200$ au), CS resides at temperatures in the range $\sim 8-12$ K and $\sim 20-30$ K, in the LkCa 15 and MWC 480 disk models, respectively. In the inner disk (i.e. for radius $\lesssim 50$ au) they reside at temperatures in the ranges $\sim 25-200$ K and $\sim 60-400$ K, respectively. Similar behaviors can be seen for \ceH2CS. While the H2CS density peaks in more embedded disk layers, the vertical temperature gradient is small in our model, and the difference in temperature when compared with CS is therefore also small (see Fig. 12). The rotational temperatures of CS and \ceH2CS are also overplotted on Fig. 14. These temperatures were derived from observations out to $\sim$190 au (MWC 480) amd $\sim$240 au (LkCa 15), and were found to be quite similar. Compared to CS observations, the models overpredict the temperature for MWC 480, and underpredict it for LkCa 15.

5 Discussion

5.1 Model predictions vs observations

In § 4.3 we compared model predictions of CS and H2CS column densities, radial profiles and excitation temperatures. In this section we summarize agreements and disagreements and discuss possible explanations for cases where there is a mismatch between observations and theory.

Column densities: The observed CS column densities of 1012–10${}^{13}{\,\rm cm^{-2}}$ are well reproduced, within a factor of a few, by our models. The best fit are found for radii in the range $\sim 150-250$ au for MWC 480, and for radii of $\sim 100$ au and in the range $\sim 200-275$ au for LkCa 15. However, we should note the discrepancies obtained for the inner disk of MWC 480 (i.e. radii $\lesssim 100$ au) between the CS column density profiles derived from the observations and the one obtained from our model. Both higher angular resolution observations and further modeling investigations are required to disentangle this issue. By comparison, the models underpredict H2CS in the outer disk of MWC 480, and incorrectly predict that H2CS should be more abundant in LkCa 15 than in MWC 480. An explanation could be related to the fact that, for simplicity, we do not include LkCa15’s wide inner cavity (e.g. Alencar et al., 2018, and reference therein) in the LkCa disk structure. This will likely impact some aspects of the chemistry model predictions, especially at radii of $<50$ au. In the outer disk, the model predictions should be more reliable, or at least not limited by this particular omission. Future disk modeling that takes into account this structural feature and its impact on radiative transport is needed to make quantitative model-data comparisons in the inner disk. A possible reason for the overall underprediction of H2CS is that grain surface formation of H2CS and subsequent photodesorption may be more important in disks than is currently implied by the model. The surface and bulk ice H2CS 2D disk densities are represented in Fig. 15, as well as those of CS, for both MWC 480 and LkCa 15 disk model. We see that both molecules are abundant in ices, especially in the disk midplanes. However, photo-irradiation laboratory experiments of \ceCH3OH ices showed an efficient dissociation of the product rather than a desorption (Bertin et al., 2016; Cruz-Diaz et al., 2016). \ceH2CS being also a rather large molecule, it could similarly efficiently be photo-dissociated in the ice. Our present model already considers ice photo-dissociation rates, however, due to the lack of laboratory studies on that topic, those are set to be equal to the gas-phase photo-dissociation ones. So, our model might misspredict the \ceH2CS ice abundance and desorption, unless additional ice formation pathways are missing in our astrochemical networks. Therefore, both, photo-irradiation of icy \ceH2CS and new \ceH2CS ice formation pathways require laboratory experiments. Moreover, the presented model does not take into account grain growth and migration, which may have depleted the outer disk of UV protection, increasing both its temperature (Cleeves, 2016) and the importance of ice photodesorption in this region (Öberg et al., 2015a). This needs to be revisited in a future model that incorporates this structural element, where the ISRF impact may increase if the outer disk is more tenuous. Interestingly, the MWC 480 continuum disk appears more compact than the LkCa 15 disk, and ice photodesorption may therefore be important in a larger portion of the MWC 480 disk, possibly explaining why it hosts more H2CS in the gas-phase compared to the LkCa 15 disk. In addition, the sulfur chemistry network is probably still incomplete, and there may be other gas- and grain-surface pathways to \ceH2CS that are not yet considered in the model. Finally, we should also note that beam dilution could be an additional explanation for the mismatch observed in between the \ceH2CS observations and its predictions from the astrochemical disk modeling. Thus, higher \ceH2CS angular resolution observations would be required toward the same disk sample to further investigate these issues.

Radial profiles: Our beam size corresponds to 80 au, and we are therefore not sensitive to radial scales below 40 au. On larger scales there is qualitative agreement between observed and modeled CS radial profiles, including the central dip seen toward MWC 480. However, the modeled predicted dip appears deeper than the one derived from the observations when using the rotational diagram method (Fig. 6), but this might in part be a question of radiative transfer. However, in all other disks we observed a central CS peak. It would thus be interesting to further investigate if this is just a spatial resolution issue or if this characterizes a specific difference in between the Herbig Ae MWC 480 disk and the other T Tauri disks of our sample. To address this issue additional observations would be required, both in the same disk sample with higher angular resolution and toward additional Herbig Ae disks to increase our statistics on the latter type of disk. Another interesting feature to further explore is the break observed in CS at the continuum edges across our whole disk sample. To address this issue, the coupling between dust and gas need to be refined in our current modeling.

Excitation temperatures: Model predictions for CS excitation temperatures are too high for MWC 480 and too low for LkCa 15. This is probably a result of poorly constrained vertical temperature structures in disk models rather than an issue with chemistry model predictions under which conditions CS is present in the gas-phase; the CS chemistry is relatively straight forward and should mainly depend on the presence of ionized S, and gas-phase carbon. We suspect that the LkCa 15 disk is modeled to be too cold. The model temperature structures of both disks were calculated using some rather simplistic model assumptions and need to be updated to fit more recent observations than the ones available. This will be the subject of a future theory-focused study. According to the model predictions, the S-bearing species presented in Fig. 12 reside at intermediate scale heights for the LkCa 15 disk and at slightly higher ones for the warmer disk MWC 480. This could result from the too simplistic physical structure modeling performed in the present study and would benefit from physical structure modeling refinement and improvement. However, as we mentioned in § 4.3, the MWC 480 disk being more massive, its visual extinction is larger (see Fig. 11), implying less UV penetration in the disk. In addition, the dramatic changes in temperature over small radii shown in Fig.14 highlight, once again, that our observations suffer from insufficient spatial resolution to effectively trace this trend, and therefore higher resolution CS observations are required to test this model result.

High spatial resolution observations of additional sulfur-bearing molecules in the same disks are needed to better constrain disk sulfur chemistry. Furthermore, observations toward other Herbig Ae disks would also help in improving our sulfur chemistry knowledge since, according to the comparison between the MWC 480 observations presented here versus the T Tauri ones, there might be some sulfur chemistry distinctions in between T Tauri and Herbig Ae disks.

Also, we should note that the impact of X-rays on disk chemistry is not included in our present model. X-rays could affect the the ion abundances (e.g. Glassgold et al., 1997; Aikawa & Herbst, 1999; Henning et al., 2010; Teague et al., 2015; Cleeves et al., 2015) such as the abundances of \ceS+, \ceCS+, \ceHCS+ as well as the SO/\ceSO2 formation, and thus affects the overall S-chemistry (e.g. Semenov et al., 2018). However, for the ALMA observations presented here, concerning disk radii beyond $\sim 10-20$  au, X-rays chemistry is not dominant (e.g. Rab et al., 2018). Nevertheless, this would be an interesting future modeling improvement.

5.2 S-carriers in disks and other astrophysical environments

Compared to other disk studies, the CS column densities derived toward our disk sample (see Fig 7) are consistent with the values found in other disks, e.g. $2.0\pm 0.16{\,\times 10^{12}}{\,\rm cm^{-2}}$ in GO Tau and $8.7\pm 1.6{\,\times 10^{12}}{\,\rm cm^{-2}}$ in LkCa 15 (Dutrey et al., 2011), $6\pm 3{\,\times 10^{12}}{\,\rm cm^{-2}}$ in DM Tau (Semenov et al., 2018) and $2.2{\,\times 10^{13}}{\,\rm cm^{-2}}$ in GG Tau (Phuong et al., 2018), which agrees with CS being the most abundant S-species detected in disks so far.

In order to place our CS and \ceH2CS disk observations in a more hollistic view of sulfur chemistry, we explore in this section what characterizes the sulfur chemistry in different interstellar environments, based on observations from the literature reported in PDR (Riviere-Marichalar et al. submitted), molecular clouds (Drozdovskaya et al., 2018) and cold dense core (Vastel et al., 2018), analogs of the main three chemical layers constituting protoplanetary disks, i.e. the disk atmosphere, intermediate molecular layers, and the disk midplane. Figure 16 presents the column density ratios of the main S-bearing species observed thus far in PPD (from this work and previous studies) versus CS, compared to their values in the different ISM environments aforementioned.

The 13CS/12CS value we derived in disks from this study is found higher ($\sim 0.06$) than those measured so far in the different ISM environments presented in Fig. 16 ($\sim 0.01-0.03$). However, the wide error bars we derived ($\pm 0.05$) do not allow us to draw strong conclusions on the plausible 13C enhancement suggested in this work. Higher resolution observations of CS isotopolgues in disks are needed to further investigate this 13C enhancement hint.

The C34S/C32S ratio we observed in disks ($\sim 0.08$) is found, within error bars, a factor of 2 higher than values found in PDRs and dense cores, and than the LISM value. The observed C34S/C32S ratio found in this work is also about an order of magnitude higher compared to the value found in the molecular envelope of the IRAS 16293-2422 Sunlike protostar. Further observational explorations of this ratio from protostar envelopes to disks could thus bring some clues on the S-chemistry, both, on the question of chemical inheritance versus chemical reset during the stellar formation and on the physical environment impact (UV flux, temperature, etc.) on it.

The \ceH2CS/CS ratio decreasing from dense to diffuse ISM, is found in disks (i.e. in LkCa 15) almost as high as in dense cores (Vastel et al., 2018). Nevertheless, the high uncertainties on our measurement do not allow us to firmly assert that the \ceH2CS/CS ratio behaves in disk more like dense core S-chemistry than diffuse ISM S-chemistry. To arbitrate on that, higher resolution observations are also needed.

Interestingly, and in opposition to the \ceH2CS/CS ratio, the \ceH2S/CS is increasing from dense to diffuse ISM and the disk value (Phuong et al., 2018) is found almost as low as the lower limit derived in dense core (Vastel et al., 2018). This could also explain why \ceH2S was so hard to detect in disks (e.g. Dutrey et al., 2011), aside from the fact that the \ceH2S detection could have been facilitate by the high mass of the GG Tau disk in which it was detected. Also, even though high uncertainties remain on value observed toward IRAS 16293-2422, \ceH2S appears as one of the main S-reservoirs in this source. A similar result is found for the OCS/CS ratio toward IRAS 16293-2422. These last two results could reflect more efficient thermal desorption in this source, that would also explain the boost in S-chemistry observed there (Drozdovskaya et al., 2018). However, and in opposition to the \ceH2S/CS ratio, the OCS/CS ratio is found lower in PDR than in dense core, suggesting that the OCS and \ceH2S molecules trace different type of environments.

As the \ceH2S/CS ratio, the SO/CS ratio increases from dense to diffuse ISM. But, only a lower limit on SO is provided toward IRAS 16293-2422 (Drozdovskaya et al., 2018). Semenov et al. (2018) provided upper limit on the SO column density in disks (i.e. toward DM Tau), giving an SO/CS ratio lower than the lower limit derived toward dense core (Vastel et al., 2018). Contrary to the SO/CS, the \ceSO2/CS ratio appears rather flat across the different environments compared here, and so do not seem to characterize any specific S-chemistry behavior that would depend on its physical environment. PDR value is found compatible with cold core trend, and with the IRAS 16293-2422 value but to a lesser extent regarding the large uncertainties found. Similarly, the upper limit set by Semenov et al. (2018) do not allow us to settle any specific trend for this ratio.

Finally, the \ceC2S/CS ratio decreases from dense to diffuse ISM, in an even more pronounced way than the OCS/CS and \ceH2S/CS ratios do. In this case, the disk upper limit (Chapillon et al., 2012) found in between the value from the dense prestellar core L1544 (Vastel et al., 2018) and the Horsehead nebula cold core position (Riviere-Marichalar et al. submitted), thus seems to be characteristic of molecular cloud environments.

On average, the S-bearing molecules observed in PPD would behave more like in cold core or molecular clouds. In particular, Fig 16 highlights that CS is one of the most abundant S-species detected so far in the different astrophysical environments discussed in this section. However, additional S-bearing PPD observations are required to test these preliminary conclusions.

6 Conclusion

We presented here new ALMA observations of sulfur-bearing molecules toward a sample of six protoplanetary disks. Our main findings are summarized below:

The CS $5-4$ rotational transition was observed and detected toward five Taurus PPDs (DM Tau, DO Tau, CI Tau, LkCa 15, MWC 480), and the CS $6-5$ transition toward three PPDs (LkCa 15, MWC 480 and V4046 Sgr). 2. 2.

We used the two main CS isotopologue rotational transitions ($6-5$ and $5-4$), toward the MWC 480 and LkCa 15 disks, to derive their respective column densities and rotational temperatures of $N_{\rm{tot}}$(CS)$\simeq(7.3\pm 0.8){\,\times 10^{12}}{\,\rm cm^{-2}}$, $T_{{\rm ex}}\simeq 23\pm 2.5$ K and $N_{\rm{tot}}$(CS)=$(1.2\pm 0.1){\,\times 10^{13}}{\,\rm cm^{-2}}$, $T_{{\rm ex}}=27\pm 2.3$ K, assuming LTE and optically thin emission. 3. 3.

We used the average CS rotational temperature derived toward MWC 480 and LkCa 15, i.e. $T_{{\rm ex}}\simeq 25\pm 5$ K, to estimate the CS column densities toward our full disk sample, resulting in an average CS column density of $N_{\rm{tot}}$(CS)$\approx 7{\,\times 10^{12}}$, with a standard deviation of $\approx 2{\,\times 10^{11}}$. 4. 4.

We report the first detections in disks of 13CS, C34S (in LkCa 15), and \ceH2CS (in MWC 480). 5. 5.

We used the CS isotopologue detections toward LkCa 15 to assess the CS disk opacity and tentatively conclude that either both 13C and 34S are enhanced by factors of a few (the former showing stronger enhancement than the latter) or the main CS isotopologue lines are optically thick, or there is substantial sub-structure that causes beam dilution. 6. 6.

From the three \ceH2CS line detections obtained toward MWC 480 we derived an \ceH2CS column density and rotational temperature of $N_{\rm{tot}}$(\ceH2CS)$\simeq 3{\,\times 10^{12}}$, $T_{{\rm ex}}\simeq 41$ K. Toward LkCa 15 we use the tentative line detections to derive a less than a factor of 2 lower column density upper limit. 7. 7.

Comparing our CS disk observations to 2D disk astrochemical modeling shows that the CS chemistry seems to be relatively well understood, with column density predictions within a factor of a few compared to the observations found in both MWC 480 and LkCa 15. 8. 8.

We also compared our \ceH2CS observations to 2D disk astrochemical modeling. For that case, our model underpredict \ceH2CS in MWC 480 by 1-2 orders of magnitude and also incorrectly predicts that it should be more abundant in LkCa 15 than in MWC 480, emphasizing the need for further S-disk theoretical investigations. 9. 9.

Concerning S-carriers in disks, our study agrees with previous work finding that CS constitutes one of the main abundant S-species in disk and that it can be produced, at least partly, in protoplanetary disk stage itself. 10. 10.

Furthermore, the main disk strata - the disk atmosphere, intermediate molecular layers and the disk midplane - are analogs to ISM astrophysical environments including PDRs, protostar molecular envelops and dense core, respectively. When comparing the S-bearing column density ratios observed so far in disks (i.e. from this study and literature ones) with their corresponding values observed in the ISM astrophysical environments aforementioned, we find that S-disk chemistry appears more closely related to the one observed in molecular cloud environments. This is consistent with model results which places the origins of the observed emission in the disk molecular layer.

Going forward, additional S-bearing PPD observations are required to increase the variety of S-species observed in disks and construct robust comparisons between S-disk chemistry and earlier ISM evolutionary stages. Such comparisons are needed both to further characterize what kind of chemistry drives the sulfur chemical evolution in disks, and to assess the degree of inheritance from the natal molecular cloud.

Higher resolution and better signal-to-noise ratio observations would also help to solve any beam dilution issue and to disentangle whether or not substantial carbon and/or sulfur isotopic enhancements occur at protoplanetary disk stage. Such studies combined with isotopic measurements of sulfur carriers in comets would further aid in addressing the inheritance question.

Additional excitation lines, observed with higher signal-to-noise ratio would allow for non-LTE modelling and therefore better constraints on both the CS and \ceH2CS column densities and line emission regions. Non-LTE modeling would also be helpful to better constrain the CS and \ceH2CS observations. For instance, it could help in better identifying in which disk layers these respective S-species dominate. On the astrochemical modeling side, there is a need both for more refined disk structures, and for more complete sulfur chemical networks.

Finally, we are ultimately aiming to address how much of the volatile sulfur will be incorporated into the nascent planets? In addition to better constraints on disk sulfur chemistry, addressing this question requires models that couple chemistry and dynamics from disk formation to planet assembly. This long-term goal is needed to address the sulfur budgets on planets and thus the planet hospitality to the origins of life (Ranjan et al., 2018).

Acknowledgments

This paper uses the following ALMA data: ADS$/$JAO.ALMA#2016.1.01046.S,

ADS$/$JAO.ALMA#2016.1.00627.S, and

ADS$/$JAO.ALMA#2013.1.01070.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), NINS (Japan), together with NRC (Canada), NSC, ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory (JAO) is operated by ESO, AUI/NRAO and NAOJ. The authors would like to thank the two anonymous referees for valuable suggestions and comments. RLG thanks Valentine Wakelam for allowing her to make use of the Nautilus code in independent research. This work was supported by an award from the Simons Foundation (SCOL # 321183, KO). J.B.B. acknowledges funding from the National Science Foundation Graduate Research Fellowship under Grant DGE1144152.

Appendix A Building rotational diagrams from disk-integrated flux densities

Assuming optically thin transitions, the disk-integrated flux densities $S_{\nu}\Delta v$, can be related to the column density of their respective upper energy state, $N_{u}$, as follows:

[TABLE]

where $S_{\nu}$ is the flux density, $\Delta v$ the line width, $A_{ul}$ the Einstein coefficient, $c$ the speed of light and $\Omega$ the solid angle subtended by the source (e.g. Bisschop et al., 2008; Loomis et al., 2018). For this analysis we use the disk flux densities $S_{\nu}\Delta v$ integrated out to 150 au and 200 au in the MWC 480 and LkCa 15 disks respectively, referred to Rσ in Table 2.

The total column density, $N_{\rm{tot}}$, and rotational temperature, $T$, can then be derived from the upper level population, $N_{u}$, following the Boltzmann distribution:

[TABLE]

with $g_{u}$ and $E_{u}$ the degeneracy and energy of the upper energy level $u$, respectively; $k_{\rm{B}}$ the Boltzmann constant; and $Q_{\rm{rot}}$ the partition function of the molecule, which for a diatomic molecule as CS can be approximated by:

[TABLE]

(Gordy & Cook, 1984). $h$ is the Planck constant and $B_{0}$ the rotational constant for CS.

In practice, this can be done by applying a linear least square regression to the logarithm of Eq. A2, as depicted in blue in Fig. 5, but this fit does not take into account that the lines may be slightly optically thick. This can be included in the fit by multiplying $N_{u}$ by the ”optical depth correction factor” for a square line profile, $C_{\tau}=\frac{\tau}{1-e^{-\tau}}$, in case $\tau~{}\cancel{\ll}~{}1$ (Goldsmith & Langer, 1999). The resulting logarithm of Eq. A2 is then:

[TABLE]

In the LTE approximation, the optical depth of a given transition at $T_{{\rm ex}}$ can be expressed as a function of its upper energy level population using Eq. A2 as follows:

[TABLE]

As mentioned in § 3.3.1, the observed CS lines are optically thin based on their brightness temperatures, i.e $\tau<1$, in which case the line profiles remain Gaussian, and the integral of the opacity is given by:

[TABLE]

where $\Delta v_{\rm{FWHM}}=\sqrt{8\ln 2}\,\sigma_{v}$ is the full width at half maximum of the observed transition. Thus, Eq. A5 can be rewritten as follows:

[TABLE]

For optically thin lines, i.e $\tau<1$, the line profiles remain Gaussian, and the integral of their opacity is given by:

[TABLE]

where $\Delta v_{\rm{FWHM}}=\sqrt{8\ln 2}\,\sigma_{v}$ is the full width at half maximum of the observed transition. Thus, Eq. A5 can be rewritten as follows:

[TABLE]

Subtracting Eq. A9 in $C_{\tau}$ allows the construction of a likelihood function $\mathcal{L}(N_{u},T_{{\rm ex}})$, that can be used for least square minimization.

Appendix B Radiative transfer equation used

Assuming a uniform excitation temperature, $T_{{\rm ex}}$, along the line of sight, the spectrum intensity at a given frequency $\nu$, corresponding to the observable source brightness temperature $T_{B}$, is given by the radiative transfer equation as follows:

[TABLE]

where $f$ represents the filling factor, $T_{bg}$ the background temperature, $\tau$ the optical depth and

[TABLE]

the Rayleigh-Jeans equivalent temperature, with $k_{\rm{B}}$ and $h$ the Boltzmann and Planck constants, respectively.

Appendix C Channel maps

In this section we present the channel maps of all the lines presented in the paper.

Appendix D Modeled fractional abundances

(v4.7.2, 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), Nautilus-v1.1 (Hersant et al., 2009; Ruaud et al., 2016; Wakelam et al., 2016)