Sulfur chemistry in protoplanetary disks: CS and H2CS
Romane Le Gal Karin I. Öberg Ryan Loomis Jamila Pegues Jennifer B. Bergner





Romane Le Gal
romane.le_gal@cfa.harvard.edu

0000-0003-1837-3772]Romane Le Gal

0000-0001-8798-1347]Karin I. Öberg
Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA

0000-0002-8932-1219]Ryan Loomis
NRAO, 520 Edgemont Rd, Charlottesville, VA 22903, USA


Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA

0000-0002-8716-0482]Jennifer B. Bergner
Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138, USA


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 ^13CS 6-5 (LkCa 15 and MWC 480), C^34S 6-5 (LkCa 15), and H2CS 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 ^13C and ^34S 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 H2CS by 1-2 orders of magnitude. Finally, comparing the H2CS/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 <cit.>. 
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 <cit.>.


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
<cit.>,
whereas in dense molecular gas it is found highly depleted: only 0.1%
of the sulfur cosmic abundance is observed in the gas phase <cit.>, implying a depletion factor of three orders of magnitude
<cit.>. This level of depletion suggests that most sulfur is locked into icy mantles coating interstellar dust grains <cit.>. However, OCS
and tentatively SO2 are yet the sole molecules detected in icy
grain mantles toward high-mass protostars <cit.> and their abundances (∼-7) make up less than 4% of the sulfur cosmic abundance. Icy H2S,
the most obvious sulfur ice chemistry product from successive
hydrogenations of atomic sulfur in solid phase, has not been detected
yet <cit.>.  


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 <cit.>, protostellar envelopes
<cit.> 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 <cit.> which, could actually even be relics from its parent ISM molecular cloud according to the Rosetta analysis of the 67P/Churyumov-Gerasimenko cometary ices <cit.>. 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 <cit.>. 

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 <cit.>, 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 <cit.>. Finally, very recently, the first detection of H2S was reported toward the T Tauri star GG Tau A <cit.>. 

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, C^34S, and H2CS 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 <ref> 
and the observational results in <ref>, including   
derived column densities and excitation temperatures for
CS and H2CS 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 <ref>. The observational and modeling results are compared and discussed in <ref> and the conclusions presented in <ref>.





§ OBSERVATIONS





 §.§ 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, <cit.>, see Table <ref>). Their ages range from ∼1 to ∼ 10 Myr (see Table <ref>), and host large (>200 au) gas-rich disks <cit.>. 

The sixth PPD is orbiting V4046 Sgr, an isolated binary system located in the 23±3Myr old β-Pictoris Moving Group <cit.>.
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 <cit.>.



 §.§ Observational details



The observed species transitions, their frequencies and spectroscopic parameters are listed in Table <ref>.

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 ∼ 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 H2CS 8_17-7_16, 9_19-8_18 and 9_18-8_17 rotational transitions, CS 6-5 and corresponding ^13CS and C^34S 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 ^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 C^34S 6-5, and H2CS 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 ∼ 19.2 min.

The H2CS 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 ∼ 13.1 min for each block.


The H2CS 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 ∼ 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 () version 4.7.2 <cit.> 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 vconsub function, the data were CLEANed <cit.> 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 <ref>. 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 <cit.> 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. <ref>-<ref>. The disk and star parameters used to build the Keplerian masks are listed Table <ref>.





§ RESULTS





 



 §.§ CS fluxes and emission morphologies



 Figure <ref> 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. <ref>. 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 <ref> 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 <ref>). 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. <ref>, 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. <ref>). 

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 <ref>. 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. <ref> and <ref>), and second out to the furthest radius, R_max, where CS emission is detected. 






 §.§ 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. <ref>). 
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 <ref>). 
The ^12CS 6-5 emission radial morphologies follow those found for ^12CS 5-4 in the respective disk (Fig. <ref>), but the breaks observed at the edge of the continuum are more pronounced.

