Sulfur chemistry in protoplanetary disks: CS and H2CS

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.


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;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 SO 2 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 (∼ 10 −7 ) make up less than 4% of the sulfur cosmic abundance. Icy H 2 S, 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. 2013Guilloteau et al. , 2016Pacheco-Vázquez et al. 2016;Sakai et al. 2016;Booth et al. 2018). Finally, very recently, the first detection of H 2 S 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 13 CS, C 34 S, and H 2 CS 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 H 2 CS 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.
The sixth PPD is orbiting V4046 Sgr, an isolated binary system located in the 23 ± 3Myr old β-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).

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. Oberg), 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 ∼ 0.5 . 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 H 2 CS 8 17 − 7 16 , 9 19 − 8 18 and 9 18 − 8 17 rotational transitions, CS 6 − 5 and corresponding 13 CS and C 34 S 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 ∼ 975 kHz, i.e. corresponding to a velocity resolution of ∼ 1 km/s. The 12 CS 6 − 5 transition was observed in three execution blocks on January 17, 2016, April 23, 2016and 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 13 CS and C 34 S 6−5, and H 2 CS 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 12 CS 6 − 5 (see above) but with 36 antennas and a total on-source integration time of ∼ 19.2 min.
The H 2 CS 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 12 CS 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 ∼ 13.1 min for each block.
The H 2 CS 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 12 CS 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 ∼ 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 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.

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. Selfcalibrations 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σ 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 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 ∼ 1.3 mm continuum emission maps. The angular resolution is ∼ 0.5 , which corresponds to ∼ 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 ∼ 1 mm continuum emission maps toward V4046 Sgr at an angular resolution of ∼ 0.5 , which corresponds to ∼ 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.

