© The Authors 2025

1 Introduction

After the main-sequence and the first red-giant phase of their evolution, low- to intermediate-mass stars ($\sim 0.8M_{\odot }<M_{}<8M_{\odot }$) become cool, luminous giants ($T_{\mathrm{e}\mathrm{f}\mathrm{f}}\le 3000K,L_{\ast }\sim {10}^3--{10}^4{L}_{\odot },R_{\ast }\sim 100R_{\odot }$; Herwig 2005) known as asymptotic giant branch (AGB) stars. These stars are characterised by dust-driven stellar winds (Habing & Olofsson 2004; Höfner & Olofsson 2018) that cause extensive mass loss (10⁻⁸ to 10⁻⁵ M_ʘ yr⁻¹; Olofsson 1999; Höfner & Olofsson 2018), forming expanding circumstellar envelopes (CSEs) of dust and gas around them. The CSEs harbour a large number of molecular species (e.g. Cernicharo et al. 2000; Patel et al. 2011; Unnikrishnan et al. 2024; Decin 2021) and are also sites of dust formation (Habing & Olofsson 2004). AGB stars contribute significantly to the chemical enrichment of the interstellar medium (ISM) and thereby to overall galactic chemical evolution (e.g. Tielens 2005; Matsuura et al. 2009; Kobayashi et al. 2011).

AGB stars with a photospheric carbon-to-oxygen abundance ratio (C/O) of greater than 1 are termed C-type (carbon-rich) stars, whereas the ones with C/O < 1 are called M-type (oxygenrich). The CSEs of C-type AGB stars exhibit richer chemistry compared to M-type stars, featuring a wide variety of molecular species, including long carbon chains (Olofsson et al. 1993; Cernicharo et al. 2000; Gong et al. 2015; Woods et al. 2003; Unnikrishnan et al. 2024) and likely also polycyclic aromatic hydrocarbons (PAHs, e.g. Cherchneff 2012; Tielens 2008; Zeichner et al. 2023; Anand et al. 2023). More than 100 molecular species have been detected so far in carbon star CSEs (McGuire 2022; Cernicharo et al. 2023a,b; Tuo et al. 2024).

Some of these species, such as CO, C2H2, HCN, CS, SiO, and SiS, are of photospheric origin, influenced by nonequilibrium shock chemistry in the extended atmosphere (e.g. Cherchneff 2006), and are injected into the circumstellar gas from the warm, dense inner regions. They are commonly referred to as parent molecules (e.g. Woods et al. 2003; Agúndez et al. 2020), in contrast to the daughter species (e.g. CN, HNC, HC3N, see Agúndez et al. 2020) that are formed farther out in the CSE, mostly by photodissociation of the parent species. The complex interplay of various processes including gas-phase reactions, dust-gas interactions, and photodissociation determines the radial variation in molecular abundances across AGB CSEs. Therefore, observations of molecular species and dust at multiple wavelengths are needed to obtain a comprehensive picture of the physical conditions and chemical processes influencing circumstellar molecular abundances.

Most observational studies of circumstellar carbon chemistry so far, the interferometric ones in particular, have almost exclusively focussed on the nearby (~190 pc), high mass-loss rate (1.5×10⁻⁵ M_ʘ yr⁻¹, MLR), C-type AGB star IRC +10 216, which is generally regarded as an archetype carbon star and has reasonably well-determined molecular abundances (e.g. Agúndez et al. 2012, 2017; Cernicharo et al. 2000; Agúndez et al. 2015; Pardo et al. 2022a; Siebert et al. 2022; Velilla Prieto et al. 2015a). Abundances derived from single-dish (SD) line observations are, however, available for a few other carbon stars as well (e.g. Woods et al. 2003; Schöier et al. 2006, 2007, 2013; Danilovich et al. 2018; Massalkhi et al. 2024). In contrast to SD observations, spatially resolved interferometric maps of circumstellar line emission (e.g. Agúndez et al. 2017; Danilovich et al. 2019; Velilla-Prieto et al. 2019) offer the possibility to directly determine the spatial distribution and radial extent of the molecular gas. The first spatially resolved spectral survey of a set of carbon stars other than IRC +10 216 was presented by Unnikrishnan et al. (2024, hereafter UK24), aimed at observationally testing the archetypal status often attributed to IRC +10 216. This study revealed similar complexity in both morphology and chemistry in the observed sources as was previously found towards IRC +10 216 (e.g. Agúndez et al. 2017).

UK24 derived molecular abundances assuming local thermodynamic equilibrium (LTE) conditions in the circumstellar gas. Though such LTE calculations are regularly used to obtain reasonable first-order abundance estimates (e.g. Smith et al. 2015; Pardo et al. 2022a; Velilla Prieto et al. 2015b), non-LTE radiative transfer (RT) modelling is required to calculate precise radial abundance profiles (e.g. Agúndez et al. 2012; Brunner et al. 2018; Van de Sande et al. 2018a; Danilovich et al. 2018, 2019; Massalkhi et al. 2024). Most of the existing RT models of molecular line emission from carbon star CSEs are limited to either sampling a range of excitation conditions using SD spectra (e.g. Woods et al. 2003; Schöier et al. 2007, 2013; Agúndez et al. 2012; Danilovich et al. 2018; Massalkhi et al. 2019, 2024) or focussing on the spatial distribution of the emission using interferometric observations of a few selected lines (e.g. Schöier et al. 2006; Agúndez et al. 2017; Velilla-Prieto et al. 2019). RT modelling with strict constraints based on the combination of the above, i.e. both excitation information from a broad range of SD data, and morphological information from spatially resolved images, as presented in this work, is the necessary next step in advancing the modelling of circumstellar line emission to estimate precise molecular abundances, as has been noted by several authors including Danilovich et al. (2018) and Massalkhi et al. (2024). This has been done for both M-type stars (e.g. Danilovich et al. 2019) and S-type stars (e.g. Brunner et al. 2018), yielding good estimates of fractional abundances of several molecular species including CS.

Among the parent species, CS is easily formed in carbonstar envelopes under thermodynamical equilibrium conditions due to the availability of large amounts of C (e.g. Agúndez et al. 2020; Danilovich et al. 2018), while SiS and SiO formation may require shock-induced chemistry (e.g. Cherchneff 2012). Being refractory species, they are expected to play a major role in circumstellar dust formation and growth as well (see Gail et al. 2013, and references therein). Factors like these affect the distribution of such species in the CSE, and can complicate their RT modelling, making it difficult to obtain simultaneous good fits to a large number of observed lines (e.g. Velilla-Prieto et al. 2019; Schöier et al. 2007). Therefore, the modelling of such species, particularly for the high-MLR C-type stars, demands a more detailed analysis than in the case of CS, including the use of non-Gaussian abundance profiles, especially when using both spatial and excitation constraints together. In the current work, we focus on CS and its isotopologues, and will present the modelling of the remaining parent species in a forthcoming paper. We present, for the first time, CS RT models constrained using spatially resolved observations for carbon stars other than IRC +10 216.

This work is the second in a series of planned publications, starting with UK24, that aim to leverage the high-angularresolution of ALMA to expand spatially resolved spectral line observations to other carbon stars beyond the well-studied IRC +10 216. The paper is organised as follows. The source sample is presented in Sect. 2. The observations used in this work are described in Sect. 3. Sect. 4 explains the RT modelling method employed, with Sect. 4.1 describing the adopted envelope model, and the dust, CO, and CS modelling procedures detailed in Sections 4.2, 4.3, and 4.4, respectively. We present the results in Sect. 5, and compare and discuss them in light of existing literature and chemical models in Sect. 6.

2 The sources

Our carbon star sample includes IRAS 15194-5115 (II Lup), IRAS 15082-4808 (V358 Lup), IRAS 07454-7112 (AI Vol), AFGL 3068 (LL Peg), and IRC +10 216 (CW Leo). The reasoning behind the source selection and the properties of the Mira variables IRAS 15194-5115, IRAS 15082-4808, and IRAS 07454-7112 were already discussed by UK24. In the current work, we also include the extreme carbon star AFGL 3068 along with the archetypical carbon star IRC +10 216, to enable direct comparisons with the rest of the sample. The basic parameters of the five stars are summarised in Table 1, and the two additions to the original sample are briefly described below.

AFGL 3068 is a Mira variable with a pulsation period of 696 days (Whitelock et al. 2006). It is an extreme carbon star, expected to be at the tip of the AGB phase. It has a well-defined large-scale spiral pattern in its CSE, thought to be caused by an embedded binary companion (e.g. Kim et al. 2017; Guerrero et al. 2020). Similar arc and shell patterns have been observed towards IRC +10 216 (Agúndez et al. 2017; Mauron & Huggins 1999; Leão et al. 2006; Dinh-V-Trung & Lim 2008), and also in molecular emission around IRAS 15194-5115 (Lykou et al. 2018, UK24).

IRC +10216 is also a Mira variable with a period of ~630-678 days (Menten et al. 2012; Groenewegen et al. 2012; Kim et al. 2015; Pardo et al. 2018; Dharmawardena et al. 2019). It is the brightest astronomical source at 5 μm excluding solar system objects (Becklin et al. 1969), and is one of the most observed evolved stars at mm/sub-mm/radio wavelengths (see e.g. Tuo et al. 2024, and references therein). All five sources in our sample are expected to be in their high mass-loss, late AGB phase (Vassiliadis & Wood 1993).

3 Observations

3.1 Photometric fluxes

To characterise the stellar parameters and circumstellar dust properties, we need to model the stellar and dust emission contributing to the spectral energy distributions (SEDs) of our sources. We used photometric flux densities compiled from the literature to constrain our dust RT models. High-quality fluxes were obtained from Gaia DR3 (Gaia Collaboration 2023) in the B_p, G, and R_p bands, where available, and from the 2MASS J, H, and K_s bands (Cutri et al. 2003) with quality flags ranging from 1 to 3. Additional fluxes were included from the DIRBE L and M bands (Price et al. 2010; Smith et al. 2004), the IRAS Point Source Catalogue at 12, 25, 60, and 100 μm (Helou & Walker 1988) with quality flag 3, and from Akari at 8.6-160 μm (Ishihara et al. 2010; Yamamura et al. 2010), also with quality flag 3, where available. Data from the Planck catalogue (Planck Collaboration VII 2011) and the SCUBA Legacy Catalogues (Di Francesco et al. 2008) were also included when available. WISE data were excluded due to well-documented saturation issues with AGB stars, which led to significant flux uncertainties (e.g. Lian et al. 2014). We note that the WISE bands W1 to W4 anyway overlap in wavelength with the DIRBE L and M bands, and IRAS at 12 and 25 μm.

The uncertainty in flux measurements can reflect the variability of the star, which is more pronounced at shorter wavelengths and diminishes progressively at longer wavelengths (e.g. Pardo et al. 2018; Ramstedt et al. 2008). Based on this, we assigned uncertainties of 75% for the optical measurements from Gaia, 60%, 50%, and 40% for the 2MASS J, H, and Ks band fluxes, respectively, 30% for the L and M band fluxes, and 20% for fluxes beyond 5 μm, in line with the estimates of Ramstedt et al. (2008). We corrected the fluxes below 3.5 μm for interstellar extinction following the method described in Andriantsaralaza et al. (2022). The photometric flux densities used in this work for all sources are listed in Table A.1.

3.2 CO line emission