Figures <ref> and <ref> show that the 6-5 line of the ^13CS and C^34S 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 C^34S 6-5 in the two disks are reported in Table <ref>. 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 <ref>.






 §.§ 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 <cit.>. This approach implies the simplifying assumptions that the lines are optically thin and at local thermal equilibrium (LTE), and are pursued in  <ref>. 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  <ref>. Finally, in  <ref>, we use the constraints on the CS rotational temperature from  <ref> 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]6 and ∼ [2.9-2.2]6 <cit.>, respectively. These are lower than expected in molecular disk layers (see e.g. Fig. <ref>) 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. <ref> and <ref> 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. <cit.>, and  <ref>), 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 <ref>. We then used the likelihood function ℒ(N_u,) described in appendix <ref> with the Python implementation  <cit.> of the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) <cit.> to fit the data and compute posterior probability distributions for  and N_u describing the range of total column densities, N_tot, and excitation temperatures, , consistent with the data. The following uninformative priors were assumed:

    (K) = 𝒰(3,300)
    
        N_tot () = 𝒰(7,20).




Random draws and results from the posterior distributions are depicted in gray in Fig. <ref>, along with the least square fit results, shown in blue. For the MWC 480 disk, the results converged toward {≃22.8±2.5 K, N_tot≃ 7.3±0.912} for a flux integrated over a 150 au radius. For the LkCa 15 disk, the results converged toward
{≃26.8±2.3 K, N_tot≃ 1.2±0.113} 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. <ref>, 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  <ref> assumed the main isotopologue CS lines are optically thin. This assumption can be explored for one disk, LkCa 15, where two CS isotpopologue lines, ^13CS and C^34S 6-5, are well detected.

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

    T_B(^13CS)/T_B(^12CS)  =  1-e^-τ_ ^13CS/1-e^-τ_ ^12CS,

 with T_B the brightness temperature (see Appendix <ref>), 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^-τ)→τ which leads to:

    T_B(^13CS)/T_B(^12CS)  =  τ_^13CS/τ_^12CS ≈  R.

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≃ 68±15. <cit.>.


For optically thin lines we can calculate the disk-integrated brightness temperatures from the disk integrated fluxes S_νΔ v
 using the Rayleigh-Jeans law:

    T_BΔ v = c^2 S_νΔ v/2 k_Bν^2 Ω.


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 <ref> this reduces the calculated flux by ∼25%. 
For the LkCa 15 disk we then obtain:


    T_B(^12CS)/T_B(^13CS) ≃  17 ≪(^12C/^13C)_LISM.


If the ^12CS/^13CS abundance ratio is equal to the ^12C/^13C LISM ratio, Eq. <ref> implies that the ^12CS 6-5 transition is optically thick.
This implication can be used to calculate the ^13CS column density and the ^12CS excitation temperature. For an optically thick line 
(1-e^-τ_^12CS)→ 1. Thus, Eq. <ref> directly provides T_ex(^12CS) from the peak of its brightness temperature:

    T_ex(^12CS)=hν/k_B/ln(1+hν/k_B/T_R(^12CS)/f+J_ν(T_bg)) .


Eq. <ref> gives an excitation temperature of T_ex(^12CS)≃ 7 K which, inserted into <ref>, gives an opacity of τ_^13CS≃ 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^-τ_^13CS)→τ_^13CS, which inserted into Eq. <ref> gives:

    T_B(^13CS)/T_B(^12CS) ≈ τ_^13CS.

Eq. <ref> gives an opacity consistent with the value obtained from Eq. <ref> and <ref>.


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

    N_tot=3h/8π^3Sμ^2Q_rot(T_ex)/g_ue^E_u/k_BT_ex/e^hν/k_B T_ex-1
    ×1/J_ν(T_ex)-J_ν(T_bg)∫ T_R dv   C_τ .

