Sulfur Chemistry in Protoplanetary Disks: CS and H₂CS

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

The sixth PPD is orbiting V4046 Sgr, an isolated binary system located in the  23 ± 3 Myr 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 (Piétu et al. 2006; Andrews et al. 2009; Rosenfeld et al. 2013).

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

Table 2.  List of Observations (Data from CDMS^a)

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

Notes.

^a https://cdms.astro.uni-koeln.de/cdms/portal/; Müller et al. (2001, 2005). ^bR_σ 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.

Download table as:  ASCIITypeset image

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 2016 December 1 with an angular resolution of ∼05. The total on-source integration times were 22–24 minutes. A total of 44 antennas were included and covered baselines from 15 to 704 m. The first block used the source J0510+1800 as its bandpass calibrator, while the second block used J0237+2848. Both blocks used sources J0423–0120 and J0426+2327 as a flux calibrator and phase calibrator, respectively.

The three H₂CS 8₁₇–7₁₆, 9₁₉–8₁₈, and 9₁₈–8₁₇ rotational transitions, CS 6–5, and corresponding ¹³CS and C³⁴S isotopologues were observed in Band 7 during ALMA Cycles 3 and 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, corresponding to a velocity resolution of ∼1 km s⁻¹.

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

The ¹³CS, C³⁴S 6–5, and H₂CS 8₁₇–7₁₆ rotational transitions were observed in a single execution block, on 2016 January 17, with the same calibrator sources as those used for the first and second execution blocks used for ¹²CS 6–5 (see above) but with 36 antennas and a total on-source integration time of ∼19.2 minutes.

The H₂CS 9₁₉–8₁₈ rotational transition was observed in two execution blocks, on 2016 December 17 and 18, with the same calibrator sources as those used for the first and second execution blocks used for ¹²CS 6–5 (see above). A total of 40 and 42 antennas were used with baseline coverages from 15 to 460 m and from 15 to 492 m, respectively. The total on-source integration times were ∼13.1 minutes for each block.

The H₂CS 9₁₈–8₁₇ rotational transition was observed in two execution blocks, on 2016 December 13 and 14, with the same calibrator sources as those used for the first and second execution blocks used for ¹²CS 6–5 (see above). A total of 36 and 39 antennas were used with baseline coverages of 15–650 m and 15–460 m, respectively. The total on-source integration times were ∼12.6 minutes for each block.

More details on the ALMA Band 7 observations presented above can be found in R. A. Loomis et al. (2019, in preparation).

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 2018 May 6 with an angular resolution of 054 and a spectral resolution of 0.1 km s⁻¹. The total on-source integration time was 29 minutes. A total of 44 antennas were included and covered baselines from 15 to 500 m. The quasar J1924–2914 was used as a bandpass and flux calibrator, and J1802–3940 was used as a phase calibrator.

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

After continuum subtraction with the CASA uvconsub function, the data were CLEANed (Högbom 1974) using a 3σ 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 (S/N) 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; see details in the Appendix of J. Pegues et al. 2019, in preparation, for our methodology) to the CLEANed data cubes to improve the S/N 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 (J. Pegues et al. 2019, in preparation), which is originating in at least as high layers as our molecules owing to high optical depths. The Keplerian masks calculated for the different lines studied here are depicted on the channel maps in the figures shown in Appendix A. The disk and star parameters used to build the Keplerian masks are listed Table 1.

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 ∼05, which corresponds to ∼80 au. The spectra of the CS 5–4 are also depicted in Figure 1. For four of our targets, there is a typical double-peaked profile tracing the Keplerian rotation of the disk, while the double-peaked signature is not obvious toward the faintest PPD, DO Tau. This could be due to the 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 us better constrain its physical structure. Figure 2 shows the CS 6–5 line and the ∼1 mm continuum emission maps toward V4046 Sgr at an angular resolution of ∼05, which corresponds to ∼32 au (see source distance in Table 1). As shown in the last panel of the figure, the spectrum presents a typical double-peaked profile corresponding to Keplerian rotation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 in the second and third rows of Figure 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 Figure 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 in 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 Figures 1 and 2), and second out to the farthest radius, R_max, where CS emission is detected.

Toward two disks in the sample, MWC 480 and LkCa 15, three additional CS lines were observed (Figure 3). The ¹²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 ¹²CS 6–5 emission radial morphologies follow those found for ¹²CS 5–4 in the respective disk (Figure 4), but the breaks observed at the edge of the continuum are more pronounced.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 3 and 4 show that the 6–5 lines of the ¹³CS and C³⁴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 emissions of CS, ¹³CS, and C³⁴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 ¹³CS and the main isotopologue is quite high, which is discussed in detail in Section 3.3.2.