CS fluxes and emission morphologies
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 ∼ 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 ν ∆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 σ , where the CS flux is ∼ 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.  a https://cdms.astro.uni-koeln.de/cdms/portal/, Müller et al. (2001Müller et al. ( , 2005 b Rσ corresponds to the radius at which the line flux is ∼ 33% of the maximum line emission peak, thus focusing on the disk from where most of the CS emission (i.e. 2/3) originates when deriving a characteristic CS column density.  [10,30,100,200,400,800,1500,1900]× the median RMS and [3,6,9,12,15]× the median RMS, respectively. Synthesized beams are shown in the lower left of each panel. Top right panel: The black cross indicates the center of the continuum. Second row: Radially deprojected and azimuthally averaged flux profiles of the continuum and the 12 CS 6 − 5 emission. Third row: Integrated intensity spectra of 12 CS 6 − 5 emission.

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 12 CS 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 12 CS 6−5 emission radial morphologies follow those found for 12 CS 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 13 CS and C 34 S 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, 13 CS and C 34 S 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 13 CS and the main isotopologue is quite high, which is discussed in detail in §3.3.2.

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 12 CS 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.

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 ∼ [1.7 − 1.3] ×10 6 and ∼ [2.9 − 2.2] ×10 6 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 ex = T 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 ∼ 0.6 − 0.8 * T kin (unless optical depths are extreme), meaning that a 9 K brightness temperature line could partly be optically thick if T kin ∼ 11 − 15 K. But, in particular for MWC 480, the lowest expected gas temperature in CS emitting regions is ∼ 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 L(N u , T 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 ex and N u describing the range of total column densities, N tot , and excitation temperatures, T ex , consistent with the data. The following uninformative priors were assumed: 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 ex 22.8±2.5 K, N tot 7.3 ± 0.9 ×10 12 cm −2 } for a flux integrated over a 150 au radius. For the LkCa 15 disk, the results converged toward {T ex 26.8 ± 2.3 K, N tot 1.2 ± 0.1 ×10 13 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 τ = 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.

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 CS and C 34 S 6 − 5, are well detected.  The ratio of the two isotopologues 12 CS and 13 CS can be expressed as: 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 13 CS 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 12 CS and 13 CS are both optically thin then we have (1 − e −τ ) → τ which leads to: where R is the abundance ratio of 12 CS to 13 CS. Assuming that their abundance ratio is equal to the 12 C/ 13 C ratio determined in the local interstellar medium (LISM), i.e. that their has been no isotopic fractionation, R 68 ± 15. (Milam et al. 2005;Asplund et al. 2009;Manfroid et al. 2009).
For optically thin lines we can calculate the diskintegrated brightness temperatures from the disk integrated fluxes S ν ∆v using the Rayleigh-Jeans law: To consistently calculate disk integrated fluxes for the main and rare CS isotopologues we applied a uv-taper to the extended baselines of the 12 CS visibilities before integrating the fluxes, to match the synthesized beams. As can be seen in Table 2 this reduces the calculated flux by ∼25%. For the LkCa 15 disk we then obtain: If the 12 CS/ 13 CS abundance ratio is equal to the 12 C/ 13 C LISM ratio, Eq. 6 implies that the 12 CS 6 − 5 transition is optically thick. This implication can be used to calculate the 13 CS column density and the 12 CS excitation temperature. For an optically thick line (1 − e −τ12 CS ) → 1. Thus, Eq. B10 directly provides T ex ( 12 CS) from the peak of its brightness temperature: Eq. 7 gives an excitation temperature of T ex ( 12 CS) 7 K which, inserted into B10, gives an opacity of τ13 CS 0.06 (this last step assumes that CS and 13 CS have the same excitation temperatures). Moreover, since we assume that 13 CS is optically thin we have (1−e −τ13 CS ) → τ13 CS , which inserted into Eq. 3 gives: Eq. 8 gives an opacity consistent with the value obtained from Eq. 7 and B10.
Once the opacity and excitation temperature derived, the column density can be computed as follows: This gives a column density of N tot ( 13 CS) 7 ×10 12 cm −2 .
Assuming that 12 CS is present at the LISM ratio, N tot ( 12 CS) 4.8 ×10 14 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 τ13 CS 0.06 thus implies that that τ12 CS > 1 (assuming LISM 12 C/ 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 C 32 S/C 34 S ratio considering the elemental 32 S/ 34 S= 24.4 ± 5.0 in the vicinity of the Sun (Chin et al. 1996). Both observed 12 CS(6−5)/C 34 S(6−5) and 12 CS(5−4)/C 34 S(6−5) give a ratio of ∼13, which is somewhat lower than the cosmic 32 S/ 34 S ratio, but not nearly as much below the cosmic ratio as 12 CS/ 13 CS. Taken at face value, the excitation temperature and opacity are T ex ( 12 CS) 7.4 K and τ34 CS = 0.07, leading to N tot (C 34 S) 6 ×10 11 cm −2 . Using the LISM ratio we thus obtain N tot ( 12 CS) 1.5 ×10 13 cm −2 , which is an order of magnitude lower than the value derived from the 12 CS/ 13 CS 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 12 CS line is optically thick, but the large difference in derived 12 CS 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 13 CS and C 34 S 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 morphologi-  Table 2) 10% calibration uncertainty and 20% uncertainty on the temperature are included.
cal differences between CS and its isotopologues. However, differences in beam-filling factors of the different isotopologues could impact this analysis. Higher signalto-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, 13 CS and C 34 S 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. Sakai et al. 2010Sakai et al. , 2013Ossenkopf et al. 2013;Smith et al. 2015;Yoshida et al. 2015;Yu et al. 2016;Huang et al. 2017). For instance, the 13 CS enhanced abundances in disk molecular layers could be explained by excessive 13 C-carrier, such as 13 CO (while 12 CO remains self-shielded). Rapid ion-molecule gas-phase reactions could thus transfer 13 C from 13 CO to hydrocarbons, and then to 13 CS. The degree of fractionation that would be needed to explain the 13 CS 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 C 34 S data within error bars, for model comparisons and calculations of column densities in disks without constraints on the excitation temperature.

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 ∼ 23 ± 3, and ∼ 27 ± 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 ±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 ≈ 7 ×10 12 cm −2 and the standard deviation ≈ 2 ×10 11 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 ≈ 5 ×10 13 cm −2 with a standard deviation of ≈ ±5 ×10 12 cm −2 for a disk temperature of 10 ± 2 K, i.e. a factor 5 higher column density.

H 2 CS in MWC 480 and LkCa 15
The three H 2 CS 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σ). The integrated flux density maps are shown in Fig. 8 and the radial profiles in Fig. 9. The emission toward MWC 480 shows no substructure, 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 13 CS and C 34 S emissions, the quality of the H 2 CS 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. H 2 CS 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 H 2 CS lines presented here are ortho lines, which can be used to derive a column density of the o-H 2 CS, 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: N tot ( cm −2 ) = U( 10 7 , 10 20 ).
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 ex 41 +38 −20 K, N tot 3.0 +2.7 −0.4 ×10 12 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 tot 2 ×10 12 cm −2 when considering the three H 2 CS line tentatively detected toward LkCa 15, a rotational temperature of T ≈ 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.

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.

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): where 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 V , becomes 1. H, defined by the midplane temperature and assuming vertical static equilibrium, can be expressed as follows: with k B the Boltzmann constant, µ = 2.4 the reduced mass of the gas, m H the proton mass, G the gravitational constant, and M the mass of the central star. The midplane temperature T mid,Rc can be estimated from the following approximation for a simple irradiated passive flared disk: with L the stellar luminosity, σ SB the Stefan-Boltzman constant and ϕ 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.

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: which gives for the isothermal case: with solution: where ρ o corresponds to the volume density of the gas in the midplane. The surface density of the disk, Σ, can thus be expressed as follows: 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 int < r < R out ): where Σ 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Ṁ = −2πrΣ(r)v r , with v r the radial velocity of the disk. For axis-symmetric disk, mass conservation thus reads: If the mass loss of the disk is constant with time, i.e. δΣ δt = 0, the disk mass can simply be expressed as follows: Thus, the surface density at a radius R reads: which highlights the link between the mass of the disk, M disk , and its outer radius, R out . The resulting 2D density profile is represented in the second column panels of Fig. 11, for each model.