This gives a column density of (^13CS)≃ 712.
Assuming that ^12CS is present at the LISM ratio, 
(^12CS)≃4.814, which is significantly higher than the rotational diagram result, by a factor of forty (see Table <ref>). 
This discrepancy can be understood when considering that τ_^13CS≃ 0.06 thus implies that that τ_^12CS > 1 (assuming LISM ^12C/^13C 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^32S/C^34S ratio considering the elemental 
^32S/^34S= 24.4±5.0 in the vicinity of the Sun <cit.>. 
Both observed ^12CS(6-5)/C^34S(6-5) and ^12CS(5-4)/C^34S(6-5) give a ratio of ∼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 (^12CS)≃ 7.4 K and τ_^34CS = 0.07, leading to (C^34S)≃ 611. Using the LISM ratio we thus obtain (^12CS)≃ 1.513, 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 <ref> and <ref> to derive CS column densities and excitation temperatures are summarized in Table <ref>. 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 <cit.> 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  <ref>).
Although Figure <ref> shows apparent asymmetries for the ^13CS and C^34S 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 C^34S may not be present at the cosmic ratios. There are many examples of isotopic fractionation in disks and the ISM, including in carbon <cit.>. 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 C^34S 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  <ref> 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. <ref> and <ref> to calculate the column density.
The results are presented in red in Fig. <ref>.
The CS disk-integrated column density varies by less than an order of magnitude across the sample of disks. The average column density is ≈ 712 and the standard deviation ≈ 211.
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 <cit.>. For colder disk temperatures the derived CS column densities increase, as expected from Eq. <ref>, leading for instance to an average column density of ≈ 513 with a standard deviation of ≈± 512 for a disk temperature of 10 ± 2 K, i.e. a factor 5 higher column density.









 §.§ H2CS in MWC 480 and LkCa 15








The three H2CS 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. <ref> and the radial profiles in Fig. <ref>. 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  <ref> for the ^13CS and C^34S emissions, the quality of the H2CS 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.


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 H2CS lines presented here are ortho lines, which can be used to derive a column density of the o-H_2CS, and a total column density assuming an o/p ratio. This is done for MWC 480 in Fig. <ref> using the rotational diagram method described Sect. <ref>. 
The following flat priors were assumed:

    (K) = 𝒰(10,100),
    
        N_tot () = 𝒰(7,20).



Random draws and results from the posterior distributions are depicted in gray in Fig. <ref>, along with the least square fit results, shown in blue. The results converged toward 
{≃41^+38_-20 K, N_tot≃ 3.0^+2.7_-0.412}, 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. <ref> and <ref> to derived an upper limit of ≲212 when considering the three H2CS 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.