We used archival SD CO observations, along with Atacama Compact Array (ACA) interferometric observations of CO J = 2-1 and J = 3-2 lines from the DEATHSTAR project (Ramstedt et al. 2020; Andriantsaralaza et al. 2021) for IRAS 07454-7112 and AFGL 3068, to constrain our gas RT models. We refer to Ramstedt et al. (2020) and Andriantsaralaza et al. (2021) for details of the DEATHSTAR ACA observations and data processing. The majority of the archival SD CO data used in this work were sourced from Olofsson et al. (1993), Schöier & Olofsson (2001), Woods et al. (2003), Ramstedt et al. (2008), De Beck et al. (2010), Danilovich et al. (2015), the APEX pointing cata-logue¹, and the James Clerk Maxwell Telescope (JCMT) cata-logue². While most of these observations are of the CO J = 1-0, J = 2-1, and J = 3-2 transitions, IRAS 15194-5115 and IRAS 07454-7112 also have several higher J transitions up to J = 1615 and J = 14-13, respectively, observed using Herschel/HIFI as part of the SUCCESS (Danilovich et al. 2015) project. The archival SD CO line intensities used in this work are listed in Table B.1.

Multiple spectra were extracted from progressively increasing apertures from the ACA cubes by convolving them to larger and larger circular beam sizes, starting at the major-axis size of the original synthesised beam. The beam size was gradually increased in steps of half the initial circular beam until the measured line flux plateaued, indicating that all detected emission is retrieved. This spectral extraction approach was employed to account for flux contributions across various spatial scales, and captures the variation in line intensities with distance from the centre, which serves as an additional constraint to the RT models. The integrated intensities of the CO spectra extracted from various apertures are listed in Table B.2.

The measurement errors on the CO line intensities were dominated by telescope calibration uncertainties, which we assumed to be 20% of the line intensity for SD data (see e.g. Schöier et al. 2007; Danilovich et al. 2018). In cases with low signal-to-noise ratio, like the archival SEST J = 3-2 data, the overall uncertainty on the intensity was estimated to be around ~30%, taking into account the high noise in the spectra as well. For the ALMA lines, we set the uncertainties to 10%, following Ramstedt et al. (2020), Andriantsaralaza et al. (2021), and UK24.

3.3 CS line emission

We have modelled the CS, ¹³CS, and C³⁴S lines from our ALMA, APEX, and Herschel/HIFI spectral surveys (see Sections 3.3.1, 3.3.2, 3.3.3) in this work. We could not model the C³³ S emission due to the lack of line detections having good signal-to-noise ratios. In some cases, we also used supplementary spectral data obtained from the literature (see Sect. 3.3.4). Our ALMA line survey of three C-type AGB stars (Sect. 3.3.1) has already been presented in detail in UK24, along with selected APEX molecular line observations that are also part of an unbiased SD spectral survey. This APEX survey and the HIFI spectral scans are discussed below, along with the supplementary archival observations used. The integrated intensities of the CS lines used in this work are listed in Table 2, and those of the ¹³CS and C³⁴S lines are given in Tables D.1 and D.2, respectively.

3.3.1 ALMA line survey

ALMA band 3 (85-116 GHz) spectral surveys (project code: 2013.1.00070.S; PI: Nyman, L. Å.) were presented by UK24 for IRAS 15194-5115, IRAS 15082-4808, and IRAS 07454-7112, and include the J = 2-1 line of CS and its isotopologues ¹³CS and C³⁴S. AFGL 3068 was also originally included in our ALMA survey, but the observations for this source were not completed, and yielded no relevant lines, leading us to rely only on SD observations (see Sections 3.3.2, 3.3.4) for this source. We extracted spectra from multiple, progressively increasing apertures for the spatially resolved ALMA CS lines, using the same procedure as presented above for the ACA CO lines (Sect. 3.2). We see no evidence of resolved-out flux in our CS J = 2-1 maps. The multi-aperture line profiles are used along with the SD spectra to constrain our RT models (see Sect. 4.4.2), ensuring that both spatial and excitation constraints are in place for the determination of the best-fit models. Following UK24, we use a standard 10% calibration uncertainty on the CS integrated intensities for lines observed with ALMA. For more details of the ALMA survey, including the specifics of the observations, data processing, and availability of processed data, we refer the reader to UK24.

3.3.2 APEX line survey

All five sources in our sample were observed using the Atacama Pathfinder Experiment (APEX, Güsten et al. 2006), in

the frequency ranges 159-211 GHz (band 5) and 200-270 GHz (band 6). Additionally, IRAS 15194-5115 was also observed at the frequencies 272-376 GHz (band 7). We used only the CS (and isotopologues) lines (Table 2) in the present work. Band 5 observations were conducted using the SEPIA receiver (Billade et al. 2012; Belitsky et al. 2018) for all sources. IRAS 15194-5115 observations in bands 6 and 7 were carried out using the SHeFI receiver (Vassilev et al. 2008), while the PI230³ instrument was used for band 6 observations of the rest of the sources. The SHeFI bands 6 and 7 observations of IRAS 15194-5115 were carried out in October 2014. The SEPIA Band 5 observations were performed in October 2016 for IRAS 15194-5115, and in June-July 2021 for the other four sources. The PI230 observations were done in July and October 2019.

All observations were performed in wobbler switching mode, with a standard wobbler throw of 1′. The SEPIA band 5 and PI230 observations were sideband separated (2SB). Singlesideband (SSB) observations with SHeFI for IRAS 15194-5115 were, when possible, optimised to avoid image sideband leakage from bright emission lines. Pointing, focus, and calibration were regularly checked with the standard catalogues employed at APEX, and we note that the science targets themselves are part of these catalogues.

The APEX spectra were originally delivered in the antenna temperature ($T_{\mathrm{A}}^{\ast }$, Ulich & Haas 1976) scale, which has already been corrected for atmospheric attenuation (see e.g. Pardo et al. 2022b, 2025, and references therein). We converted all spectra used in this work to the main beam temperature (T_MB) scale, using $T_{\mathrm{M}\mathrm{B}}=T_{\mathrm{A}}^{\ast }/{\eta }_{\mathrm{M}\mathrm{B}}$, where η_MB is the main beam efficiency. The adopted values for APEX main beam efficiency³ are 0.8 for SEPIA band 5, 0.75 for PI230 and SHeFI band 6, and 0.73 for SHeFI band 7. The APEX beam⁴ varies from ~39″ at 159 GHz to ~16″ at 376 GHz. Data reduction was performed using the GILDAS/CLASS⁵ software. A first-order polynomial baseline, calculated on the line-free regions, was subtracted from the spectra. As in UK24, we assume a standard 20% calibration uncertainty for the APEX data.

3.3.3 HIFI survey of II Lup

We present observations of CS, ¹³CS, and C³⁴S line emission towards IRAS 15194-5115 from the high-sensitivity, highresolution spectral scan (PI: De Beck, E.) carried out between 5 and 19 September 2012 using the Heterodyne Instrument for the Far Infrared (HIFI, de Graauw et al. 2010) instrument on board the Herschel Space Observatory (Pilbratt et al. 2010). The observations cover the frequency range 479.5-1121.9 GHz, spanning bands 1a, 1b, 2a, 2b, 3a, 3b, 4a, and 4b. All observations were carried out in spectral scan mode, using a dual beam switch calibration scheme with off-source positions 3′ away from the science target’s position. To improve baseline stability in bands 3 and 4, the fast chop mode was selected; this was not needed for bands 1 and 2. As in the case of the APEX data, we used a first-order polynomial for baseline subtraction.

For a detailed description of the different observing modes, the in-orbit performance of HIFI, and the data processing schemes, we refer to Roelfsema et al. (2012), Mueller et al. (2014), and Shipman et al. (2017). The conversion of the HIFI spectral scans from double-sideband spectra to single-sideband spectra follows the conjugate gradient method described by Comito & Schilke (2002) taking into consideration all relevant sideband gain ratios. The data used here were extracted from the highly processed data product (HPDP⁶,⁷) made available in the Herschel Science Archive as part of an effort to uniformly process 500 spectral scan observations towards different objects. The HPDP consists of one spectrum per polarisation (horizontal and vertical) for each band (1a through 4b). For our purposes, we averaged the spectra obtained for the horizontal and vertical polarisations.

The nominal spectral resolution of HIFI is 1.1 MHz, corresponding to a velocity resolution of 0.3-0.7 km s⁻¹ across the survey. Beam widths vary from 44″.2 at 479.5 GHz to 18″.9 at 1121.9 GHz (Roelfsema et al. 2012, see their Eq. (3)). The absolute pointing accuracy is ~0″.8, with a pointing stability of 0″.3⁸. HIFI main-beam efficiencies (~0.54-0.66 in our frequency range) for $T_{\mathrm{A}}^{\ast }$ to T_MB conversion, and calibration uncertainty estimates (~10%) are adopted from Mueller et al. (2014). The overall uncertainty on the line intensities used in this work can be larger, especially for the higher J lines due to relatively low signal-to-noise ratios, as can be seen from Fig. 1. The full molecular inventory from the APEX and HIFI surveys will be presented in a future work.

3.3.4 Supplementary CS observations

In addition to the ALMA, APEX, and HIFI data described above, we used archival CS lines available from the literature to provide better constraints to our RT models. For IRAS 07454-7112, we used APEX spectra of CS J = 4-3 and 6-5 lines from Danilovich et al. (2018). We also obtained the same lines for IRAS 15194-5115, for which we already had these lines from our APEX survey. In this case, we compared the line profiles from both APEX datasets and found them to be identical within calibration uncertainties.

For AFGL 3068, we used an ACA CS J = 7-6 data cube from the DEATHSTAR project (Ramstedt et al. 2020; Andriantsaralaza et al. 2021), and an OSO 20 m CS J = 2-1 integrated intensity reported by Woods et al. (2003). For IRC +10 216, we obtained an interferometric data cube of the CS J = 2-1 line, produced by combining ALMA data with an IRAM 30 m on-the-fly map, from Velilla-Prieto et al. (2019). This helped to recover all the emission for this very extended source (see e.g. Cernicharo et al. 2015), without being restricted by the maximum recoverable scale of the ALMA observations. We also received a Yebes 40 m spectrum of the CS J = 1-0 line, along with IRAM 30 m spectra of the J = 2-1, 3-2, 4-3, 5-4, 6-5, and 7-6 lines from Massalkhi et al. (2024). From these, we used all lines except the CS J = 7-6, for which Massalkhi et al. (2024) mention potential observational issues, in our RT modelling. The calibration uncertainty for the IRAM 30 m varies from ~10%-30% depending on the observation frequency of observation (see Massalkhi et al. 2024). The integrated intensities of the supplementary CS lines used in this work are listed in Table 2.

Fig. 1 Observed (black) and modelled (red) CS line profiles for IRAS 15194-5115. The transition quantum numbers, telescope used for the observation, and the corresponding beam size in arcseconds are listed at the top right corner of each panel. The line fluxes are in units of Jansky for the ALMA spectra, whereas they are given in the main-beam temperature (TMB [K]) scale for all SD spectra shown.

4 Radiative transfer modelling

4.1 Envelope model

We adopted a spherically symmetric, smooth CSE model, assuming an isotropic, constant mass-loss rate (MLR, Ṁ) and a constant dust-to-gas mass ratio (Ψ). Though some of our sources display asymmetric features in their CSEs (see Sect. 2), we were still able to reproduce the overall line profiles and radial intensity profiles from the ALMA observations with these assumptions, similar to those regularly used in the literature to model similar sources (e.g. Massalkhi et al. 2024; Danilovich et al. 2015, 2018, 2019; Velilla-Prieto et al. 2019; Schöier et al. 2007, 2013; Agúndez et al. 2012, 2017). A full treatment of the circumstellar asymmetries and density variations is beyond the scope of the modelling presented in this paper.