Emission-line fluxes encode information about the molecular column density and excitation conditions. For the two disks we have observed in two ¹²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 they are pursued in Section 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 Section 3.3.2. Finally, in Section 3.3.3, we use the constraints on the CS rotational temperature from Section 3.3.1 to estimate the CS column density in our entire disk sample.

3.3.1. Rotational Diagram Analysis

The first assumption for rotational diagram analysis is that the lines are in LTE. Lines can be assumed to be in LTE when their critical densities are surpassed by the gas densities in their emitting regions. The critical densities of the CS 5–4 and 6–5 transitions for temperatures of 20–50 K are $\sim [1.7\mbox{--}1.3]\times {10}^{6}$ and $\sim [2.9\mbox{--}2.2]\times {10}^{6}\,{\mathrm{cm}}^{-3}$, respectively (Shirley 2015). These are lower than expected in molecular disk layers (see, e.g., Figure 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. In the present study, the peaks of the CS brightness temperatures are found to be below ∼5 K. Thermalized lines usually have brightness temperatures that are ∼[0.6–0.8] × T_kin (unless optical depths are extreme), meaning that a 5 K brightness temperature line could partly be optically thick if T_kin ∼ 6–9 K. But, in particular for MWC 480, the lowest expected gas temperature in CS-emitting regions is ∼20–30 K (see Section 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 B. We then used the likelihood function ${ \mathcal L }({N}_{u},{T}_{\mathrm{ex}})$ described in Appendix B 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}_{\mathrm{ex}}$ and N_u describing the range of total column densities, N_tot, and excitation temperatures, ${T}_{\mathrm{ex}}$, consistent with the data. The following uninformative priors were assumed:

$$\begin{eqnarray}&&{T}_{\mathrm{ex}}({\rm{K}})={ \mathcal U }(3,300)\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{N}_{\mathrm{tot}}({\mathrm{cm}}^{-2})={ \mathcal U }({10}^{7},{10}^{20}).\end{eqnarray} \tag{ 2 }$$

Random draws and results from the posterior distributions are depicted in gray in Figure 5, along with the least-squares fit results, shown in blue. For the MWC 480 disk, the results converged toward {${T}_{\mathrm{ex}}\simeq 22.8\pm 2.5$ K, ${N}_{\mathrm{tot}}\simeq (7.3\pm 0.9)\times {10}^{12}\,{\mathrm{cm}}^{-2}$} for a flux integrated over a 150 au radius. For the LkCa 15 disk, the results converged toward {${T}_{\mathrm{ex}}\simeq 26.8\pm 2.3$ K, ${N}_{\mathrm{tot}}\,\simeq (1.2\pm 0.1)\times {10}^{13}\,{\mathrm{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.

Zoom In Zoom Out Reset image size

We next calculate the column density of CS as a function of radius, using the same rotational diagram analysis to the azimuthally deprojected radial intensities. The results are presented in Figure 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.

Zoom In Zoom Out Reset image size

3.3.2. Isotopologue Ratio Analysis

The analysis in Section 3.3.1 assumed that the main isotopologue CS lines are optically thin. This assumption can be explored for one disk, LkCa 15, where two CS isotopologue lines, ¹³CS and C³⁴S 6–5, are well detected.

The ratio of the two isotopologues ¹²CS and ¹³CS can be expressed as

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{{\rm{B}}}{(}^{13}\mathrm{CS})}{{T}_{{\rm{B}}}{(}^{12}\mathrm{CS})}=\displaystyle \frac{1-{e}^{-{\tau }_{{}^{13}\mathrm{CS}}}}{1-{e}^{-{\tau }_{{}^{12}\mathrm{CS}}}},\end{eqnarray} \tag{ 3 }$$

with T_B the brightness temperature (see Appendix C), 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 ¹³CS 6–5, and the latter is a reasonable approximation if the main isotopologue line does not become optically thick high up in the disk atmosphere. Furthermore, if ¹²CS and ¹³CS are both optically thin, then we have $(1-{e}^{-\tau })\to \tau$, which leads to

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{{\rm{B}}}{(}^{13}\mathrm{CS})}{{T}_{{\rm{B}}}{(}^{12}\mathrm{CS})}=\displaystyle \frac{{\tau }_{{}^{13}\mathrm{CS}}}{{\tau }_{{}^{12}\mathrm{CS}}}\,\approx R,\end{eqnarray} \tag{ 4 }$$

where R is the abundance ratio of ¹²CS to ¹³CS. Assuming that their abundance ratio is equal to the ¹²C/¹³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 disk-integrated brightness temperatures from the disk-integrated fluxes S_νΔv using the Rayleigh–Jeans law:

$$\begin{eqnarray}&&{T}_{{\rm{B}}}{\rm{\Delta }}v=\displaystyle \frac{{c}^{2}{S}_{\nu }{\rm{\Delta }}v}{2{k}_{{\rm{B}}}{\nu }^{2}{\rm{\Omega }}}.\end{eqnarray} \tag{ 5 }$$

To consistently calculate disk-integrated fluxes for the main and rare CS isotopologues, we applied a uv-taper to the extended baselines of the ¹²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

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{{\rm{B}}}{(}^{12}\mathrm{CS})}{{T}_{{\rm{B}}}{(}^{13}\mathrm{CS})}\ \simeq \ 17\ll {\left(\displaystyle \frac{{}^{12}{\rm{C}}}{{}^{13}{\rm{C}}}\right)}_{\mathrm{LISM}}.\end{eqnarray} \tag{ 6 }$$

If the ¹²CS/¹³CS abundance ratio is equal to the ¹²C/¹³C LISM ratio, Equation (6) implies that the ¹²CS 6–5 transition is optically thick. This implication can be used to calculate the ¹³CS column density and the ${}^{12}\mathrm{CS}$ excitation temperature. For an optically thick line $(1-{e}^{-{\tau }_{{}^{12}\mathrm{CS}}})\to 1$. Thus, Equation (33) directly provides T_ex(¹²CS) from the peak of its brightness temperature:

$$\begin{eqnarray}&&{T}_{\mathrm{ex}}{(}^{12}\mathrm{CS})=\displaystyle \frac{h\nu /{k}_{{\rm{B}}}}{\mathrm{ln}\left(1+\tfrac{h\nu /{k}_{{\rm{B}}}}{{T}_{{\rm{R}}}{(}^{12}\mathrm{CS})/f+{J}_{\nu }({T}_{\mathrm{bg}})}\right)}.\end{eqnarray} \tag{ 7 }$$

Equation (7) gives an excitation temperature of T_ex(¹²CS) ≃ 7 K, which, inserted into (33), gives an opacity of ${\tau }_{{}^{13}\mathrm{CS}}\simeq 0.06$ (this last step assumes that CS and ¹³CS have the same excitation temperatures). Moreover, since we assume that ¹³CS is optically thin, we have $(1-{e}^{-{\tau }_{{}^{13}\mathrm{CS}}})\to {\tau }_{{}^{13}\mathrm{CS}}$, which inserted into Equation (3) gives

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{{\rm{B}}}{(}^{13}\mathrm{CS})}{{T}_{{\rm{B}}}{(}^{12}\mathrm{CS})}\ \approx \ {\tau }_{{}^{13}\mathrm{CS}}.\end{eqnarray} \tag{ 8 }$$