§ 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  <cit.>, tuned to the physical structures of MWC 480 and LkCa 15.




 §.§ 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 <ref>. We use the previously calculated temperature and density profiles from <cit.> and <cit.>,assuming that the temperature and surface density, governing line emissions, vary as power laws of the radii <cit.>.
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 <cit.> used in <cit.> and <cit.>:


    T(z) = {[ T_mid+(T_atm-T_mid)[ sin(π z/2z_q) ]^2δ                                   z<z_q;                                   T_atm                                 z≥ z_q, ].

where 

    T_mid=T_mid,R_c (r/R_c)^-q,
    
    T_atm=T_atm,R_c (r/R_c)^-q,

vary as power law of the radii <cit.>. 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, , becomes ≲ 1. H, defined by the midplane temperature and assuming vertical static equilibrium, can be expressed as follows:


    H=√(%s/%s)k_B  T_mid r^3μ m_H  G  M_⋆,

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,R_c can be estimated from the following approximation for a simple irradiated passive flared disk:

    T_mid(r)≈(φ L_⋆/8π r^2 σ_SB)^1/4,

with L_⋆ the stellar luminosity, σ_SB the Stefan-Boltzman constant and φ the flaring angle (e.g., <cit.>.

The dust temperature is assumed equal to that of the gas. This is not a valid assumption in the uppermost layers of the disk <cit.>, but is a reasonable simplification in this study focusing on S-chemistry in the 
molecular disk layers <cit.>. 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 <cit.>. The resulting 2D temperature profiles for the MWC 480 and LkCa 15 disks are shown in the first column panels of Fig. <ref>.




  §.§.§ Density profile





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

    δlnρ(z)/δ z=-[G M_⋆ z/r^3μ  m_H/k_B T + δln T/δ z],

which gives for the isothermal case:

    δρ(z)/ρ(z)=-[G M_⋆ z/r^3μ  m_H/k_B T_mid]δ z,

with solution:

    ρ(z)=ρ_oexp(-z^2/2H^2),

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:

    Σ=ρ_0 ∫_-∞^+∞exp(-z^2/2H^2)=√(2π)ρ_0H.


Furthermore, the surface density of the disk is assumed to follow a simple power law varying as r^-3/2 <cit.> with sharp cutoff at specified inner and outer radii (i.e. R_int< r <R_out):

    Σ (r) = Σ_c ( r/R_c)^-3/2,

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:

    2π r δΣ/δ t=-Ṁ/δ r,

If the mass loss of the disk is constant with time, i.e. δΣ/δ t=0, the disk mass can simply be expressed as follows:

    M_disk=∫_R_int^R_out 2 π r Σ(r) dr.


Thus, the surface density at a radius R reads:

    Σ_R=M_disk R^-3/2/4π√(R_out),

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. <ref>, 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.621 <cit.>, with N_H=N(H)+2N(H2) 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. <ref>, 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 <cit.>'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) <cit.>. 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 <cit.>. 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,R_c (i.e. the UV flux factor at a radius R_c, see Table <ref>), is divided by two <cit.>. 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:

    f_UV=f_UV,R_c/2/(r/R_c)^2+(4H/R_c)^2.



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. <ref>, for each model.






















































 §.§ Chemical modeling



Once built, the 
1+1D physical structure described above (and depicted in Fig. <ref> for the MWC 480 and LkCa 15 cases) is fed as an input in the pseudo 

time-dependent astrochemical modeling code  <cit.> used in three-phase mode <cit.>. This rate-equation gas-grain code – based on <cit.> and <cit.> – 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 <cit.>, and non thermal ones, such as cosmic-ray induced desorption <cit.>, photodesorption and chemical desorption <cit.>. For photodesorption we use the prescriptions discussed in <cit.>, where we followed the recommendations of <cit.> 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-4 molecule/photon <cit.> 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 μ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 <cit.>, 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 16 years <cit.>. 
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(H2)=24 and a cosmic ionization rate of 1.3-17; 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 <ref>)  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 <ref>. 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 <ref>), 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 <cit.>. 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. <ref>). 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-8 <cit.>,
 i.e. ∼2 orders of magnitude below the cosmic value ∼1.5-5 <cit.>. 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 <cit.>. For the present study, this approximation seems reasonable since CS and H2CS are supposed to reside in disk molecular layers, which are similar to dense astrophysical molecular environments.










 §.§ Model results








The resulting 
number densities of CS, H2CS and their main precursors, S+ and S, 
at 1 Myr 
are presented in Fig. <ref> for the MWC 480 and LkCa 15 disk models (the corresponding fractional abundances

are shown in Appendix <ref>). 
Fig. <ref> 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 <ref>). 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. <ref>), 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. <ref> with Fig. <ref>, we see that the layers where CS and H_2CS reside have similar densities of ∼6-8.


In both disks the S and H2CS 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_2CS 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_yH (with x=1-4 and y=2-3). This produces carbonated S-ions including HCS+, CS+, HC3S+, and C2S+, 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+, H3+, HCO+) and with He+, and freeze-out onto grain surfaces.

H2CS is mainly produced by the gas-phase neutral-neutral reaction of atomic S and CH3 and the electronic dissociative recombination of H3CS+, whose formation begins with the reaction between S+ and CH4. 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-H2CS releasing some (∼ 1%) of its H2CS product in the gas-phase by chemical reactive desorption. Similarly to CS, H2CS is mainly destroyed by photodissociation processes in the disk atmosphere, reactions with protons and protonated ions (i.e. with H+, H3+, 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. <ref>. 

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 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  <ref>). Running the same models with undepleted S-element abundance results in an overprediction of CS by 1-2 orders of magnitude, while H2CS 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_2CS column densities are 1–2 orders of magnitude below the predicted CS column densities in each disks.
Similarly to CS, the H_2CS 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_2CS precursors, with the possible exception of S toward LkCa 15. 