The central AGB star is characterised by its luminosity (L_⋆) and temperature (T_⋆). Dust is assumed to be present only beyond the dust condensation radius (R_c ), which is assumed to be at three times the stellar radius (e.g. Höfner 2009; Bladh & Höfner 2012; Höfner & Olofsson 2018) in our models (see Sect. 4.2). We assume a smoothly accelerating wind with a terminal radial expansion velocity ν_∞ (Table 1). Between the surface of the star and the dust condensation radius, the gas expansion velocity (ν_exp) is assumed to be a constant (ν₀ = 3 km s⁻¹), approximately equal to the sound speed at R_c (e.g. Danilovich et al. 2015). As the gas volume in this region (1-3 R_⋆) is very low compared to the rest of the CSE, any variation in the velocity profile here will not significantly affect our results. Beyond R_c , the velocity profile can be expressed as (Höfner & Olofsson 2018) $${\nu }_{\text{exp}}(r)={\nu }_0+({\nu }_{\mathrm{\infty }}-{\nu }_0){(1-\frac{R_{\text{c}}}{r})}^{\beta },$$(1)

where the parameter β controls the wind acceleration. We also adopted a constant turbulent velocity of 1.5 km s⁻¹ and assumed that dust and gas have the same velocity. The circumstellar gas density was assumed to decrease as r⁻², given by $$n_{{\text{H}}_2}=\frac{{\dot{M}}_g}{4\pi r^2\nu (r)}\frac{1}{m_{{\text{H}}_2}},$$(2)

where Ṁ_g is the gas MLR, ν(r) the wind velocity at radius r, and m_H2 is the molecular mass of hydrogen. The circumstellar gas temperature profile was assumed to be given by a power law of the format $$T(r)=T_0{\left(\frac{r}{R_{\text{in}}}\right)}^{\alpha },$$(3)

where T₀ is the gas temperature at the inner radius (R_in) of the models, and the slope, α, is a free parameter. We tried setting a lower limit of 10 K for the gas temperature in the outermost regions of the models, but it made no impact on the results. The fractional abundance of molecular species is defined with respect to molecular hydrogen (H₂), as $$f_X(r)=\frac{n_X(r)}{n_{{\text{H}}_2}(r)},$$(4)

where n_X (r) and n_H2 are the number density distributions of the molecule and H₂, respectively. It was implicitly assumed here that there is no photodissociation of H₂ within the spatial extents of the other molecules, which is valid since H₂ has the most spatially extended abundance distribution among all circumstellar species (Huggins & Glassgold 1982).

4.2 Dust modelling

We modelled the SEDs of our sources using MCMax (Min et al. 2009), a Monte Carlo RT code that takes into account absorption, re-emission, and scattering to solve the radiative transfer of dust. In addition to standard Monte Carlo RT, where photon packets and their interaction with dust grains are traced across the medium, MCMax includes an implementation of the diffusion approximation (Min et al. 2009) to incorporate energy transport by radiative diffusion (e.g. Wehrse et al. 2000) as well, which improves the accuracy of the estimated temperature profile in regions of high optical depth, and can also increase computational speed. The code assumes axial symmetry for the medium.

Key input parameters include the radiation field from the central star, the dust density distribution, and the optical properties of the dust. We used the spherical, hydrostatic COMARCS model atmospheres (Aringer et al. 2019) as input stellar spectra for IRAS 15194-5115, IRAS 15082-4808, and IRAS 07454-7112, with a surface gravity (log g) of −0.5, and a microturbulence velocity of 0.5 km s⁻¹. The COMARCS programme generates synthetic spectra, assuming local thermodynamic and chemical equilibrium, taking into account atomic and molecular lines (Aringer et al. 2016, 2019). These models thus serve as better representations of the atmospheric spectra of AGB stars than simple black bodies. The input COMARCS spectra were extrapolated and their resolution was reduced following the same procedure as described by Andriantsaralaza et al. (2022). For AFGL 3068 and IRC +10 216, we had to resort to black body spectra instead of the COMARCS models, as we could not obtain good fits using models in the available temperature range of COMARCS (2500-3300 K) and had to assume cooler stars.

We modelled the dusty envelope with a fixed inner radius of 3R_⋆, in line with both dynamical models and interferometric measurements, that place the dust condensation radius (R_c) for carbon stars in the range 2-3 R_⋆ (e.g. Andersen et al. 2003; Nowotny et al. 2011; Wittkowski et al. 2007; Karovicova et al. 2011; Norris et al. 2012; Bladh & Höfner 2012; Höfner & Olofsson 2018). We used amorphous carbon dust grains with a grain density of 2 g cm⁻³, typical of carbon stars (Höfner & Olofsson 2018). The optical properties of the dust were adopted from Suh (2000).

Testing different assumed grain sizes usually expected for AGB circumstellar dust (see e.g. Mattsson et al. 2010; Höfner & Olofsson 2018, and references therein), we concluded that only 1 μm grains led to satisfying results. Models with smaller grains all gave poor fits to the short-wavelength part of the SEDs, and/or led to unreasonably high dust condensation temperatures (T_c > 2000 K), as long as R_c was fixed at 3R_⋆. When we attempted to leave R_c as a free parameter instead of fixing it to 3R_⋆, models with grains smaller than 1 μm all returned unrealistically large R_c values of up to ~16-31 R_⋆, incompatible with dust-driven wind models. This is further discussed in Sect. 6.1.

The stellar luminosity (L_⋆), stellar temperature (T_⋆), and dust-MLR (Ṁ_d) were the free parameters in the SED fitting. For each star, an initial grid was constructed based on previously reported values of the free parameters from the literature, when available. The grids were then refined in the vicinity of local minima. A dust-MLR step of 0.5 dex was employed. The luminosity step was initially set to 1000 L_⊙ and subsequently refined to 500 L_⊙ when necessary. Distance estimates for all five stars were obtained from Andriantsaralaza et al. (2022). MCMax computes the SED, the optical depth at 10 μm (τ₁₀), and the temperature distribution of the dust for each grid point. The best-fit model for each source was identified by minimising the reduced chi-squared (${\chi }_{\text{red}}^2$), using $${\chi }_{\text{red}}^2=\frac{1}{(N-p)}\sum _{i=1}^N\frac{(F_{\text{mod},\lambda }-F_{\text{obs},\lambda }{)}^2}{{\sigma }_{\lambda }^2},$$(5)

where N is the number of flux measurements used, F_mod,λ is the modelled flux density at wavelength λ, F_obs,λ denotes the corresponding observed flux density (Table A.1) with uncertainty σλ (see Sect. 3.1), and p is the number of free parameters, equal to 3 in our models. The observed SEDs and the corresponding best-fit models are shown in Fig. A.1.

4.3 CO modelling

To ensure consistency in our modelling across the sources, we derived gas mass-loss rates by CO emission line modelling for each star. We used the one-dimensional (1D), non-LTE, Monte Carlo programme (MCP) RT code (Schöier & Olofsson 2001) to determine the CO molecular level populations across the CSEs. MCP has been extensively used for RT modelling of AGB CSEs (e.g. Danilovich et al. 2014, 2015; Ramstedt et al. 2008; Schöier et al. 2002, 2006, 2007, 2013). With the results from the dust modelling (T_⋆, L_⋆, τ₁₀ and T_dust(r); see Sect. 4.2) as input, along with basic source parameters (Table 1) and molecular data, MCP calculates the molecular excitation across different radial bins by using the Monte Carlo method to solve the equations of statistical equilibrium for each model in a grid of gas MLRs and temperatures.

We included rotational energy levels from J = 0-40 in both the ground (3 = 0) and first vibrationally excited (3 = 1) states of CO in our modelling, covering upper-level energies (E_up) from 5.5 K (3 = 0, J = 1) up to −7578 K (3 = 1, J = 40). We use rate coefficients for inelastic CO-H2 collisions provided by Yang et al. (2010), weighted assuming the statistical ortho-topara H₂ ratio of 3 (Danilovich et al. 2018; Massalkhi et al. 2024). The inner radius of the models was set to the dust condensation radius (R_c), equal to 3R_⋆. The fractional abundance of CO was calculated using $$f(r)=f_0\exp (-\text{ln}2{\left(\frac{r}{R_{1/2}}\right)}^s),$$(6)

where f₀ is the initial photospheric abundance of CO with respect to H₂, R_1/2 is the photodissociation radius defined as the radius where the CO abundance has dropped to f_0/2, and s is the slope of the profile. We adopted a typical canonical value of f₀ = 1 × 10⁻³ for all stars (e.g. Danilovich et al. 2015). For the sake of consistency, we did not leave the CO radial extent as a free parameter, as ALMA data that can constrain this parameter were only available for two of the sources in our sample. Instead, we used the values of R1/2 and s tabulated by Saberi et al. (2019), which are dependent on the gas mass-loss rate (M_g) and the terminal expansion velocity (ν_∞, see Table 1). The gas expansion velocity is described by Eq. (1), with ν_∞ taken from Woods et al. (2003). We assumed no drift velocity between the gas and the dust.

The kinetic temperature profile of the gas, assumed to be a power law (Eq. (3)), is varied in the grid with T₀ and α as free parameters, along with the gas-MLR, Ṁ_g. We created individual grids for each source, with T₀ in the range 500-2250 K, with a step-size of 250 K, α varying from −1.4 to −0.4, in steps of 0.2. We centred the Ṁ_g grids around previously determined gas-MLR values for the respective sources from the literature, scaled to the new distances from Andriantsaralaza et al. (2022), with a step of 0.2 dex, overall spanning over 1.6 units in log scale.