Equation (8) gives an opacity consistent with the value obtained from Equations (7) and (33).

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

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\mathrm{tot}} & = & \displaystyle \frac{3h}{8{\pi }^{3}S{\mu }^{2}}\displaystyle \frac{{Q}_{\mathrm{rot}}({T}_{\mathrm{ex}})}{{g}_{u}}\displaystyle \frac{{e}^{{E}_{u}/{k}_{{\rm{B}}}{T}_{\mathrm{ex}}}}{{e}^{h\nu /{k}_{{\rm{B}}}{T}_{\mathrm{ex}}}-1}\\ & & \times \displaystyle \frac{1}{{J}_{\nu }({T}_{\mathrm{ex}})-{J}_{\nu }({T}_{\mathrm{bg}})}\int {T}_{R}{d}v\ {C}_{\tau }.\end{array}\end{eqnarray} \tag{ 9 }$$

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

Table 3.  Column Densities and Excitation Temperatures Derived from the Methods Presented in Section 3 when Applied

Species	Source	Radius	${N}_{\mathrm{tot}}$	${T}_{\mathrm{ex}}$	Method
		(au)	($\,{\mathrm{cm}}^{-2}$)	(K)
12CS	MWC 480	150	$(7.3\pm 0.8)\times {10}^{12}$	22.8 ± 2.5	Rotational diagram
	LkCa 15	200	$(1.2\pm 0.1)\times {10}^{13}$	26.8 ± 2.3	Rotational diagram
			$[1.5\times {10}^{13}\mbox{--}5\times {10}^{14}]$	∼7	Isotopologue ratio analysis
H2CS	MWC 480	300	${3.0}_{-0.4}^{+2.7}\times {10}^{12}$	${40.5}_{-20.2}^{+38.4}$	Rotational diagram
	LkCa 15	320	$\lesssim 2\times {10}^{12}$	40	Using Equations (26) and (27) and a fixed ${T}_{\mathrm{ex}}$

Download table as:  ASCIITypeset image