Compared with observations toward MWC 480, the model underpredicts H_2CS by several orders of magnitude. Furthermore the model predicts that more H_2CS should be present toward LkCa 15 compared to MWC 480, which is opposite to what we observed. Finally the model predicts that H_2CS should be centrally very peaked, while observations suggest more extended emission.






Figure <ref> presents the calculated mean disk gas temperature weighted by the density of CS and H_2CS, 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 H2CS. While the H_2CS 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. <ref>). The rotational temperatures of CS and H2CS are also overplotted on Fig. <ref>. 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.







§ DISCUSSION










 §.§ Model predictions vs observations









In  <ref> we compared model predictions of CS and H_2CS 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 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_2CS in the outer disk of MWC 480, and incorrectly predict that H_2CS 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 <cit.> 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 H_2CS is that grain surface formation of H_2CS and subsequent photodesorption may be more important in disks than is currently implied by the model. The surface and bulk ice H_2CS 2D disk densities are represented in Fig. <ref>, 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 CH3OH ices showed an efficient dissociation of the product rather than a desorption <cit.>. H2CS 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 H2CS ice abundance and desorption, unless additional ice formation pathways are missing in our astrochemical networks. Therefore, both, photo-irradiation of icy H2CS and new H2CS 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 <cit.> and the importance of ice photodesorption in this region <cit.>. 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_2CS 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 H2CS 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 H2CS observations and its predictions from the astrochemical disk modeling. Thus, higher H2CS 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. <ref>), 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. <ref> 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  <ref>, the MWC 480 disk being more massive, its visual extinction is larger (see Fig. <ref>), implying less UV penetration in the disk. In addition, the dramatic changes in temperature over small radii shown in Fig.<ref> 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 <cit.> such as the abundances of S+, CS+, HCS+ as well as the SO/SO2 formation, and thus affects the overall S-chemistry <cit.>. However, for the ALMA observations presented here, concerning disk radii beyond ∼10-20  au, X-rays chemistry is not dominant <cit.>. 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 <ref>) are consistent with the values found in other disks, e.g. 2.0±0.1612 in GO Tau and 8.7±1.612 in LkCa 15 <cit.>,
6±312 in DM Tau <cit.> and 2.213 in GG Tau <cit.>, which agrees with CS being the most abundant S-species detected in disks so far.




In order to place our CS and H2CS 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 <cit.> and cold dense core <cit.>, analogs of the main three chemical layers constituting protoplanetary disks, i.e. the disk atmosphere, intermediate molecular layers, and the disk midplane. Figure <ref> 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 (∼ 0.06) than those measured so far in the different ISM environments presented in Fig. <ref> (∼ 0.01-0.03). However, the wide error bars we derived (± 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 C^34S/C^32S 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^34S/C^32S 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 H2CS/CS ratio decreasing from dense to diffuse ISM, is found in disks (i.e. in LkCa 15) almost as high as in dense cores <cit.>. Nevertheless, the high uncertainties on our measurement do not allow us to firmly assert that the H2CS/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 H2CS/CS ratio, the H2S/CS is increasing from dense to diffuse ISM and the disk value <cit.> is found almost as low as the lower limit derived in dense core <cit.>. This could also explain why H2S was so hard to detect in disks <cit.>, aside from the fact that the H2S 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, H2S 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 <cit.>. However, and in opposition to the H2S/CS ratio, the OCS/CS ratio is found lower in PDR than in dense core, suggesting that the OCS and H2S molecules trace different type of environments. 


As the H2S/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 <cit.>. <cit.> 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 <cit.>. Contrary to the SO/CS, the SO2/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 <cit.> do not allow us to settle any specific trend for this ratio. 

Finally, the C2S/CS ratio decreases from dense to diffuse ISM, in an even more pronounced way than the OCS/CS and H2S/CS ratios do. In this case, the disk upper limit <cit.> found in between the value from the dense prestellar core L1544 <cit.> 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 <ref> 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. 













§ 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). 
    
  * 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 (CS)≃ (7.3±0.8)12, ≃23±2.5 K and (CS)=(1.2±0.1)13, =27±2.3 K, assuming LTE and optically thin emission.
    
  * We used the average CS rotational temperature derived toward MWC 480 and LkCa 15, i.e. ≃25±5 K, to estimate the CS column densities toward our full disk sample, resulting in an average CS column density of (CS)≈712, with a standard deviation of ≈211.
    
  * We report the first detections in disks of ^13CS, C^34S (in LkCa 15), and H2CS (in MWC 480). 
    
  * 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.
    
  * From the three H2CS line detections obtained toward MWC 480 we derived an H2CS column density and rotational temperature of (H2CS)≃ 3 12, ≃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.
    
    
  * 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. 
    
  * We also compared our H2CS observations to 2D disk astrochemical modeling. For that case, our model underpredict H2CS 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.
    
  * 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.
    
  * 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 H2CS column densities and line emission regions.