Model line profiles were generated for all the SD observations and for all beam sizes used to extract spectra from the ALMA line cubes (see Sect. 3.2). For each source, models with ${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2$ ≤ 5 were first identified. From among these, the best-fit models were then determined by χ² minimisation coupled with visual inspection of the modelled line profile shapes, as commonly done in the literature (see e.g. Ramstedt et al. 2008). This is done since the χ² values are based exclusively on the integrated line intensities and do not take into account the shape of the line profiles. The observed CO line profiles along with the best-fit modelled spectra are shown in Figs. B.1–B.4. We estimate an overall uncertainty of a factor of ~2-3 in our reported MLR values (see Sect. 6.1).

4.4 CS modelling

4.4.1 Molecular data

The CS molecular description employed in this work is the same as that used by Danilovich et al. (2018), comprising of rotational energy levels ranging from J = 0-40 within both the ground (ν = 0) and first excited (ν = 1) vibrational states, amounting to a total of 82 energy levels. The ν = 0 levels include 820 collisional transitions, and the ν = 0 and ν = 1 levels collectively account for 160 radiative transitions. The energy levels and radiative transition data were sourced from the Cologne Database for Molecular Spectroscopy (CDMS⁹; Müller et al. 2005; Endres et al. 2016). The collisional rate coefficients were adapted from CO-H2 collision calculations presented by Yang et al. (2010), assuming an H2 ortho-to-para ratio of 3. We also tested the more recent CS-H2 collisional rates from Denis-Alpizar et al. (2018), and found no significant difference between the modelled line intensities in the two cases. For the isotopologue models of ¹³CS and C³⁴ S, we used energy levels and transitions for the same range of states (J = 0-40, 3 = 0, 1) as ¹²CS, taken from the CDMS. We employed the same set of collisional rates as for the main CS isotopologue for both ¹³CS and C³⁴S.

4.4.2 Modelling procedure

We used a 1D accelerated lambda iteration (ALI) non-LTE spectral line RT code to model CS line emission. The code is based on the ALI method (Rybicki & Hummer 1991), and is described in detail by Maercker et al. (2008). ALI provides reliable convergence even at high optical depths, and has been widely used to model circumstellar line emission (e.g. Justtanont et al. 2005; Maercker et al. 2008, 2009; Schöier et al. 2011; Brunner et al. 2018; Danilovich et al. 2016, 2017, 2018, 2019). The code has been well benchmarked against other similar non-LTE line RT codes (see van Zadelhoff et al. 2002; Maercker et al. 2008). In this work, we use the latest available version of the ALI code, the same used by Danilovich et al. (2018, 2019) to model circumstellar molecular line emission.

ALI takes the results from the dust and CO RT models (see Sections 4.2, 4.3), i.e. the stellar luminosity and temperature, the kinetic temperature profiles of the circumstellar dust and gas, the dust optical depth at 10 μm, and the gas mass-loss rate, as inputs, and for each modelled line, calculates the molecular excitation at each radial point in the CSE. We used 40 radial points distributed logarithmically between Rin and Rout . New level populations are calculated at each iteration, and we assume that the model has converged once the average population change at each radial point between two successive iterations is less than 10⁻³. The above choices of the number of radial shells and the convergence threshold are well justified from the literature (see e.g. Maercker et al. 2008). Considering that CS is a parent species originating close to the stellar surface, we assumed a centrally-peaked Gaussian distribution to model the circumstellar CS fractional abundance (Eq. (4)) profile, given by $$f(r)=f_0\exp (-{\left(\frac{r}{R_e}\right)}^2),$$(7)

where f₀ is the peak abundance at the inner edge of the envelope, and R_e is the e-folding radius, at which the CS abundance has dropped to f₀/e. f₀ and R_e as the only free parameters in our modelling. The outer radius (R_out) of the models was set to five times the respective R_e, ensuring that the model covers a sufficient volume to recover all observed emission.

We designed individual grids for each source. We periodically checked the model outputs to ensure the grid coverage was sufficiently refined and extended over an adequate region of the parameter space. We used integrated line intensities from a number of progressively increasing apertures for the ALMA lines (see Fig. 2), along with line intensities from SD observations (Table 2) to constrain the models, as we did for the CO modelling. This way, we employed both the excitation information (from the different J transitions) and the azimuthally averaged radial distribution of the emission (from the ALMA maps) to constrain the CS models. The small apertures from the ALMA maps and the high-J SD lines give constraints on the inner wind, while the larger ALMA apertures and low-J SD lines constrain more the radial extent.

CS can be photodissociated by radiation having a broad range of frequencies from across the broadband continuum, unlike the case of CO where photodissociation can occur only due to radiation from specific UV lines (van Dishoeck 1988). Hence, CS has no differential self-shielding between its isotopologues (e.g. Hrodmarsson & van Dishoeck 2023), there is no reason to expect that the radial extent of the abundance profiles, which is determined by photodissociation in the outer envelope, should vary between the different isotopologues. So, for ¹³CS and C³⁴S, we used the same modelling procedure as described above for CS, with the exception that we now have only f₀ as a free parameter, with R_e being fixed to that of the corresponding best-fit CS model. The f₀ grids for isotopologue modelling were generated by scaling down the peak-abundance grids used for ¹²C³²S with the expected isotopic ratios for each source (see Sect. 6).

For each star, the uncertainties on the line intensities of spectra extracted from different ALMA apertures (see Sect. 3.3.1) are inherently correlated as they are taken from the same data cube. To account for this in the uncertainty estimates for the fit parameters (f₀ and R_e), we included these correlations in the χ² calculations, treating the ALMA spectra as correlated and the SD lines as independent observations. If we sort the lines with all the ALMA spectra first, followed by the SD ones, the resulting covariance matrix, C, is defined as $$\begin{array}{r}{C}_{ij}=\{\begin{array}{ll}{\sigma }_i^2,& i=j,\\ {\sigma }_i{\sigma }_j\rho ,& i\ne j,\text{\hspace{0.17em}with\hspace{0.17em}}i\text{\hspace{0.17em}and\hspace{0.17em}}j\le n_{\text{ALMA}},\\ 0,& i\ne j,\text{\hspace{0.17em}with\hspace{0.17em}}i\text{\hspace{0.17em}or\hspace{0.17em}}j>n_{\text{ALMA}},\end{array}\end{array}$$

where σ_i represents the uncertainty in the integrated line intensities, ρ is the correlation coefficient, and n_ALMA denotes the number of spectra extracted from the ALMA cube. We assume a correlation coefficient of 0.9, indicating that the spectra from different ALMA apertures are strongly correlated. We did not set ρ exactly equal to 1, since the radial substructure and complex morphology of the line emission would lead to a small degree of decorrelation among the spectra extracted from different apertures. We checked that minor changes to the correlation coefficient around the chosen value do not affect the results significantly, and only lead to marginal changes in the estimated uncertainties as well. The off-diagonal elements of C were set to zero for the SD lines. It was also ensured that the spatial information from the ALMA spectra and the excitation information from the different J lines both had the same total weight, to prevent the χ² minimisation from being strongly biased towards the spatial information simply due to the larger number of ALMA apertures compared to the number of available SD lines. We note that while such equal weighting is a sensible scheme in our context, there is no one true way to assign these relative weights in a general case, and an informed choice should be made based on the number of available ALMA and SD spectra.

For each model in the grid, the residual vector, r, is defined as the difference between the modelled (m) and observed (o) integrated line intensity vectors (r = m - o). The inverted covariance matrix (C⁻¹) was then used to calculate the χ² values as $${\chi }^2={\mathbf{r}}^{\mathrm{\top }}{\mathbf{C}}^{-1}\mathbf{r}.$$(8)

Line profiles from models with the lowest χ² values in the grid were visually inspected before finalising the choice of the best-fit model, for the same reasons as described in Sect. 4.3. For a fit with two free parameters, as is the case here, the 1σ, 2σ, and 3σ uncertainty levels for the best-fit parameters are given by ${\chi }_{\text{min}}^2+2.3,{\chi }_{\text{min}}^2+6.18$, and ${\chi }_{\text{min}}^2+11.83$, respectively, where ${\chi }_{\text{min}}^2$ is the minimum value of χ². We note here that the uncertainties estimated for f₀ and R_e are assuming the gas MLRs listed in Table 3, and depend on the choice of relative weighting between the different lines. The results from the CS modelling are presented in Sections 5.3 and 5.4.

Fig. 2 Comparison of modelled (red) and observed (black) integrated intensities of CS J = 2-1 spectra extracted from successively larger beams (see Sect. 4.4.2). The blue-shaded region represents the uncertainty in the observed line intensities.

5 Results

5.1 Dust

For each star, the best-fit model corresponds to the model with the lowest reduced chi-square (${\chi }_{\mathrm{r}\mathrm{e}\mathrm{d}}^2$) values. Furthermore, to ensure consistency with a dust-driven wind scenario (see Sect. 4.2), we only considered models that satisfy an upper limit constraint on the condensation temperature (Tc < 1900 K) of amorphous carbon dust (see Nanni et al. 2016; Höfner & Olofsson 2018). The SEDs of our sources peak at approximately 10 μm, as shown in Fig. A.1. The models reproduce the photometric data reasonably well for all stars. The most notable deviation between data and model occurs at the shortest wavelengths, corresponding to the Gaia data, which are consistently overestimated by the models for IRAS 15194-5115, IRAS 15082-4808, and IRAS 07457-7112, for which we used the COMARCS models as input stellar spectra (Sect. 4.2). However, due to the intrinsic variability of these stars, which is most prominent at short wavelengths, and the large contribution of dust emission in the SEDs at such relatively high dust mass-loss rates, the Gaia data points are not the most critical constraints. This is reflected in the high uncertainties assigned to them (see Sect. 3.1).

The three free parameters (Ṁ_d, T_⋆, and L_⋆) are well constrained for all stars except for IRAS 15194-5115, where T_⋆ and L_⋆ remain poorly constrained, resulting in degenerate viable models for that star. The dust mass-loss rates (Ṁ_d) derived across our sample span over an order of magnitude, from 1.26 × 10⁻⁸ M_ʘ yr⁻¹ for IRAS 07454-7112 to 1.26 × 10⁻⁷ M_ʘ yr⁻¹ for AFGL 3068, with an average value of 6.60 × 10⁻⁸ M_ʘ yr⁻¹. The results of our dust modelling are listed in Table 3.

5.2 CO

Our best-fit models to the CO line emission (see Figs. B.1–B.5) yielded well-constrained mass-loss rates for the studied sources, spanning from 8.3 × 10⁻⁶ to 4.2 × 10⁻⁵ M_⊙ yr⁻¹ (Table 3). To define the accelerating gas velocity profiles (Eq. (1)), we initially used β = 1 for all sources, as is commonly done in the literature (e.g. Massalkhi et al. 2024; Danilovich et al. 2015). This led to good models for all stars except IRAS 15194-5115, for which β = 1 resulted in very weak modelled lines for the high-J CO transitions J = 9-8 and J = 10-9. Hence, we tested different β values for this source, and found that only β = 5, implying slower acceleration compared to β = 1, produced significant emission for the high- J transitions (J ≥ 6-5). We then tested β = 5 for the other sources as well, and found no significant difference in the line intensities beyond the calibration uncertainties of the observations, compared to the β=1 case.

For IRAS 15194-5115, only SD observations were available, and a total of 11 transitions were used to constrain the models. The best-fit model for this source (with β = 5, see above) fits the data reasonably well overall, though some lines are slightly overestimated while others are underestimated (Fig. B.1). The most problematic case is the APEX J = 7-6 line, where the intensity ratio between model and data is 1.77. This line might be affected by calibration issues, and is also less reliable overall because of the relatively high noise level. High J-lines are systematically underestimated, indicating that the inner wind is not well constrained.

Similarly, only SD observations were available for IRAS 15082-4808, covering four transitions in total, up to the CO J = 4-3 line. In general, the fit is reasonable (Fig. B.2), as the APEX data and the SEST J = 1-0 lines match well with the model, although two of the SEST lines are not well reproduced. Since the APEX data are expected to be more reliable compared to SEST, we prioritised achieving better agreement with APEX rather than SEST.

IRAS 07454-7112 has the lowest gas MLR in the sample. Both SD and ALMA data were available for this star. The analysis includes transitions up to J = 14-13, with a total of six lines. The best model successfully reproduces the general shape of the lines (Fig. B.3). Although the fits are generally acceptable for the ALMA lines, with some lines being underestimated (J = 32) and others overestimated (J = 2-1), no model satisfactorily reproduces the J = 2-1 APEX line. For the CO J = 3-2 ALMA line, the integrated intensity decreases beyond 14.3″, likely due to resolved-out flux. For this star, high-J transitions are also systematically underestimated, suggesting that the inner wind is not very well constrained.

For IRC+10,216, only SD data were used as constraints, and the analysis included four transitions. While the line shapes are well reproduced, the modelled line strengths do not agree consistently with observations. For example, the model overestimates the NRAO CO 1-0 line intensity while reproducing the OSO and SEST CO 1-0 lines well. We note that these 1-0 SD observations should, however, be treated with caution as the telescope intensity calibrations can be highly uncertain for that transition (Fig. B.4).

Both ALMA and SD data were available for AFGL 3068. A total of five transitions were used, with the highest being J = 6-5. This star has the highest MLR in the sample, and its CO line profiles exhibit an irregular triangular shape, characteristic of high-MLR AGB stars. The best-fit model reproduces the line widths well. The CO J = 3-2 ALMA line shows signs of resolved-out flux, as beyond 15″, the integrated intensity decreases with radial distance. While ALMA lines show reasonably good agreement with the model, some APEX lines could not be reconciled with ALMA observations (see Fig. B.5).

5.3 CS

Models that reproduce the observed CS emission reasonably well were obtained for all five sources in the sample, yielding well-constrained abundance profiles (see Fig. 7). Figure 1 and Figures C.1–C.3 show the best-fit models and the observations, where we have chosen to present only a selection of the many extracted ALMA spectra. The comparison between modelled and observed integrated intensities for all extracted ALMA spectra is shown in Fig. 2. For the five modelled stars, we find f₀ in the range 1-4 × 10⁻⁶ and Re in the range 1.8-6.8 × 10¹⁶ cm (see Table 3).

For each grid of models, maps showing χ² as a function of the two adjustable parameters, f₀ and R_e, with the best-fit model and corresponding 1σ, 2σ, and 3σ uncertainty levels marked are shown in Fig. C.5. We verified that the best-fit models obtained by considering only the spatial information from the different ALMA apertures, and only the excitation information from the different transitions (the largest aperture ALMA spectrum and the SD lines), are consistent within their respective uncertainties, and show no signs of significant divergence (see Figs. 3, C.4). The results for each source are summarised below.

Overall, we obtain good fits to the CS line profiles for IRAS 15194-5115 (Fig. 1), though the quality of the line intensity fits is not as good in comparison to those we get for the other stars in the sample. The precision with which the CS abundance can be constrained in the innermost regions of the CSE is limited by the large uncertainties on the high J HIFI line intensities (Sect. 3.3.3). We note that the underlying velocity profile used is also different from that of the other sources for this star (see Sect. 5.2). Our best-fit model reproduces the line profiles well, and gives a peak abundance of 3.1X10⁻⁶ and an e-folding radius of 6.8X10¹⁶ cm. However, this model underestimates the CS J = 2-1 line intensities from apertures smaller than ~7″ by ~10-50%, as seen from Fig. 2. Another model with f₀ = 4.4×10⁻⁶ and R_e = 4.6×1016 cm, that falls within the 1σ interval of the above-mentioned best-fit model manages to reproduce the line intensities from the inner part of the envelope better, but slightly underestimates the larger aperture ALMA spectra, while also significantly overestimating the higher- J HIFI spectra (J = 1817 to 22-21). Our best-fit model also overestimates the APEX CS J = 7-6 line and a few of the high excitation HIFI lines, mainly J = 18-17 and J = 21-20 (see Fig. 1). However, we obtain very good fits to the observed line profiles for ALMA apertures from ~7″ onwards, up to around 35″, which encompasses all the detected emission for the J = 2-1 line (Fig. 2), as well as most of our other APEX and HIFI lines (Fig. 1). The complexity in modelling this source is further discussed in light of previous related works in Sect. 6.

For IRAS 15082-4808 and IRAS 07454-7112, we obtain models that reproduce the line profiles from both the various ALMA apertures and the SD lines (see Figs. 2, C.1a, C.1b), with the exception of the J = 7-6 line for IRAS 07454-7112, which our model underestimates (Fig. C.1b). Though our model for IRAS 15082-4808 also underestimates the J = 4-3 line (Fig. C.1a), the modelled spectra for this line falls within the uncertainty range of the observed spectrum (Sect. 3.3.2).

For AFGL 3068, we are limited to SD spectra (Fig. C.2), as we do not have spatially resolved data available. Though we have the CS J = 7-6 line for this source observed with the ACA, we cannot extract much information about the radial distribution of the emission, as the CS emitting region is similar in size to the ALMA synthesised beam of the observation. We consequently treated this line as an SD observation in this work. The lack of spatial information and the small number of SD lines available imply that we cannot place very rigorous constraints on the abundance profile for this source as we have done in the case of the other four stars. This is reflected in the larger uncertainties of the estimated parameters (Table 3) for this source.

For IRC +10 216, we discarded the smallest ALMA aperture (see Fig. C.3) from the χ² minimisation, as it displays an asymmetric, secondary peak around the systemic velocity, which could be due to small scale asymmetries in the emission that cannot be modelled by our smooth, 1D, spherically symmetric RT code. We obtained a good fit to the rest of the spectra and were able to simultaneously fit the cumulative radial emission profile from the ALMA CS J = 2-1 line (Fig. 2) and the SD line intensities. While the APEX CS J = 5-4 line from our spectral survey and the IRAM 30 m J = 3-2, 4-3, 5-4, and 6-5 lines yield modelled spectra that match the observations within the corresponding uncertainties, our model overestimates the Yebes J = 1-0 line. We also note that the calibration uncertainties for the IRAM 30 m telescope rise significantly with frequency (see Sect. 3.3.4).

5.4 ¹³CS and C³⁴S

Our models reproduce the observed line profiles for ¹³CS and C³⁴S for all five stars in the sample. For IRAS 15194-5115, IRAS 15082-4808, and IRAS 07454-7112, we have spatially resolved J = 2-1 ALMA maps for both isotopologues, so we were able to also model multi-aperture spectra (Sect. 3.3.1) for these sources (Figs. D.1, D.6), along with the available SD lines. For AFGL 3068 and IRC +10 216, we only had SD spectra for the isotopologues. The best-fit models and observed spectra are shown in Figs. D.2–D.5 for ¹³CS, and Figs. D.7–D.11 for C³⁴S. The best-fit f₀ values obtained are given in Table 4.

We note that for AFGL 3068, though we get good fits to the isotopologue line profiles using the peak abundance listed in Table 4 and the R_e value from the best-fit ¹²CS model (Table 3), there are also other models, around f₀ = 2.3 × 10⁻⁸ and R_e = 6.0 × 10¹⁶ cm (larger than the best-fit ¹²CS R_e by a factor of ~3), which yield equally good fits. Some models with this R_e also fall within the 1σ uncertainty range from our best-fit ¹²CS model (see Fig. C.5d), but the quality of our ¹²CS fits there is considerably lower than that of our best-fit model. This indicates that our R_e estimate for AFGL 3068 may be underestimated, and the f₀ value correspondingly overestimated, which is also evident from the rather degenerate χ² contours in Fig. C.5d.

We also calculated the ¹²CS/¹³CS and C³²S/C³⁴S ratios from the peak abundances of the best-fit models (Table 4). The estimated isotopic ratios for AFGL 3068 should be regarded with caution as they suffer from the fact that the models are substantially less constrained compared to the other sources (Sect. 5.3).

6 Discussion

6.1 Dust and CO

As noted in Sect. 4.2, we fixed the dust condensation radius (R_c) at 3R_⋆ in our dust modelling, since, if left as a free parameter, its value becomes unreasonably large for a dust-driven wind. This issue can be seen in several previous works as well. Ramstedt & Olofsson (2014) and Ramstedt et al. (2008) report values up to ~46 au for M- and S-type stars and up to ~2θ au for C-type stars from dust RT models where the inner radius was a free parameter. In general, using a grain size distribution instead of a single fixed grain size may help alleviate this problem, while also approximating the circumstellar dust more realistically. This is, however, outside the scope of this work. Our interest in the dust here amounts to only finding a reasonable estimate of the dust radiation field to serve as input to our CO and CS RT models, rather than undertaking a comprehensive analysis focussed on the dust properties themselves. We obtained reasonable dust condensation temperatures in the range of ~1500-1800 K, which is also in line with the ~300 K uncertainty in T_c estimated by Nanni et al. (2016) for amorphous carbon dust.

Our estimates of the gas MLRs from CO RT modelling show very good agreement with those from the literature for all five sources, in most cases within a factor of ~2-2.5 (e.g. Ramstedt & Olofsson 2014; Danilovich et al. 2015; Schöier et al. 2007; Woods et al. 2003). Considering this and the general predictions on the uncertainties of AGB MLR estimates (factor of ~3, see e.g. Ramstedt et al. 2008) from the literature, we estimate an uncertainty of a factor of ~2-3 in our reported MLR values, assuming the input dust properties and CO abundance. The gas kinetic temperature profiles obtained are also well constrained, with all reasonable models in the grid falling within ±1 step size (~250 K for T₀, 0.2 dex for α, see Sect. 4.3) of the best-fit model for all sources except IRAS 15194-5115, and within ±2 step sizes of the best-fit model for IRAS 15194-5115. We note that these uncertainties on the MLRs and gas temperature profiles do not take into account the errors in the underlying dust properties and assumed canonical CO abundance, which cannot be realistically quantified.

As was noted in Sect. 5.2, we are unable to strictly constrain the gas velocity profile exponent, β, for sources other than IRAS 15194-5115. Observations of high-J line emission are needed to properly characterise the gas velocity in the innermost regions of the CSE. Likewise, the gas kinetic temperature profiles are also uncertain in the innermost regions of the winds. For IRAS 15194-5115 and IRAS 07454-7112, where we have sufficiently high-J HIFI observations, up to J = 1615 and J = 14-13, respectively, all otherwise reasonable models in our grids systematically underestimate these lines. For the other sources, no high-excitation data are available to constrain the gas temperature in these regions. Nevertheless, the existing observations provide good constraints on the gas temperature profiles across most of the CSE, except for the innermost regions. High-angular-resolution observations that resolve the inner wind (e.g. Khouri et al. 2024; Fonfría et al. 2014, 2019), and deep observations of narrow molecular lines from vibrationally excited levels, originating in regions close to the star (e.g. Patel et al. 2008) are needed to constrain the physical properties of the innermost parts of these CSEs.

Andriantsaralaza et al. (2021) measured the extent of the CO emission using ALMA observations for several carbon stars including IRAS 07454-7112 and AFGL 3068. The measured extent of CO emission and the CO photodissociation radius (see Sect. 4.3) can differ based on the excitation conditions in the CSE. These quantities will be compared for a number of carbon stars from the DEATHSTAR sample (Ramstedt et al. 2020; Andriantsaralaza et al. 2021) in an upcoming paper (Andriantsaralaza et al., in prep.).

6.2 CS

In this work, we employ the combined modelling of interferometric and SD data, to put rigorous constraints on the CS models using information about both the radial emission distribution and the conditions traced by different transitions simultaneously, for all sources in our sample except AFGL 3068 for which we have no spatially resolved data available. As was mentioned in Sect. 5.3, we reproduce the observed CS line profiles for all sources, and are able to robustly constrain their CS abundance profiles, with the exception of AFGL 3068, which presents comparatively large uncertainties. We find that CS maintains a significant gas-phase abundance relatively far out into the CSE, up to ~3-7×10¹⁶ cm for all five sources. This is consistent with previous predictions from observations of IRC +10 216 (Agúndez et al. 2012; Massalkhi et al. 2019; Velilla-Prieto et al. 2019). No significant depletion onto the dust is seen. The e-folding radius increases with MLR as expected. The size of the molecular envelope of CS is set by its photodissociation into C and S by the ambient interstellar UV radiation field. We find no discernible trend in the CS peak abundances of our sources with envelope density beyond the estimated uncertainties. A correlation between CS abundance and CSE density was reported for S-type stars by Danilovich et al. (2018), who also, however, did not detect any such trend for carbon stars.

The CS envelope sizes and peak abundances estimated by Woods et al. (2003) for our five sources using SD observations of CS J = 2-1 and J = 5-4 and simple photodissociation models match fairly well with those from this work. Our CS peak-abundance values are also consistent with estimates derived from large samples of C-type stars (e.g. Massalkhi et al. 2019). We also note that the LTE CS abundances calculated by UK24 using population diagrams are a factor of 2.5-6.5 higher than the non-LTE abundances estimated in this work. This highlights that LTE calculations can serve as simple, useful tools for obtaining first-order abundance estimates, but stricter constraints from non-LTE RT models are required for more precise calculations.

Our best-fit model for IRC +10 216 has a peak CS abundance approximately double than that calculated by Massalkhi et al. (2024), while the e-folding radii from the two estimates match within the 1σ uncertainty level. The difference found in peak abundance is due to the fact that Massalkhi et al. (2024) uses an MLR greater than ours by a factor of ~1.6. We have tested the dependence of the derived peak abundances on the input gas MLR for our models, and found that varying the MLR by a factor of two on either side leads to roughly a factor of two change in the best-fit abundance in the opposite direction, indicating an inverse linear relationship as expected. CS in two of our sources (IRAS 07454-7112 and IRAS 15194-5115) was also modelled by Danilovich et al. (2018) using APEX observations. For IRAS 07454-7112, we reduce the relative uncertainty in the CS e-folding radius by a factor of ~2.5 with respect to Danilovich et al. (2018), and further provide well-constrained uncertainties on the peak abundance as well.

IRAS 15194-5115 is the source with the largest number of observed CS lines among the sample, owing to the lines from the HIFI survey (Sect. 3.3.3). It is known from the literature to be a difficult source to model, possessing a complex circumstellar morphology, with extensive substructure seen in molecular line emission (see UK24), possibly caused by the presence of a binary companion (e.g. Feast et al. 2003; Lykou et al. 2018). Danilovich et al. (2018) reported that no satisfactory fit was found for the CS and SiS line profiles of this star using their 1D, constant MLR, smoothly accelerating RT model, which yielded good fits for the same species towards a number of other C-, M-, and S-type AGB stars. Our models with an underlying slow wind acceleration (β = 5, instead of 1, see Sect. 5.2) gave acceptable fits to the CS line profiles of IRAS 15194-5115, though as mentioned in Sect. 5.3, the overall quality of the fits is lower compared to those obtained for the other sources, with difficulties in fitting the inner wind/high-J and extended/low-J emission simultaneously, similar to that reported by Danilovich et al. (2018).

However, owing to the large number of lines available to model for IRAS 15194-5115, it is clear that even though a few of the SD spectra are sporadically overestimated (Sect. 5.3), there appears to be no trends showing deviations consistently correlated with line excitation. Hence, the intermittent mismatches found between observations and models may be attributed, as indicated by Danilovich et al. (2018), to the complex morphology of the gas distribution in the CSE and its effects on line excitation, which cannot be accommodated in our RT model. The availability of spatially resolved observations across a large range of energy levels to constrain RT models alone will not be sufficient to tackle this issue. Also, using more complicated step functions in the inner envelope instead of typical Gaussian abundance profiles, as done by Danilovich et al. (2019) for modelling SD observations of M-type stars, might not be completely straightforward in a case where spatially resolved observations sensitive to inner abundance changes are also available. To properly model this source would require RT models employing a more realistic representation of the underlying observed complex gas density structure, and chemical models that take into account the effects due to discrete variations in density and potential binarity (e.g. Cordiner & Millar 2009; Van de Sande & Millar 2022). The possibility of employing 3D RT modelling (e.g. Homan et al. 2017; Khouri et al. 2024), constrained by morphological information from the ALMA lines also remains to be explored.

Fig. 3 χ2 values for our grid of CS models for IRC +10 216. The contours denote 1σ, 2σ, and 3σ ranges, when constrained using only the spatial information from the CS J = 2-1 line (green) vs just the excitation information from the different J lines (blue). Also shown are the contours when using both the spatial and excitation information simultaneously (red) to constrain the RT models. The green, blue, and red stars mark the best-fit models estimated from the above three cases, respectively.

6.2.1 Spatial versus excitation constraints

Fig. 3 shows the difference in χ² values obtained for our grid of CS RT models for IRC + 10 216, when constrained using only the spatial information and only the excitation information separately, and also when using both simultaneously. This shows how well the spectra from various transitions and the spectra from different spatial apertures from the interferometric maps manage to independently constrain the CS abundance profiles. Similar figures for the other sources are shown in Fig. C.4. The general pattern is similar for all sources that while the best-fit models obtained from both cases fall sufficiently close to each other, the excitation constraints do not manage to strongly break the degeneracy between the two free parameters, resulting in the classic diagonally elongated confidence contours in the χ² maps (see Figs. 3, C.5d). The spatial constraints from the resolved observations, however, limit the overall fit to a narrower range of the parameter space, leading to better estimates of the uncertainties on the best-fit values. This highlights the advantage of incorporating constraints based on the spatial distribution of the line emission to refine RT models.

Fig. 4 Isotopic ratios from this work compared to estimates from the literature.

6.2.2 Isotopic ratios

Fig. 4 shows the ¹²CS/¹³CS and C³²S/C³⁴S ratios obtained from this work (see Sect. 5.4 and Table 4) in comparison to various ¹²C/¹³C and ³²S/³⁴S isotopic ratio estimates from the literature. For all sources except AFGL 3068, the isotopologue abundance ratios we obtain match the isotopic ratios from the literature (Cernicharo et al. 2000; Woods et al. 2003; Ramstedt & Olofsson 2014; He et al. 2008; Danilovich et al. 2018, UK24). The derived ¹²CS/¹³CS ratios are hence good measures of the ¹²C/¹³C isotopic ratio, suggesting that isotopologue-selective effects are likely not significant in CS chemistry. We likely overestimate the ¹²C/¹³C and ³²S/³⁴S ratios for AFGL 3068, due to the poorly constrained abundance profiles, which may be overestimating the ¹²C³²S peak-abundance (Sect. 5.4).

IRAS 15194-5115 stands out from the rest of the sample with its exceptionally low 12C/¹³C ratio of ~7 (see Fig. 4), which is reproduced by our ¹²CS/¹³CS abundance ratio as well. Ramstedt & Olofsson (2014) note that it is possible for J-type AGB stars with low MLRs to have ¹²C/¹³C ratios as low as 3. Though IRAS 15194-5115 is a J-type star (Smith et al. 2015), it has a high MLR (see Table 3; and Woods et al. 2003). A reason for the very low ¹²C/¹³C ratio can be hot-bottom burning (HBB), a process in which the bottom of the convective envelope becomes hot enough to start nuclear burning, destroying ¹²C and producing small amounts of ¹³C. Stars in the luminosity range −6.4 < M_bol < −6.3 may retain enough ¹²C to remain C-type even after undergoing HBB, resulting in such low ¹²C/¹³C ratios (Boothroyd et al. 1993). From this work, we however find a bolometric magnitude of −5.3 (Table 3) for IRAS 15194-5115, in line with previous estimates (e.g. Woods et al. 2003; Lykou et al. 2018) as well. It is hence unclear why IRAS 15194-5115 has a peculiarly low ¹²C/¹³C ratio, though possible explanations could include differential photodissociation due to the presence of a hot binary companion (Vlemmings et al. 2013), extra mixing during the early-AGB phase or during the helium-core flash, or processes like cool bottom processing (see Ramstedt & Olofsson 2014, and references therein).

Our estimated C³²S/C³⁴S ratios are very close to the solar ³²S/34S ratio (~22, Asplund et al. 2009). This is expected, as ³²S and ³⁴S are not nucleosynthesised at any stage of low/intermediate-mass stellar evolution up to the AGB, and hence their ratio does not change significantly before or during the AGB phase (e.g. Karakas & Lugaro 2016). These isotopes are mostly produced in supernovae, and also by oxygen-burning (e.g. Hughes et al. 2008). It is also interesting to note that the CS line intensity ratios (see Fig. 4) are always less than the corresponding abundance ratios from RT modelling (with the exception of ¹²CS/¹³CS for IRC+10 216). This is due to the high optical depths in the ¹²C³²S lines, which lead to suppressed line intensities.

6.2.3 A note on CS excitation: IR pumping

IR pumping refers to the phenomenon where molecules are excited to higher vibrational states by the absorption of IR photons from the dust or the central star, and then spontaneously decay to rotational levels in the ground vibrational state, increasing the populations in these levels. This can increase the line intensities of the ν = 0 rotational lines, and has been shown to contribute significantly to the excitation of several species in CSEs (e.g. Agúndez & Cernicharo 2006; Agúndez et al. 2012). Using only the ν = 0 levels is equivalent to not including IR pumping in the models. To check the importance of infrared (IR) pumping in the excitation of CS in our models, we compared the ratio of the excitation temperatures to the gas kinetic temperature and the line emitting regions obtained when both ν = 0 and 1 levels are considered in the level population calculations, to when only the ν = 0 levels are used. The T_ex/T_kin ratio is shown in Fig. 5 for the best-fit CS model for IRAS 15194-5115, for various low- and high-J transitions, along with their respective emitting regions. The emitting region for a line in the RT model is calculated by multiplying the brightness emitted from each radial shell (see Sect. 4.4.2) by the square of the radius at which the corresponding shell is located. The emitting region moves inwards as we go to higher-J lines as expected, as these lines are excited at higher temperatures and densities compared to their lower-J counterparts.

It is seen that IR pumping affects both the excitation and radial emission distribution in our models. The emitting regions become more compact with their peaks shifted inwards when IR pumping is turned off (Fig. 5). Not considering IR pumping also leads to a reduction in the observed line intensity for all the modelled CS lines by ~50-70%. This would lead to overestimating the peak abundances and e-folding radii, showing that it is crucial to include vibrationally excited levels in the calculation of level populations when modelling CS emission. The fact that the modelled line intensities did not change significantly upon testing different collisional rates (see Sect. 4.4.1) can also be indicative of the dominant role of radiative excitation in the formation of these lines. In the warm and dense inner part of the envelope, we have T_ex ≈ Tkin, i.e. the rotational levels are ther-malised (Fig. 5). Suprathermal excitation (T_ex > T_kin) is found for all modelled lines as we move farther away from the star, which can be attributed to IR pumping becoming significant as the IR-emitting dust content builds up. At even larger radii, as the gas density decreases, the excitation eventually becomes subthermal (T_ex < T_kin) for all the studied lines. These results are in line with the findings of Massalkhi et al. (2019) from modelling CS J = 3-2 line emission for a large set of carbon stars.

Fig. 5 Ratio of the excitation temperature (Tex) and the kinetic temperature (Tkin) for various low- and high-J CS transitions (top), and the corresponding line emitting regions (bottom), normalised to the peak of the J = 21-20 emitting region, as functions of radius across the CSE, from the best-fit model for IRAS 15194-5115. See text for a more detailed description.

6.2.4 Comparisons with chemical models

The correlation between the radial extent of molecular abundance profiles and envelope density has been analysed for AGB CSEs by several authors based on RT models (Danilovich et al. 2018; Massalkhi et al. 2019, 2024) and chemical models (Maes et al. 2023). This has led to several empirical relations describing the increase in the R_e of various molecules including CS with increasing circumstellar density, represented by Ṁ/ν_∞, where Ṁ is the gas MLR and ν_∞ is the terminal expansion velocity of the gas. In Fig. 6, we explore how the R_e estimates from our models compare to empirical trends with Ṁ/ν_∞ from the literature. We also show the uncertainty in the CS e-folding radii obtained from the chemical model, calculated by performing an uncertainty analysis where the uncertainties of all reaction rates in the RATE22 chemical network (Millar et al. 2024) were propagated through the chemical model using a Monte Carlo sampling method for $\dot{M}=1.0\times {10}^{-5}{M}_{\odot }y{r}^{-1}$ and ${\nu }_{\mathrm{\infty }}=15$ km s⁻¹ (Van de Sande et al., in prep.). Within the overall uncertainties, our results follow the trends from the literature well, and agree with the general trend where denser outflows show larger R_e, indicating larger molecular envelopes. Based on these trends as seen from Fig. 6, we expect the e-folding radius for the CS abundance profile of AFGL 3068 to increase when constrained better using additional observations, as also noted before (Sect. 5.4).

In order to compare our retrieved radial extents of CS abundance profiles directly with predictions from chemical models, we ran chemical kinetic models specific to our sources using the publicly available, state-of-the-art UMIST Database for Astrochemistry (UDfA¹⁰, Millar et al. 2024) circumstellar chemical model¹¹. Previous versions of this chemical model have been extensively employed to model various aspects of AGB circumstellar chemistry (e.g. Van de Sande & Millar 2019, 2022; Van de Sande et al. 2019, 2020, 2021). In this work, we used the recent RATE22 UDfA update (Millar et al. 2024), that includes updated rate coefficients for a total of 8767 reactions among 737 species. The 1D chemical kinetics code uses these rates to model the gas-phase chemistry in AGB outflows and derive molecular abundance profiles. It assumes a constant mass-loss rate and gas velocity, leading to a uniformly expanding, spherically symmetric CSE. The gas density in the envelope is assumed to fall as 1/r², where r is the distance from the centre of the star. For a detailed description of the chemical model, we refer to Millar et al. (2000); Cordiner & Millar (2009) and Van de Sande et al. (2018b).

We adopted parent species abundances compiled by Agúndez et al. (2020) except for CS, with the addition of metals (Na, Mg, Fe) with initial abundances relative to H2 set to 1.0×10⁻⁸, as a source of electrons. For each of our sources, we used the peak CS abundance obtained from our best-fit RT models as the initial abundance of CS in the chemical model. Other relevant source properties, including gas MLR and ν_∞ were taken from Tables 1 and 3. We tested both the Pattillo et al. (2018) CS photodissociation rate, included in RATE22 and also the new, lower, rate calculated from the cross sections from Hrodmarsson & van Dishoeck (2023), based on Xu et al. (2019). The resulting CS abundance profiles corresponding to the two different photodissociation rates are shown in Fig. 7, along with those obtained from RT modelling. We found that varying the input initial CS abundance to the chemical model within the 1σ range obtained from the RT models did not cause any significant change to the predicted R_e values.

We find that the chemical model abundance profiles match our CS abundance profiles from RT modelling within the estimated uncertainties. For three out of the five sources (IRAS 15194-5115, IRAS 15082-4808, and IRC+10 216), both the photodissociation rates match the retrieved RT abundance profiles within 1σ intervals. The deviation between the RT and chemical model Re values is large in the case of AFGL 3068, owing to its poorly constrained RT models, that likely underestimate Re (Sect. 5.4). In the case of IRAS 07454-7112, we find that the chemical model using the CS photodissociation rate based on Hrodmarsson & van Dishoeck (2023) matches our retrieved abundance profile, while that with the higher Pattillo et al. (2018) rate overestimates the CS radial extent. Fig. 7 also shows the uncertainty on the chemical model envelope sizes due to uncertainties in the kinetic data, calculated by Van de Sande et al. (in prep.), as described above. We note that this uncertainty on the predicted envelope sizes does not vary significantly with density, though lower densities show slightly larger error bars. This uncertainty calculation was done using the Pattillo et al. (2018) rate. Assuming that the error bar on the predicted envelope size does not vary significantly when changing the CS photodissociation rate, we can say that abundance profiles based on the two rates both do a reasonable job overall at reproducing CS abundance profiles obtained by RT modelling, and match within the estimated uncertainty on the chemistry. Further, at high envelope densities, as is the case with our sources, photodissociation is often attenuated by the high dust extinction, irrespective of the employed photodissociation rate (e.g. Massalkhi et al. 2024).

The overestimation of CS radial extents by their photodissociation model compared to their RT modelling results described by Massalkhi et al. (2024) using the Pattillo et al. (2018) rate could at least partly be due to the fact that their empirical fit to the CS e-folding radii includes both C-rich and O-rich outflows together. CS may have a slightly larger radial extent in C-rich envelopes compared to the O-rich ones, due to additional self-shielding against photodissociation due to its higher abundance in C-rich CSEs (e.g. Danilovich et al. 2018). We note that while the location of R_e is determined by photodissociation, its uncertainty also depends on the chemistry reforming the parent after photodissociation (Van de Sande et al., in prep.). To conclude, this implies that in the case of a parent species such as CS, where different photodissociation rates exist in the literature (e.g. Pattillo et al. 2018; Xu et al. 2019), comparisons between the different photodissociation models should be done with caution, until proper uncertainties have been estimated.

Finally, we note that this study utilises observations from a large number of telescopes, each of which has its own calibration accuracy and uncertainties. This heterogeneous origin of the data used poses an inherent limitation to the precision with which conclusions can be drawn from the presented analysis. We have attempted to accommodate this as far as possible, by taking into account the estimated calibration uncertainties of the respective telescopes when assigning error bars on the line intensities. These uncertainties indeed propagate to the final estimates of the physical parameters and molecular abundances as well. The only way to make these calculations more precise is to avoid uncertainties in relative calibration by using the same telescope as far as possible, though this comes with obvious practical limitations, particularly given the crucial need for spatially resolved data.

Fig. 6 Variation in e-folding radius of the CS abundance profiles with circumstellar gas density. The dashed cyan line (Maes et al. 2023) is from chemical modelling results, while the other trends plotted are calculated from RT models. The magenta point is obtained from chemical modelling results (Van de Sande et al., in prep.), for a star of MLR 1.0×10-5 Mʘ yr−1 and a ν∞ of 15 km s−1, and the corresponding uncertainty range is estimated from sensitivity calculations on chemical models (see text). The rightmost blue circle is for AFGL 3068, and is likely underestimated (see text).

Fig. 7 Comparison of CS abundance profiles obtained from RT modelling (blue) and chemical modelling: red - using CS photodissociation rate from Pattillo et al. (2018), green: using CS photodissociation rate calculated from the cross sections from Hrodmarsson & van Dishoeck (2023), based on Xu et al. (2019). The shaded blue region represents the 1σ uncertainty range in the RT model abundance profiles. The circles denote the Re of the respective models. The red and green error bars on the chemical model Re denote the estimated uncertainty in the chemical models (Van de Sande et al., in prep.). It is assumed that this uncertainty does not vary significantly with the outflow density and photodissociation rate (see text).

7 Summary and conclusions

We modelled the SEDs and CO line emission of our sample of five carbon-rich AGB stars, IRAS 15194-5115, IRAS 15082-4808, IRAS 07454-7112, AFGL 3068, and IRC +10216, to consistently estimate their stellar parameters, dust radiation fields and mass-loss rates. With these results as inputs, we employed detailed non-LTE radiative transfer modelling to constrain the abundance profiles of CS across the five sources. We utilised both spatially resolved ALMA observations of the CS J = 2-1 line as well as a broad range of SD lines from APEX, Herschel/HIFI, and the IRAM 30 m telescopes to simultaneously set spatial and excitation constraints on our RT models. In general, our models reproduce the observed CS line emission well and yield robust estimates of the CS abundance profiles. The heterogeneous origin of the observations used is a major source of uncertainty in the analysis. The calculated peak abundance and e-folding radius have comparatively larger error bars for AFGL 3068, due to the lack of spatially resolved observations, leading to less rigorous constraints on the models. We also determined the ¹²CS/¹³CS and C³²S/C³⁴S isotopologue abundance ratios, which are in good agreement with the expected ¹²C/¹³C and ³²S/³⁴S isotopic ratios for our sources.

Our estimates of the radial extents of the CS abundance profiles follow the expected trends with envelope density. We find no significant difference beyond their respective uncertainties between the calculated CS peak abundances for the studied sources. Our results are also in good agreement with predictions of the radial extent of CS abundance profiles from advanced chemical models, though we cannot discern between the different predicted photodissociation rates of CS beyond the uncertainties involved. These results highlight the importance of leveraging both spatial and excitation information simultaneously in fine-tuning radiative transfer models. A similar analysis of other molecules, especially the parent species, informed by high-resolution observations spanning a large range of excitation conditions is the next step required to advance our understanding of the chemistry in AGB CSEs.

Acknowledgements

The authors sincerely thank the anonymous referee for their constructive feedback, which improved the quality and clarity of this paper. We wish to thank P. Bergman for his support in setting up the ALI RT code. We thank M. Agúndez for providing Yebes 40 and IRAM 30 m CS line profiles of IRC +10216, L. Velilla-Prieto for answering our queries regarding the combined ALMA - IRAM 30 m CS J = 2-1 data, and J. H. Black for discussions about collisional rates and molecular data files. We are grateful to D. Tafoya for his continuous support in ALMA data related tasks. RU acknowledges data reduction support from the Nordic ALMA Regional Centre (ARC) node based at Onsala Space Observatory (OSO), Sweden. The Nordic ARC node is funded through Swedish Research Council grant No. 2017-00648. MA acknowledges support from the Olle Engkvist Foundation under project 229-0368. EDB acknowledges financial support from the Swedish National Space Agency. TD is supported in part by the Australian Research Council through a Discovery Early Career Researcher Award (DE230100183). MVdS acknowledges support from the Oort Fellowship at Leiden Observatory. TJM’s research at QUB is supported by grant ST/T000198/1 from the STFC. The work of MGR is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00070.S, ADS/JAO.ALMA#2015.1.01271.S, ADS/JAO.ALMA#2017.1.00595.S, ADS/JAO.ALMA#2013.1.00432.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This paper is based on observations with the Atacama Pathfinder EXperiment (APEX) telescope. APEX is a collaboration between the Max Planck Institute for Radio Astronomy, the European Southern Observatory, and the Onsala Space Observatory. Swedish observations on APEX are supported through Swedish Research Council grant No. 2017-00648. The APEX observations were obtained under project numbers O-0107.F-9310 (SEPIA/B5), O-0104.F-9305 (PI230), and O-087.F-9319, O-094.F-9318, O-096.F-9336, and O-098.F-9303 (SHeFI). HIFI has been designed and built by a consortium of institutes and university departments from across Europe, Canada, and the United States (NASA) under the leadership of SRON, Netherlands Institute for Space Research, Groningen, The Netherlands, and with major contributions from Germany, France and the US. Consortium members are Canada: CSA, U. Waterloo; France: CESR, LAB, LERMA, IRAM; Germany: KOSMA, MPIfR, MPS; Ireland: NUI Maynooth; Italy: ASI, IFSI-INAF, Osservatorio Astrofísico di Arcetri-INAF; The Netherlands: SRON, TUD; Poland: CAMK, CBK; Spain: Observatorio Astronómico Nacional (IGN), Centro de Astrobiología (INTA-CSIC); Sweden: Chalmers University of Technology - MC2, RSS & GARD, Onsala Space Observatory, Swedish National Space Board, Stockholm University - Stockholm Observatory; Switzerland: ETH Zurich, FHNW; USA: CalTech, JPL, NHSC. This work is based on observations carried out with the IRAM 30 m and the Yebes 40 m telescopes. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). The Yebes 40 m telescope at Yebes Observatory is operated by the Spanish Geographic Institute (IGN, Ministerio de Transportes, Movilidad y Agenda Urbana). This work has made use of data from the European Space Agency (ESA) mission Gaia¹², processed by the Gaia Data Processing and Analysis Consortium¹³ (DPAC). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Two Micron All Sky Survey¹⁴ (2MASS), which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work has made use of observations with AKARI¹⁵, a JAXA project with the participation of ESA. This research uses data from the Diffuse Infrared Background Experiment¹⁶ (DIRBE) instrument on the Cosmic Background Explorer (COBE). NASA Goddard Space Flight Center developed the DIRBE data sets under the direction of the COBE Science Working Group. This publication makes use of data from the IRAS point source catalogue¹⁷. The Infrared Astronomical Satellite (IRAS) was a joint project of the US, UK and the Netherlands. This research has made use of the International Variable Star Index¹⁸ (VSX) database, operated at AAVSO, Cambridge, Massachusetts, USA. This research has made use of NASA’s Astrophysics Data System¹⁹ (ADS). This project has made use of the VizieR²⁰ catalogue access tool and the SIMBAD²¹ database, operated at CDS, Strasbourg, France. This work made use of Astropy:²² a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022). This work has made use of GILDAS²³ and CASA²⁴ software to reduce and analyse data. The computations in this work were run in the Vera HPC cluster using resources provided by C3SE, the Chalmers e-Commons e-Infrastructure group at Chalmers University of Technology, Gothenburg, Sweden.

Appendix A Photometric fluxes and best-fit SEDs

The photometric flux densities (see Sect. 3.1) for the five stars used for SED fitting are listed in Table A.1. The dust radiative transfer modelling procedure adopted is explained in Sect. 4.2, and the results are presented in Sect. 5.1. See Table 3 for the stellar and dust parameters calculated from the SED fitting. The best-fit SEDs obtained from the modelling are shown below in Fig. A.1.

Fig. A.1 Photometric data (black, Sect. 3.1) and modelled SEDs (red, Sect. 4.2).

Appendix B CO line intensities and best-fit models

This appendix presents the velocity-integrated CO line intensities used to constrain our gas radiative transfer models. The line intensities from the SD observations (Sect. 3.2) used are listed in Table B.1 for all five sources, and the integrated intensities for the spectra from different apertures extracted from the ACA lines (Sect. 3.2) for IRAS 07454-7112 and AFGL 3068 are given in Table B.2. The CO modelling is described in Sect. 4.3, and the corresponding results are presented in Sect. 5.2.

Fig. B.1 Observed (black, Sect. 3.2) and modelled (red, Sect. 4.3) CO line profiles for IRAS 15194-5115. The transition quantum numbers, telescope used for the observation, and the corresponding beam size in arcseconds are listed at the top right corner of each panel. The lines are in the main-beam temperature (TMB [K]) scale for all SD spectra shown, whereas the flux is given in units of Jansky for the ALMA lines where available.

Fig. B.2 Same as Fig. B.1, but for IRAS 15082-4808.

Fig. B.3 Same as Fig. B.1, but for IRAS 07454-7112.

Fig. B.4 Same as Fig. B.1, but for IRC +10 216.

Fig. B.5 Same as Fig. B.1, but for AFGL 3068.

Appendix C CS line fits and χ² maps

This appendix shows the observed and modelled best-fit CS line profiles for the stars in the sample other than IRAS 15194-5115 (Figs. C.1–C.3). The best-fit model for IRAS 15194-5115 is shown in Fig. 1. The χ² maps for the grid of models we ran for each of the five sources, with the CS peak abundance (f₀) and the e-folding radii (R_e) of the CS abundance profile as free parameters are also shown (Figs. C.4, C.5). The modelling is described in detail in Sect. 4.4 and the results obtained are presented in Sect. 5.3.

Fig. C.1 Observed (black) and modelled (red) CS line profiles for (a) IRAS 15082-4808 and (b) IRAS 07454-7112. The transition quantum numbers, telescope used for the observation, and the corresponding beam size in arcseconds are listed at the top right corner of each panel. The line fluxes are in units of Jansky for the ALMA spectra, whereas they are given in the main-beam temperature (TMB [K]) scale for all SD spectra shown.

Fig. C.2 Same as Fig. C.1, but for AFGL 3068.

Fig. C.3 Same as Fig. C.1, but for IRC +10 216.

Fig. C.4 Same as Fig. 3, but for (a) IRAS 15194-5115, (b) IRAS 15082-4808, and (c) IRAS 07454-7112. The figures show χ2 values for our grid of CS models for the above sources. The contours denote 1σ, 2σ, and 3σ ranges, when constrained using only the spatial information from the CS J = 2 - 1 line (green) vs just the excitation information from the different J lines (blue). Also shown are the contours when using both the spatial and excitation information simultaneously (red) to constrain the RT models. The green, blue and red stars mark the best-fit models as estimated from the above three cases, respectively.

Fig. C.5 χ2 maps for the CS RT models. The free parameters f0 and Re are on the x- and y-axes, respectively. The red + sign denotes the best-fit model, and the white contours correspond to the 1σ, 2σ, and 3σ ranges.

Appendix D ¹³CS and C³⁴S line profile fits and intensity profiles

Tables D.1 and D.2 list the observed integrated intensities of the ¹³CS and C³⁴S lines modelled in this work respectively. The line profile fits for ¹³CS and C³⁴S are also shown in this appendix. The lines of these molecules are modelled by only varying the peak abundance, keeping the e-folding radius fixed to the value obtained from CS modelling for each source. See Sect. 4.4.2 for details.

Fig. D.1 Comparison of modelled (red) and observed (black) integrated intensities of 13CS J = 2 - 1 spectra extracted from successively larger beams (see Sect. 4.4.2). The blue-shaded region represents the uncertainty in the observed line intensities.

Fig. D.2 Observed (black, Sect. 3.3) and modelled (red, Sect. 4.4) 13CS line profiles for IRAS 15194-5115. The transition quantum numbers, telescope used for the observation, and the corresponding beam size in arcseconds are listed at the top right corner of each panel. The line fluxes are in units of Jansky for the ALMA spectra, and in main-beam temperature (TMB [K]) scale for all SD spectra.

Fig. D.3 Same as Fig. D.2, but for IRAS 15082-4808.

Fig. D.4 Same as Fig. D.2, but for IRAS 07454-7112.

Fig. D.5 Same as Fig. D.2, but for AFGL 3068 (left) and IRC +10 216 (right).

Fig. D.6 Same as Fig. D.1, but for C34S.

Fig. D.7 Observed (black, Sect. 3.3) and modelled (red, Sect. 4.4) C34S line profiles for IRAS 15194-5115. The transition quantum numbers, telescope used for the observation, and the corresponding beam size in arcseconds are listed at the top right corner of each panel. The line fluxes are in units of Jansky for the ALMA spectra, whereas they are given in the main-beam temperature (TMB [K]) scale for all SD spectra shown.

Fig. D.8 Same as Fig. D.7, but for IRAS 15082-4808.

Fig. D.9 Same as Fig. D.7, but for IRAS 07454-7112.

Fig. D.10 Same as Fig. D.7, but for AFGL 3068.

Fig. D.11 Same as Fig. D.7, but for IRC +10 216.

References

All Tables

All Figures

	Fig. 1 Observed (black) and modelled (red) CS line profiles for IRAS 15194-5115. The transition quantum numbers, telescope used for the observation, and the corresponding beam size in arcseconds are listed at the top right corner of each panel. The line fluxes are in units of Jansky for the ALMA spectra, whereas they are given in the main-beam temperature (TMB [K]) scale for all SD spectra shown.
In the text

	Fig. 2 Comparison of modelled (red) and observed (black) integrated intensities of CS J = 2-1 spectra extracted from successively larger beams (see Sect. 4.4.2). The blue-shaded region represents the uncertainty in the observed line intensities.
In the text

Fig. 3 χ2 values for our grid of CS models for IRC +10 216. The contours denote 1σ, 2σ, and 3σ ranges, when constrained using only the spatial information from the CS J = 2-1 line (green) vs just the excitation information from the different J lines (blue). Also shown are the contours when using both the spatial and excitation information simultaneously (red) to constrain the RT models. The green, blue, and red stars mark the best-fit models estimated from the above three cases, respectively.

In the text

	Fig. 4 Isotopic ratios from this work compared to estimates from the literature.
In the text

	Fig. 5 Ratio of the excitation temperature (Tex) and the kinetic temperature (Tkin) for various low- and high-J CS transitions (top), and the corresponding line emitting regions (bottom), normalised to the peak of the J = 21-20 emitting region, as functions of radius across the CSE, from the best-fit model for IRAS 15194-5115. See text for a more detailed description.
In the text

Fig. 6 Variation in e-folding radius of the CS abundance profiles with circumstellar gas density. The dashed cyan line (Maes et al. 2023) is from chemical modelling results, while the other trends plotted are calculated from RT models. The magenta point is obtained from chemical modelling results (Van de Sande et al., in prep.), for a star of MLR 1.0×10-5 Mʘ yr−1 and a ν∞ of 15 km s−1, and the corresponding uncertainty range is estimated from sensitivity calculations on chemical models (see text). The rightmost blue circle is for AFGL 3068, and is likely underestimated (see text).

In the text

Fig. 7 Comparison of CS abundance profiles obtained from RT modelling (blue) and chemical modelling: red - using CS photodissociation rate from Pattillo et al. (2018), green: using CS photodissociation rate calculated from the cross sections from Hrodmarsson & van Dishoeck (2023), based on Xu et al. (2019). The shaded blue region represents the 1σ uncertainty range in the RT model abundance profiles. The circles denote the Re of the respective models. The red and green error bars on the chemical model Re denote the estimated uncertainty in the chemical models (Van de Sande et al., in prep.). It is assumed that this uncertainty does not vary significantly with the outflow density and photodissociation rate (see text).

In the text

	Fig. A.1 Photometric data (black, Sect. 3.1) and modelled SEDs (red, Sect. 4.2).
In the text

Fig. B.1 Observed (black, Sect. 3.2) and modelled (red, Sect. 4.3) CO line profiles for IRAS 15194-5115. The transition quantum numbers, telescope used for the observation, and the corresponding beam size in arcseconds are listed at the top right corner of each panel. The lines are in the main-beam temperature (TMB [K]) scale for all SD spectra shown, whereas the flux is given in units of Jansky for the ALMA lines where available.

In the text

	Fig. B.2 Same as Fig. B.1, but for IRAS 15082-4808.
In the text

	Fig. B.3 Same as Fig. B.1, but for IRAS 07454-7112.
In the text

	Fig. B.4 Same as Fig. B.1, but for IRC +10 216.
In the text

	Fig. B.5 Same as Fig. B.1, but for AFGL 3068.
In the text

Fig. C.1 Observed (black) and modelled (red) CS line profiles for (a) IRAS 15082-4808 and (b) IRAS 07454-7112. The transition quantum numbers, telescope used for the observation, and the corresponding beam size in arcseconds are listed at the top right corner of each panel. The line fluxes are in units of Jansky for the ALMA spectra, whereas they are given in the main-beam temperature (TMB [K]) scale for all SD spectra shown.

In the text

	Fig. C.2 Same as Fig. C.1, but for AFGL 3068.
In the text

	Fig. C.3 Same as Fig. C.1, but for IRC +10 216.
In the text

Fig. C.4 Same as Fig. 3, but for (a) IRAS 15194-5115, (b) IRAS 15082-4808, and (c) IRAS 07454-7112. The figures show χ2 values for our grid of CS models for the above sources. The contours denote 1σ, 2σ, and 3σ ranges, when constrained using only the spatial information from the CS J = 2 - 1 line (green) vs just the excitation information from the different J lines (blue). Also shown are the contours when using both the spatial and excitation information simultaneously (red) to constrain the RT models. The green, blue and red stars mark the best-fit models as estimated from the above three cases, respectively.

In the text

	Fig. C.5 χ2 maps for the CS RT models. The free parameters f0 and Re are on the x- and y-axes, respectively. The red + sign denotes the best-fit model, and the white contours correspond to the 1σ, 2σ, and 3σ ranges.
In the text

	Fig. D.1 Comparison of modelled (red) and observed (black) integrated intensities of 13CS J = 2 - 1 spectra extracted from successively larger beams (see Sect. 4.4.2). The blue-shaded region represents the uncertainty in the observed line intensities.
In the text

	Fig. D.2 Observed (black, Sect. 3.3) and modelled (red, Sect. 4.4) 13CS line profiles for IRAS 15194-5115. The transition quantum numbers, telescope used for the observation, and the corresponding beam size in arcseconds are listed at the top right corner of each panel. The line fluxes are in units of Jansky for the ALMA spectra, and in main-beam temperature (TMB [K]) scale for all SD spectra.
In the text

	Fig. D.3 Same as Fig. D.2, but for IRAS 15082-4808.
In the text

	Fig. D.4 Same as Fig. D.2, but for IRAS 07454-7112.
In the text

	Fig. D.5 Same as Fig. D.2, but for AFGL 3068 (left) and IRC +10 216 (right).
In the text

	Fig. D.6 Same as Fig. D.1, but for C34S.
In the text

Fig. D.7 Observed (black, Sect. 3.3) and modelled (red, Sect. 4.4) C34S line profiles for IRAS 15194-5115. The transition quantum numbers, telescope used for the observation, and the corresponding beam size in arcseconds are listed at the top right corner of each panel. The line fluxes are in units of Jansky for the ALMA spectra, whereas they are given in the main-beam temperature (TMB [K]) scale for all SD spectra shown.

In the text

	Fig. D.8 Same as Fig. D.7, but for IRAS 15082-4808.
In the text

	Fig. D.9 Same as Fig. D.7, but for IRAS 07454-7112.
In the text

	Fig. D.10 Same as Fig. D.7, but for AFGL 3068.
In the text

	Fig. D.11 Same as Fig. D.7, but for IRC +10 216.
In the text