A similar computation can be done for the C³²S/C³⁴S ratio considering the elemental ³²S/³⁴S = 24.4 ± 5.0 in the vicinity of the Sun (Chin et al. 1996). Both observed ¹²CS(6–5)/C³⁴S(6–5) and ¹²CS(5–4)/C³⁴S(6–5) give a ratio of ∼13, which is somewhat lower than the cosmic ³²S/³⁴S ratio, but not nearly as much below the cosmic ratio as ¹²CS/¹³CS. Taken at face value, the excitation temperature and opacity are ${T}_{\mathrm{ex}}$(¹²CS) ≃ 7.4 K and ${\tau }_{{}^{34}\mathrm{CS}}=0.07$, leading to ${N}_{\mathrm{tot}}$(C³⁴S) $\simeq 6\times {10}^{11}\,{\mathrm{cm}}^{-2}$. Using the LISM ratio, we thus obtain ${N}_{\mathrm{tot}}$(¹²CS) $\simeq 1.5\,\times {10}^{13}\,{\mathrm{cm}}^{-2}$, which is an order of magnitude lower than the value derived from the ¹²CS/¹³CS ratio, and similar to the rotational diagram result.

The results obtained from the different methods presented in Sections 3.3.1 and 3.3.2 to derive CS column densities and excitation temperatures are summarized in Table 3. Both isotopologue analyses suggest that the ¹²CS line is optically thick, but the large difference in derived ¹²CS opacities and column densities, as well as the low CS excitation temperatures, casts 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, 2018; Andrews et al. 2018; Favre et al. 2019), this is not obviously true. A beam filling factor of 20%–40% would increase the excitation temperature to 20–30 K, in better agreement with our model expectations (see Section 4.3). Although Figure 3 shows apparent asymmetries for the ¹³CS and C³⁴S emissions, these are likely due to noise in the image that propagates to the zeroth-moment map. Hence, 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-S/N 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 the disk midplane, but LTE may not hold if CS is instead emitting from the disk atmosphere.

Third, and perhaps most importantly, ¹³CS and C³⁴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., Woods & Willacy 2009; Sakai et al. 2010, 2013; Ossenkopf et al. 2013; Smith et al. 2015; Yoshida et al. 2015; Yu et al. 2016; Huang et al. 2017). For instance, the ¹³CS-enhanced abundances in disk molecular layers could be explained by an excessive ¹³C-carrier, such as ¹³CO (while ¹²CO remains self-shielded). Rapid ion–molecule gas-phase reactions could thus transfer ¹³C from ¹³CO to hydrocarbons, and then to ¹³CS. The degree of fractionation that would be needed to explain the ¹³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³⁴S data within error bars, for model comparisons and calculations of column densities in disks without constraints on the excitation temperature.

3.3.3. CS Column Density Estimates for the Disk Sample