Visual extinction profile
The vertical visual extinction gradient is computed from the hydrostatic density profile using the conversion factor N H /A V = 1.6 ×10 21 (Wagenblast & Hartquist 1989), with N H = N (H) + 2N (H 2 ) the vertical hydrogen column density of hydrogen nuclei. This conversion factor assumes a typical mean grain radius size of 0.1 µ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.

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 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 UV,Rc (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: 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.

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 timedependent astrochemical modeling code Nautilus v1.1 ) used in three-phase mode ). This rate-equation gas-grain code -based on Hasegawa et al. (1992) and Hasegawa & Herbst (1993) -simulates the time-dependent chemistry of ∼ 1100 species (half in the gas phase and half in solid phase) linked together via more than ∼ 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 ×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 neutralneutral 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 µ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 ×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 H = n(H) + n(H 2 ) = 2 ×10 4 cm −3 and a cosmic ionization rate of 1.3 ×10 −17 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 5) 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 5. 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   Note-(1) Wakelam & Herbst (2008); (2)  model described above and considers a "low" sulfur elemental abundance of 8 ×10 −8 (see Table 5, Graedel et al. 1982;Wakelam & Herbst 2008), i.e. ∼ 2 orders of magnitude below the cosmic value ∼ 1.5 ×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 H 2 CS are supposed to reside in disk molecular layers, which are similar to dense astrophysical molecular environments.

Model results
The resulting number densities of CS, H 2 CS and their main precursors, S + 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 H 2 CS reside have similar densities of ∼ 10 6 − 10 8 cm −3 .
In both disks the S and H 2 CS layers reside deeper toward the disk midplane than S + and CS. The depth of the S/S + transition is set by a balance of, on the one hand, production of S + from S through photo-ionization and charge transfer reaction with C + and, on the other hand, recombination of S + with electrons to form S. The main formation pathway to form CS starts with S + , while H 2 CS 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 S + and small hydrocarbons CH x and C y H (with x = 1−4 and y = 2−3). This produces carbonated S-ions including HCS + , CS + , HC 3 S + , and C 2 S + , that subsequently recombine with electrons to form neutral S-containing species. CS is thus a daughter molecule of S + . 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 H + , H + 3 , HCO + ) and with He + , and freeze-out onto grain surfaces.
H 2 CS is mainly produced by the gas-phase neutralneutral reaction of atomic S and CH 3 and the electronic dissociative recombination of H 3 CS + , whose formation begins with the reaction between S + and CH 4 . 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 s−H + s−HCS −−→ s−H 2 CS releasing some (∼ 1%) of its H 2 CS product in the gas-phase by chemical reactive desorption. Similarly to CS, H 2 CS is mainly destroyed by photodissociation processes in the disk atmosphere, reactions with protons and protonated ions (i.e. with H + , H + 3 , HCO + ), 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 ∼ 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 ∼ 2 − 10, consistent with ob- servations. 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 13 CS 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 H 2 CS 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 H 2 CS column densities are 1-2 orders of magnitude below the predicted CS column densities in each disks. Similarly to CS, the H 2 CS 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 H 2 CS precursors, with the possible exception of S toward LkCa 15. Compared with observations toward MWC 480, the model underpredicts H 2 CS by several orders of magnitude. Furthermore the model predicts that more H 2 CS should be present toward LkCa 15 compared to MWC 480, which is opposite to what we observed. Finally the model predicts that H 2 CS 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 H 2 CS, 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 ∼ 8 − 12 K and ∼ 20 − 30 K, in the LkCa 15 and MWC 480 disk models, respectively. In the inner disk (i.e. for radius 50 au) they reside at temperatures in the ranges ∼ 25−200 K and ∼ 60−400 K, respectively. Similar behaviors can be seen for H 2 CS. While the H 2 CS 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 H 2 CS are also overplotted on Fig. 14. These temperatures were derived from observations out to ∼190 au (MWC 480) amd ∼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.