Non-LTE modeling would also be helpful to better constrain the CS and H2CS 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 <cit.>. 

 
 

  










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







§ BUILDING ROTATIONAL DIAGRAMS FROM DISK-INTEGRATED FLUX DENSITIES



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

    N_u= 4π S_νΔ v/A_ulΩ hc,

where S_ν is the flux density, Δ v the line width, A_ul the Einstein coefficient, c the speed of light and
Ω the solid angle subtended by the source <cit.>. For this analysis we use the disk flux densities S_νΔ v integrated out to 150 au and 200 au in the MWC 480 and LkCa 15 disks respectively, referred to R_σ in Table <ref>.

The total column density, N_tot, and rotational temperature, T, can then be derived from the upper level population, N_u, following the Boltzmann distribution:

    N_u/g_u= N_tot/Q_rot(T)e^-E_u/k_BT,

with g_u and E_u the degeneracy and energy of the upper energy level u, respectively; k_B the Boltzmann constant; and Q_rot the partition function of the molecule, which for a diatomic molecule as CS can be approximated by:

    Q_rot(T) ≈k_B T/h B_0+1/3,

<cit.>. 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. <ref>, as depicted in blue in Fig. <ref>, 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_τ =τ/1-e^-τ, in case τ ≪ 1 <cit.>.
The resulting logarithm of Eq. <ref> is then: 

    lnN_u/g_u= lnN_tot-lnQ_rot(T)-lnC_τ-E_u/k_BT.


In the LTE approximation, the optical depth of a given transition at  can be expressed as a function of its upper energy level population using Eq. <ref> as follows:

    ∫τ_ν dv =  N_u A_ul  c^3/8πν^3(e^hν/k_B -1).


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

    ∫τ_ν dv = 
        τ_νΔ v_FWHM√(π/4ln2),

where Δ v_FWHM=√(8ln 2) σ_v is the full width at half maximum of the observed transition. Thus, Eq. <ref> can be rewritten as follows:

    τ_ν = √(4ln2/π) N_u A_ul  c^3/Δ v_FWHM  8πν^3(e^hν/k_B -1).



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

    ∫τ_ν dv = 
        τ_νΔ v_FWHM√(π/4ln2),

where Δ v_FWHM=√(8ln 2) σ_v is the full width at half maximum of the observed transition. Thus, Eq. <ref> can be rewritten as follows:

    τ_ν = √(4ln2/π) N_u A_ul  c^3/Δ v_FWHM  8πν^3(e^hν/k_B -1).


Subtracting Eq. <ref> in C_τ allows 
the construction of a likelihood function ℒ(N_u,), that can be used for least square minimization. 



§ RADIATIVE TRANSFER EQUATION USED


Assuming  a  uniform  excitation  temperature, , 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:


    T_B  =   f   [J_ν(T_ex)-J_ν(T_bg)]   [1-e^-τ_ν],

where f represents the filling factor, T_bg the background temperature, τ the optical depth and

    J_ν(T)=hν/k_B/e^hν/k_BT-1,

the Rayleigh-Jeans equivalent temperature, with k_B and h the Boltzmann and Planck constants, respectively.





§ CHANNEL MAPS

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










































§ MODELED FRACTIONAL ABUNDANCES














 <cit.>,  <cit.>,  <cit.>,  <cit.>,   <cit.>,   <cit.>,  <cit.>,  <cit.>