In 4/6 disks we observed a single CS line and therefore have no constraint on the CS excitation temperature. To estimate column densities in these sources, we use the constraints on the CS rotational temperature in the LkCa 15 and MWC 480 disks and apply them to the disk sample. In Section 3.3.1 we found disk-averaged rotational temperatures of ∼23 ± 3 K 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, as well as an uncertainty of ±5 K. We then use Equations (26) and (27) to calculate the column density. The results are presented in Figure 7. The CS disk-integrated column density varies by less than an order of magnitude across the sample of disks. The average column density is $\approx 7\times {10}^{12}\,{\mathrm{cm}}^{-2}$, and the standard deviation is $\approx 2\times {10}^{11}\,{\mathrm{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 Equation (27), leading, for instance, to an average column density of $\approx 5\times {10}^{13}\,{\mathrm{cm}}^{-2}$ with a standard deviation of $\approx \pm 5\times {10}^{12}\,{\mathrm{cm}}^{-2}$ for a disk temperature of 10 ± 2 K, i.e., a factor of 5 higher column density.

Zoom In Zoom Out Reset image size

The three H₂CS ortho lines (8₁₇–7₁₆, 9₁₉–8₁₈, and 9₁₈–8₁₇) are detected toward MWC 480 and tentatively detected toward LkCa 15 (2σ–3σ). The integrated flux density maps are shown in Figure 8 and the radial profiles in Figure 9. The emission toward MWC 480 shows no substructure, including no central depression. However, the resolution and S/N are too low to rule out a similar central dip as is seen in CS. As mentioned in Section 3.3.2 for the ¹³CS and C³⁴S emissions, the quality of the H₂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-S/N and higher-resolution data would help in disentangling potential differences in beam filling factors for the different transitions that could impact the analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{T}_{\mathrm{ex}}({\rm{K}})={ \mathcal U }(10,100),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{N}_{\mathrm{tot}}({\mathrm{cm}}^{-2})={ \mathcal U }({10}^{7},{10}^{20}).\end{eqnarray} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

Random draws and results from the posterior distributions are depicted in gray in Figure 10, along with the least-squares fit results, shown in blue. The results converged toward {${T}_{\mathrm{ex}}\simeq {41}_{-20}^{+38}$ K, ${N}_{\mathrm{tot}}\simeq {3.0}_{-0.4}^{+2.7}\times {10}^{12}\,{\mathrm{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 S/N and closely spaced energy levels. Still, it is noteworthy that the best-fit rotational temperature is close to the rotational temperature of CS, assuming that CS is optically thin.

We used Equations (26) and (27) to derive an upper limit of ${N}_{\mathrm{tot}}\lesssim 2\times {10}^{12}\,{\mathrm{cm}}^{-2}$ when considering the three H₂CS lines tentatively detected toward LkCa 15, a rotational temperature of $T\approx 40$ K, and an ortho-to-para ratio of 3. This corresponds to less than a factor of 2 lower than the value observed toward MWC 480.

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.

Table 4.  Physical Parameters Used for Our Disk Models

Parameters	MWC 480a	LkCa 15b
Stellar mass, M⋆ (M⊙)	1.8	1.0
Disk mass, Md (M⊙)	0.18	0.03
Characteristic radius, Rc (au)	100	100
Outer cutoff radius, Rout (au)	500	550
Density power-law index, γ	1.5	1.5
Midplane temperature at Rc,c T100 au (K)	30	15
Atmosphere temperature at Rc, T100 au (K)	48	20
Surface density at Rc	5.7	0.9
Temperature power-law index, q	0.5	0.5
Vertical temperature gradient index, β	2	2
UV flux of reference, ${f}_{\mathrm{UV},{{\rm{R}}}_{{\rm{c}}}}$ (in units from Draine 1978)	8500d	2550e

Notes.

^aGuilloteau et al. (2011). ^bPiétu et al. (2007). ^cThe midplane temperatures are estimated from Equation (16), the luminosities listed in Table 1, and a typical flaring angle φ = 0.05. ^dFrom Dutrey et al. (2011), originally computed from the Kurucz (1993) ATLAS9 of stellar spectra. ^eFrom Dutrey et al. (2011), originally scaled from Bergin et al. (2004).

Download table as:  ASCIITypeset image

4.1.1. Temperature Profile

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

$$\begin{eqnarray}T(z)=\left\{\begin{array}{ll}{T}_{\mathrm{mid}}+({T}_{\mathrm{atm}}-{T}_{\mathrm{mid}}){\left[\sin \left(\displaystyle \frac{\pi z}{2{z}_{q}}\right)\right]}^{2\delta } & {\rm{if}}\ z\lt {z}_{q}\\ {T}_{\mathrm{atm}} & {\rm{if}}\ z\geqslant {z}_{q},\end{array}\right.\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&{T}_{\mathrm{mid}}={T}_{\mathrm{mid},{R}_{c}}{\left(\displaystyle \frac{r}{{R}_{{\rm{c}}}}\right)}^{-q},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{atm}}={T}_{\mathrm{atm},{R}_{c}}{\left(\displaystyle \frac{r}{{R}_{{\rm{c}}}}\right)}^{-q}\end{eqnarray} \tag{ 14 }$$

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

$$\begin{eqnarray}&&H=\sqrt{\displaystyle \frac{{k}_{{\rm{B}}}\ {T}_{\mathrm{mid}}\ {r}^{3}}{\mu \ {m}_{{\rm{H}}}\ G\ {M}_{\star }}},\end{eqnarray} \tag{ 15 }$$

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

$$\begin{eqnarray}&&{T}_{\mathrm{mid}}(r)\approx {\left(\displaystyle \frac{\varphi {L}_{\star }}{8\pi {r}^{2}{\sigma }_{\mathrm{SB}}}\right)}^{1/4},\end{eqnarray} \tag{ 16 }$$

with L_⋆ the stellar luminosity, σ_SB the Stefan–Boltzmann constant, and φ the flaring angle (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 it 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 of panels of Figure 11.

Zoom In Zoom Out Reset image size

4.1.2. Density Profile

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

$$\begin{eqnarray}&&\displaystyle \frac{\delta \mathrm{ln}\rho (z)}{\delta z}=-\left[\displaystyle \frac{{{GM}}_{\star }z}{{r}^{3}}\displaystyle \frac{\mu \,{m}_{{\rm{H}}}}{{k}_{{\rm{B}}}T}+\displaystyle \frac{\delta \mathrm{ln}T}{\delta z}\right],\end{eqnarray} \tag{ 17 }$$

which gives for the isothermal case

$$\begin{eqnarray}&&\displaystyle \frac{\delta \rho (z)}{\rho (z)}=-\left[\displaystyle \frac{{{GM}}_{\star }z}{{r}^{3}}\displaystyle \frac{\mu \,{m}_{{\rm{H}}}}{{k}_{{\rm{B}}}{T}_{\mathrm{mid}}}\right]\delta z,\end{eqnarray} \tag{ 18 }$$

with solution

$$\begin{eqnarray}&&\rho (z)={\rho }_{o}\exp \left(-\displaystyle \frac{{z}^{2}}{2{H}^{2}}\right),\end{eqnarray} \tag{ 19 }$$

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:

$$\begin{eqnarray}&&{\rm{\Sigma }}={\rho }_{0}{\int }_{-\infty }^{+\infty }\exp \left(-\displaystyle \frac{{z}^{2}}{2{H}^{2}}\right)=\sqrt{2\pi }{\rho }_{0}H.\end{eqnarray} \tag{ 20 }$$

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 a sharp cutoff at specified inner and outer radii (i.e., R_int < r < R_out):

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{c}{\left(\displaystyle \frac{r}{{R}_{c}}\right)}^{-3/2},\end{eqnarray} \tag{ 21 }$$

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 $\dot{M}=-2\pi r{\rm{\Sigma }}(r){v}_{r}$, with v_r the radial velocity of the disk. For an axis-symmetric disk, mass conservation thus reads

$$\begin{eqnarray}&&2\pi r\displaystyle \frac{\delta {\rm{\Sigma }}}{\delta t}=-\displaystyle \frac{\dot{M}}{\delta r}.\end{eqnarray} \tag{ 22 }$$

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

$$\begin{eqnarray}&&{M}_{\mathrm{disk}}={\int }_{{R}_{\mathrm{int}}}^{{R}_{\mathrm{out}}}2\pi r{\rm{\Sigma }}(r){dr}.\end{eqnarray} \tag{ 23 }$$

Thus, the surface density at a radius R reads as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{R}=\displaystyle \frac{{M}_{\mathrm{disk}}{R}^{-3/2}}{4\pi \sqrt{{R}_{\mathrm{out}}}},\end{eqnarray} \tag{ 24 }$$

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 of panels of Figure 11, for each model.

4.1.3. Visual Extinction Profile

4.1.4. UV Flux Distribution

In the disk model the simplifying assumption of parameterizing the UV flux from the star as a multiple of the standard ISM radiation field in the units of Draine (1978) 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 spectra can be dominated by accretion luminosity and therefore Lyα, instead. Thus, for T Tauri stars, we are probably underestimating UV penetration into the disk, since Lyα 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 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}_{\mathrm{UV},{{\rm{R}}}_{{\rm{c}}}}$ (i.e., the UV flux factor at a radius R_c; see Table 1), is divided by two (e.g., Wakelam et al. 2016). This zeroth-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 as

$$\begin{eqnarray}&&{f}_{\mathrm{UV}}=\displaystyle \frac{{f}_{\mathrm{UV},{R}_{{\rm{c}}}}/2}{{\left(\tfrac{r}{{R}_{c}}\right)}^{2}+{\left(\tfrac{4H}{{R}_{c}}\right)}^{2}}.\end{eqnarray} \tag{ 25 }$$

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 of panels of Figure 11, for each model.

Once built, the 1+1D physical structure described above (and depicted in Figure 11 for the MWC 480 and LkCa 15 cases) is fed as an input in the pseudo-time-dependent astrochemical modeling code Nautilus v1.1 (Wakelam et al. 2016) used in three-phase mode (Ruaud et al. 2016). This rate-equation gas-grain code—based on Hasegawa et al. (1992) and Hasegawa & Herbst (1993)—simulates the time-dependent chemistry of ∼1100 species (half in the gas phase and half in solid phase) linked together via more than ∼12,000 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 nonthermal 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 photodesorption yield for all the molecules rather than individual ones, only experimentally determined for a few species. Thus, we set a constant generic photodesorption yield of $1\times {10}^{-4}$ molecule photon^–1 (Andersson & van Dishoeck 2008) for all the species contained in our chemical network. In the gas phase typical bi-molecular ion–neutral and neutral–neutral reactions are considered, as well as cosmic-ray-induced processes, photoionizations, and photodissociations caused by both stellar and interstellar UV photons. A dust-to-gas 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 (Fuente et al. 2017; Le Gal et al. 2017; Vidal et al. 2017) that are relevant for the purpose of this work.

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

Table 5.  Initial Elemental Abundances

Species	${n}_{i}/{n}_{{\rm{H}}}$	Reference
H2	0.5
He	9.0 × 10−2	1
C+	1.7 × 10−4	2
N	6.2 × 10−5	2
O	2.4 × 10−4	3
S+	8.0 × 10−8	4
Si+	8.0 × 10−9	4
Fe+	3.0 × 10−9	4
Na+	2.0 × 10−9	4
Mg+	7.0 × 10−9	4
P+	2.0 × 10−10	4
Cl+	1.0 × 10−9	4
F+	6.7 × 10−9	5

References. (1) Wakelam & Herbst (2008); (2) Jenkins (2009); (3) Hincelin et al. (2011); (4) Graedel et al. (1982); (5) Neufeld et al. (2015).

Download table as:  ASCIITypeset image

We ran the disk chemistry models for 10⁶ yr. 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 as young or twice as old. 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 simulating the LkCa 15 disk (see Figure 11). Both models start from the same initial chemical composition, which results from the zero-dimensional model described above and considers a "low" sulfur elemental abundance of $8\times {10}^{-8}$ (see Table 5; Graedel et al. 1982; Wakelam & Herbst 2008), i.e., ∼2 orders of magnitude below the cosmic value $\sim 1.5\times {10}^{-5}$ (Asplund et al. 2009). This "low" sulfur elemental abundance corresponds to the typical sulfur-depleted value used in models to reproduce the abundances observed in dense astrophysical molecular environments (e.g., Tieftrunk et al. 1994; Wakelam et al. 2004; Vidal et al. 2017; Vastel et al. 2018). For the present study, this approximation seems reasonable since CS and H₂CS are supposed to reside in disk molecular layers, which are similar to dense astrophysical molecular environments.

The resulting number densities of CS, H₂CS, and their main precursors, S⁺ and S, at 1 Myr are presented in Figure 12 for the MWC 480 and LkCa 15 disk models (the corresponding fractional abundances are shown in Appendix D). Figure 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 Figure 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 Figure 12 with Figure 11, we see that the layers where CS and H₂CS reside have similar densities of $\sim {10}^{6}-{10}^{8}\,{\mathrm{cm}}^{-3}$.

Zoom In Zoom Out Reset image size

In both disks the S and H₂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 photoionization 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₂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_yH (with x = 1–4 and y = 2–3). This produces carbonated S-ions, including HCS⁺, CS⁺, HC₃S⁺, and C₂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⁺, ${{\rm{H}}}_{3}^{+}$, and HCO⁺) and with He⁺, and freeze-out onto grain surfaces.

H₂CS is mainly produced by the gas-phase neutral–neutral reaction of atomic S and CH₃ and the electronic dissociative recombination of H₃CS⁺, whose formation begins with the reaction between S⁺ and CH₄. 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₂CS releasing some (∼1%) of its H₂CS product in the gas phase by chemical reactive desorption. Similarly to CS, H₂CS is mainly destroyed by photodissociation processes in the disk atmosphere, reactions with protons and protonated ions (i.e., with H⁺, ${{\rm{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 Figure 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 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 ¹³CS observations toward LkCa 15 being an order of magnitude higher is thus inconsistent with our model predictions. However, we note that our disk models were run with an initially depleted S-element abundance (see Section 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₂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.

Zoom In Zoom Out Reset image size

The predicted H₂CS column densities are 1–2 orders of magnitude below the predicted CS column densities in each disk. Similarly to CS, the H₂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₂CS precursors, with the possible exception of S toward LkCa 15. Compared with observations toward MWC 480, the model underpredicts H₂CS by several orders of magnitude. Furthermore, the model predicts that more H₂CS should be present toward LkCa 15 compared to MWC 480, which is opposite to what we observed. Finally, the model predicts that H₂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₂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 of ∼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 range of ∼25–200 K and ∼60–400 K, respectively. Similar behaviors can be seen for H₂CS. While the H₂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 Figure 12). The rotational temperatures of CS and H₂CS are also overplotted on Figure 14. These temperatures were derived from observations out to ∼190 au (MWC 480) and ∼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.

Zoom In Zoom Out Reset image size

In Section 4.3 we compared model predictions of CS and H₂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¹²–10¹³$\,{\mathrm{cm}}^{-2}$ are well reproduced, within a factor of a few, by our models. The best fits are found for radii in the range of ∼150–250 au for MWC 480 and for radii of ∼100 au and in the range of ∼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₂CS in the outer disk of MWC 480 and incorrectly predict that H₂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 15 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₂CS is that grain-surface formation of H₂CS and subsequent photodesorption may be more important in disks than is currently implied by the model. The surface and bulk ice H₂CS 2D disk densities are represented in Figure 15, as well as those of CS, for both MWC 480 and LkCa 15 disk models. We see that both molecules are abundant in ices, especially in the disk midplanes. However, photo-irradiation laboratory experiments of CH₃OH ices showed an efficient dissociation of the product rather than a desorption (Bertin et al. 2016; Cruz-Diaz et al. 2016). H₂CS being also a rather large molecule, it could similarly efficiently be photodissociated in the ice. Our present model already considers ice photodissociation rates; however, due to the lack of laboratory studies on that topic, those are set to be equal to the gas-phase photodissociation ones. Hence, our model might mispredict the H₂CS ice abundance and desorption, unless additional ice formation pathways are missing in our astrochemical networks. Therefore, both photo-irradiation of icy H₂CS and new H₂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₂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₂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₂CS observations and its predictions from the astrochemical disk modeling. Thus, higher H₂CS angular resolution observations would be required toward the same disk sample to further investigate these issues.

Zoom In Zoom Out Reset image size

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 (Figure 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 needs 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 straightforward 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 Figure 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 Section 4.3, the MWC 480 disk being more massive, its visual extinction is larger (see Figure 11), implying less UV penetration in the disk. In addition, the dramatic changes in temperature over small radii shown in Figure 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 of the MWC 480 observations presented here versus the T Tauri ones, there might be some sulfur chemistry distinctions in between T Tauri and Herbig Ae disks.

Also, we should note that the impact of X-rays on disk chemistry is not included in our present model. X-rays could affect the the ion abundances (e.g., Glassgold et al. 1997; Aikawa & Herbst 1999; Henning et al. 2010; Cleeves et al. 2015; Teague et al. 2015) such as the abundances of S⁺, CS⁺, HCS⁺, as well as the SO/SO₂ formation, and thus affect 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-ray chemistry is not dominant (e.g., Rab et al. 2018). Nevertheless, this would be an interesting future modeling improvement.

Compared to other disk studies, the CS column densities derived toward our disk sample (see Figure 7) are consistent with the values found in other disks, e.g., $(2.0\pm 0.16)\times {10}^{12}\,{\mathrm{cm}}^{-2}$ in GO Tau (Dutrey et al. 2011), $(8.7\pm 1.6)\times {10}^{12}\,{\mathrm{cm}}^{-2}$ in LkCa 15 (Dutrey et al. 2011), $(6\pm 3)\times {10}^{12}\,{\mathrm{cm}}^{-2}$ in DM Tau (Semenov et al. 2018), and $2.2\times {10}^{13}\,{\mathrm{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₂CS disk observations in a more holistic 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 PDRs (Rivere-Marichalar et al. 2019), molecular clouds (Drozdovskaya et al. 2018), and cold dense core (Vastel et al. 2018), analogs of the main three chemical layers constituting PPDs, 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 PPDs (from this work and previous studies) versus CS, compared to their values in the different ISM environments aforementioned.

Zoom In Zoom Out Reset image size

The ¹³CS/¹²CS value we derived in disks from this study is found to be higher (∼0.06) than those measured so far in the different ISM environments presented in Figure 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 ¹³C enhancement suggested in this work. Higher-resolution observations of CS isotopologues in disks are needed to further investigate this ¹³C enhancement hint.

The C³⁴S/C³²S ratio we observed in disks (∼0.08) is found, within error bars, to be a factor of 2 higher than values found in PDRs and dense cores and than the LISM value. The observed C³⁴S/C³²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 Sun-like 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 H₂CS/CS ratio decreasing from dense to diffuse ISM is found in disks (i.e., in LkCa 15) to be 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₂CS/CS ratio behaves in disks 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 H₂CS/CS ratio, the H₂S/CS is increasing from dense to diffuse ISM and the disk value (Phuong et al. 2018) is found to be almost as low as the lower limit derived in dense cores (Vastel et al. 2018). This could also explain why H₂S was so hard to detect in disks (e.g., Dutrey et al. 2011), aside from the fact that the H₂S detection could have been facilitated by the high mass of the GG Tau disk in which it was detected. Also, even though high uncertainties remain on values observed toward IRAS 16293–2422, H₂S appears as one of the main abundant S-species 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, which would also explain the boost in S-chemistry observed there (Drozdovskaya et al. 2018). However, and in opposition to the H₂S/CS ratio, the OCS/CS ratio is found to be lower in PDRs than in dense cores, suggesting that the OCS and H₂S molecules trace different types of environments.

As with the H₂S/CS ratio, the SO/CS ratio increases from dense to diffuse ISM. However, only a lower limit on SO is provided toward IRAS 16293–2422 (Drozdovskaya et al. 2018). Semenov et al. (2018) provided an 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 cores (Vastel et al. 2018). Contrary to the SO/CS ratio, the SO₂/CS ratio appears rather flat across the different environments compared here and so does not seem to characterize any specific S-chemistry behavior that would depend on its physical environment. The PDR value is found to be compatible with the 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) does not allow us to settle any specific trend for this ratio.

Finally, the C₂S/CS ratio decreases from dense to diffuse ISM, in an even more pronounced way than the OCS/CS and H₂CS/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 (Rivere-Marichalar et al. 2019) thus seems to be characteristic of molecular cloud environments.

On average, the S-bearing molecules observed in PPDs would behave more like in cold cores or molecular clouds. In particular, Figure 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.