Model predictions vs observations
In § 4.3 we compared model predictions of CS and H 2 CS 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 10 12 -10 13 cm −2 are well reproduced, within a factor of a few, by our models. The best fit are found for radii in the range ∼ 150 − 250 au for MWC 480, and for radii of ∼ 100 au and in the range ∼ 200−275 au for LkCa 15. However, we should note the discrepancies obtained for the inner disk of MWC 480 (i.e. radii 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 H 2 CS in the outer disk of MWC 480, and incorrectly predict that H 2 CS 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 modeldata comparisons in the inner disk. A possible reason for the overall underprediction of H 2 CS is that grain surface formation of H 2 CS and subsequent photodesorption may be more important in disks than is currently implied by the model. The surface and bulk ice H 2 CS 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 CH 3 OH ices showed an efficient dissociation of the product rather than a desorption (Bertin et al. 2016;Cruz-Diaz et al. 2016). H 2 CS 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 H 2 CS ice abundance and desorption, unless additional ice formation pathways are missing in our astrochemical networks. Therefore, both, photo-irradiation of icy H 2 CS and new H 2 CS 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 H 2 CS 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 H 2 CS 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 H 2 CS observations and its predictions from the astrochemical disk modeling. Thus, higher H 2 CS 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 gasphase; 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 theoryfocused 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. Xrays 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 abun-dances of S + , CS + , HCS + as well as the SO/SO 2 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 ∼ 10 − 20 au, X-rays chemistry is not dominant (e.g. Rab et al. 2018). Nevertheless, this would be an interesting future modeling improvement.

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 ± 0.16 ×10 12 cm −2 in GO Tau and 8.7 ± 1.6 ×10 12 cm −2 in LkCa 15 ), 6 ± 3 ×10 12 cm −2 in DM Tau (Semenov et al. 2018) and 2.2 ×10 13 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 H 2 CS 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 13 CS/ 12 CS value we derived in disks from this study is found higher (∼ 0.06) than those measured so far in the different ISM environments presented in Fig. 16 (∼ 0.01 − 0.03). However, the wide error bars we derived (±0.05) do not allow us to draw strong conclusions on the plausible 13 C enhancement suggested in this work. Higher resolution observations of CS isotopolgues in disks are needed to further investigate this 13 C enhancement hint.
The C 34 S/C 32 S ratio we observed in disks (∼ 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 C 34 S/C 32 S 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 Schemistry, 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 H 2 CS/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 H 2 CS/CS ratio behaves in disk more like dense core S-chemistry than diffuse ISM Schemistry. To arbitrate on that, higher resolution observations are also needed.
Interestingly, and in opposition to the H 2 CS/CS ratio, the H 2 S/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 H 2 S was so hard to detect in disks (e.g. Dutrey et al. 2011), aside from the fact that the H 2 S 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, H 2 S 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 H 2 S/CS ratio, the OCS/CS ratio is found lower in PDR than in dense core, suggesting that the OCS and H 2 S molecules trace different type of environments.
As the H 2 S/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 SO 2 /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 C 2 S/CS ratio decreases from dense to diffuse ISM, in an even more pronounced way than the OCS/CS and H 2 S/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.
3. We used the average CS rotational temperature derived toward MWC 480 and LkCa 15, i.e. T ex 25 ± 5 K, to estimate the CS column densities toward our full disk sample, resulting in an average CS column density of N tot (CS)≈ 7 ×10 12 , with a standard deviation of ≈ 2 ×10 11 . 4. We report the first detections in disks of 13 CS, C 34 S (in LkCa 15), and H 2 CS (in MWC 480).
5. We used the CS isotopologue detections toward LkCa 15 to assess the CS disk opacity and tentatively conclude that either both 13 C and 34 S 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. From the three H 2 CS line detections obtained toward MWC 480 we derived an H 2 CS column density and rotational temperature of N tot (H 2 CS) 3 ×10 12 , T ex 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. 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. We also compared our H 2 CS observations to 2D disk astrochemical modeling. For that case, our model underpredict H 2 CS 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 Sdisk theoretical investigations.
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. 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 protoplan-etary 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 H 2 CS column densities and line emission regions. Non-LTE modeling would also be helpful to better constrain the CS and H 2 CS 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).

B. RADIATIVE TRANSFER EQUATION USED
Assuming a uniform excitation temperature, T ex , along the line of sight, the spectrum intensity at a given frequency ν, corresponding to the observable source brightness temperature T B , is given by the radiative transfer equation as follows: where f represents the filling factor, T bg the background temperature, τ the optical depth and J ν (T ) = hν/k B e hν/kBT − 1 , the Rayleigh-Jeans equivalent temperature, with k B and h the Boltzmann and Planck constants, respectively.

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