© The Authors 2023

1. Introduction

The asymptotic giant branch (AGB) is a late stage of the stellar evolution of stars with masses between ∼0.8 M_⊙ and 8 M_⊙. Stars in this phase are characterised by their high mass-loss rates of 10⁻⁸ M_⊙ yr⁻¹ to 10⁻⁴ M_⊙ yr⁻¹. As a result of this mass loss, AGB stars have a circumstellar envelope (CSE) that is formed from the escaping material and which is cool and dense enough for several types of molecules and dust to efficiently form. There are three chemical types of CSEs, depending on their carbon-to-oxygen ratio: carbon-rich (C/O > 1), oxygen-rich (C/O < 1), and S-type (C/O ∼ 1).

Sodium chloride (NaCl) has been detected towards several evolved stars, including the C-rich AGB star CW Leonis (CW Leo; Cernicharo & Guelin 1987; Agúndez et al. 2012; Quintana-Lacaci et al. 2016), the O-rich AGB star IK Tau (Milam et al. 2007; Velilla Prieto et al. 2017; Decin et al. 2018), and the red supergiants (RSGs) VY Canis Majoris (VY CMa; Milam et al. 2007; Decin et al. 2016) and NML Cyg (Singh et al. 2022). It is a diatomic molecule with a high dipole moment (μ = 9.0 D; De Leeuw & Dymanus 1970), which facilitates its detection because it allows for efficient radiative excitation and produces bright lines. Like most metal-bearing molecules, NaCl is highly refractory (Decin et al. 2016), and it is expected to easily condense onto dust grains (Höfner 2008; Höfner et al. 2016).

Cernicharo & Guelin (1987) first detected Na³⁵Cl (hereafter NaCl) and Na³⁷Cl towards an AGB star, the carbon star CW Leo. Based on observations from the Institut de Radio Astronomie Millimétrique (IRAM) 30 m telescope, Agúndez et al. (2012) modelled the NaCl abundance around CW Leo, assuming a 1D number density profile, and found a total NaCl abundance (both ³⁵Cl and ³⁷Cl) of 1.8 × 10⁻⁹ between 3R_⋆ and ∼100R_⋆ before a drop-off, suggesting that NaCl exists in both the hot inner layers and in the cooler regions of the CSE of CW Leo. They also found an isotopic ratio of Na³⁵Cl/Na³⁷Cl = 2.9 ± 0.3, which is smaller than but close to the solar ratio of 3.127 ± 0.006 (Rosman & Taylor 1998; Asplund et al. 2021). This was confirmed by the results of Quintana-Lacaci et al. (2016), who used the Atacama Large Millimetre/sub-millimetre Array (ALMA) to observe spatially resolved metal-bearing molecules in the inner regions of the CSE of CW Leo. They detected NaCl up to 83R_⋆ from the star and determined its distribution to be either a spiral or a torus. Decin et al. (2016) used ALMA to observe spatially resolved NaCl in the CSE of the RSG VY CMa and found it to be present out to a distance of 220R_⋆, which is beyond the dust condensation radius, despite the highly refractory nature of the molecule. These results suggest that there is a chemical process preventing all NaCl from condensing onto dust grains, such as thermal desorption or shock-induced sputtering. As with CW Leo, Decin et al. (2016) detected non-spherically symmetric structures in the NaCl emission towards VY CMa, but unlike the case of CW Leo, these take the form of separated clumps.

IK Tau is an O-rich, Mira variable AGB star with a moderately high mass-loss rate of ∼5 × 10⁻⁶ M_⊙ yr⁻¹ (Decin et al. 2018). Decin et al. (2018) created a spectral atlas of all species detected in an ALMA spectral scan between 335 and 362 GHz, including the relevant parameters of the detected features. They found a complex wind morphology and both large- and small-scale inhomogeneities in the form of clumps and/or arcs, which were observed in the emission of NaCl, CO, HCN, SiS, and CS in the inner wind.

Thus far, NaCl towards IK Tau has been studied using low spatial resolution and unresolved observations. Milam et al. (2007) modelled the abundance of NaCl in the oxygen-rich CSEs of IK Tau and the RSG VY CMa using observations at 1 mm and 2 mm of the Arizona Radio Observatory (ARO) Submillimeter Telescope (SMT) and the ARO 12m telescope, respectively. They found an abundance relative to H₂ of ∼4 × 10⁻⁹ for IK Tau, assuming a spherical distribution of NaCl, a low expansion velocity of 3 km s⁻¹, and a source size of 0.3″. Velilla Prieto et al. (2017) observed NaCl towards IK Tau with the IRAM 30m telescope and used population diagrams to determine its abundance relative to H₂, finding a value of 3.1 × 10⁻⁷, around two orders of magnitude higher. They used a 1D radiative transfer model and assumed local thermodynamic equilibrium (LTE), an optically thin CSE, and a source size of 0.3″. Decin et al. (2018) detected NaCl and Na³⁷Cl towards IK Tau with ALMA, but not towards R Dor, an O-rich AGB star with a low mass-loss rate. The observed NaCl is not centred on the continuum peak but is observed in irregular clumps around the stellar position; none of the other molecules they observed exhibited this behaviour. Gobrecht et al. (2016) provide a theoretical prediction of the abundance of NaCl around IK Tau that they determined by modelling the dust formation in the inner wind of IK Tau using a 1D, shock-induced chemical model. They predict an abundance of NaCl between 3.7 × 10⁻¹⁰ and 1 × 10⁻⁸ relative to H₂.

In this paper we aim to study the 3D distribution of NaCl around IK Tau to determine the abundance of NaCl relative to H₂ for each individual clump. Section 2 presents the observations. In Sect. 3.2 we de-project the observed channel maps to obtain the 3D distribution of NaCl. In Sects. 3.3 and 3.4 we use the radiative transfer modelling code Magritte to estimate the abundance of NaCl relative to H₂ for each of the observed clumps. In Sect. 4 these results are discussed, and in Sect. 5 we give our conclusions.

2. Observations

IK Tau was observed by ALMA in Band 7 in August 2015 (proposal 2013.1.00166.S, PI L. Decin). The measurements consist of an unbroken spectral scan between 335–362 GHz, with a spatial resolution of ∼120 × 150 mas and a frequency resolution of 1.95 MHz. The data reduction is outlined in Decin et al. (2018) and was not repeated for the present work. Two NaCl and two Na³⁷Cl emission lines were covered in the observed frequency range. In Table 1 we give the quantum numbers and rest frequencies of the covered NaCl lines, along with the velocity resolution of the observations and the σ_rms. The lines correspond to transitions between rotational energy levels of the ground vibrational state.

The spectra of these lines, extracted for a circular aperture with radius 320 mas, are shown in Fig. 1. While NaCl (26 → 25) and (27 → 26) are clearly detected above the noise, Na³⁷Cl (27 → 26) and (28 → 27) are only tentatively detected. This is partly owing to the fact that both Na³⁷Cl lines lie in noisier parts of the observed spectrum (see Fig. 3 in Decin et al. 2018). There were some blends with NaCl lines reported in Decin et al. (2018), which we find are not significant for the following reasons. NaCl (27 → 26) is reported to be blended with Si¹⁸O (J = 9 → 8, v = 5) at 350.959 GHz; however, this identification is tentative. With the covered Si¹⁶O (J = 8 → 7, v = 5) line only having a peak flux of 0.059 Jy and ¹⁶O/¹⁸O > 400 for IK Tau (Danilovich et al. 2017a), the contribution from the Si¹⁸O line to the total flux should be negligible. The only other possible blend is between Na³⁷Cl (27 → 26) and SO₂ (57_15, 43 → 58_14, 44) at 343.477 GHz, but a detailed study of the SO₂ emission in this dataset reports that this line was not detected towards IK Tau (Danilovich et al. 2020b).

Fig. 1. Four NaCl lines covered in the spectral scan of IK Tau by Decin et al. (2018). The lines were extracted using a circular aperture with a radius of 320 mas, centred on the star. Spectra for additional extraction apertures are shown in Fig. A.1.

To more carefully examine the Na³⁷Cl emission, we plotted all four NaCl line profiles as extracted from different regions of the channel maps. As well as the circular extraction aperture with a radius of 320 mas that we used for Fig. 1, we also used a circular aperture with radius 800 mas, which results in noisier spectra but captures all of the NaCl flux, some of which was outside of the smaller aperture. We also extract spectra from two irregular apertures centred on the brightest two regions of NaCl flux, which we refer to as clumps A and C (see below for the complete description of clump labels). As can be seen in Fig. A.1, clump A is characteristically blue-shifted while clump C is characteristically red-shifted. Although the Na³⁷Cl lines are still not clearly detected above the noise, the clump A spectra tend to show more blue emission, while the clump C spectra tend to show more red emission. This suggests that Na³⁷Cl behaves in broadly the same way as Na³⁵Cl. We also note that the 800 mas NaCl line profiles appear to be narrower than the 320 mas profiles. This is because clump A dominates the 800 mas line profiles while being partly cropped out of the 320 mas profiles, and because the less intense redder emission is partly obscured by the noise in the 800 mas profiles. The line widths of the more extensive NaCl transitions covered by Velilla Prieto et al. (2017) (with levels from J = 7 to J = 24) with the IRAM 30m antenna are generally reported to be narrower than other molecular lines such as CO, SiS, and SiO. Although those observations are not spatially resolved, we suggest that the narrower line profiles are a result of clump A dominating the NaCl emission generally. By comparing the integrated spectra of the (27 → 26) line for both Na³⁵Cl and Na³⁷Cl, we found a lower limit for the ratio Na³⁵Cl/Na³⁷Cl > 2.1, taking the different line frequencies into account as described by Danilovich et al. (2020b). This limit is comparable to the ratio found for the carbon star CW Leo (Kahane et al. 2000; Agúndez et al. 2012). More sensitive observations would allow this ratio to be constrained with more precision.

Figure 2 shows the channel maps of the NaCl J = 26 → 25 transition, and the channel maps of transitions NaCl J = 27 → 26 can be found in Fig. A.2. The location of the star on the channel maps is indicated by a red star, and is assumed to be at the location of the continuum peak. The systemic velocity of the star, υ_LSR = 34 km s⁻¹ (Maercker et al. 2016), has been subtracted from the channel map velocities. The emission in the Na³⁷Cl channel maps is generally not distinguishable from the noise. Hence, we only show the spectra for these tentative lines.

Fig. 2. Channel maps of the NaCl J = 26 → 25 line. The red star indicates IK Tau’s position, and the white ellipse shows the beam size. The white contours are at 3σrms, 5σrms, and 10σrms (σrms = 2.4 mJy). The υLSR = 34 km s−1 has been subtracted from the velocities.

As reported by Decin et al. (2018), the NaCl emission is not centred on the star, nor is it spherically symmetric. Instead it is unevenly distributed in clumps around the star, with no clumps present to the west of the star. We identify up to five clumps comprising the NaCl emission. In Fig. 2 we plot the NaCl J = 26 → 25 channel maps, and label five different clumps with the letters A to E. However, only three of these clumps (A, B and C) are detected in the NaCl J = 27 → 26 channel maps (Fig. A.2). The lower signal-to-noise of the NaCl J = 27 → 26 channel maps can be explained by the slightly higher noise at those frequencies and the slightly weaker signal (e.g., compare the spectra in Fig. A.1, where the peak flux of NaCl J = 26 → 25 is generally 5–20% brighter than the peak flux of NaCl J = 27 → 26), possibly owing to the excitation conditions. We also note that the spectrum of clump C is the only one where NaCl J = 27 → 26 is 5% brighter than NaCl J = 26 → 25 (which can also be seen by comparing Figs. 2 and A.2), indicating that the excitation conditions vary between clumps. Clumps are considered in our analysis if their 3σ_rms contours are larger than or roughly equal in size to the synthetic beam, or are in a channel map and spatial position directly adjacent to such a clump. It should also be noted that clumps A and B are counted as different clumps, despite them being connected to each other, because the connection between them has a smaller cross-section than the main regions of those clumps, seeming to form a bridge between the two clumps. Additionally, clump E is only slightly larger than the synthetic beam, making it an uncertain detection, especially since it is not detected in the NaCl J = 27 → 26 channel maps. The existence of clump D is less straightforward to determine. It is only detected in one channel of NaCl J = 26 → 25 (at −0.5 km s⁻¹; Fig. 2) but that detection is above 5σ. The 3σ contour is close to, but slightly larger than, the beam size, meaning that the clump emission spans several pixels within the −0.5 km s⁻¹ channel. The non-detection of clump D in NaCl J = 27 → 26 could be because the velocity of the clump is very localised in the line-of-sight direction. The somewhat coarse velocity resolution of our data (1.7 km s⁻¹) could, for example, result in clump D being detected in NaCl J = 26 → 25 if it is (close to) centred on the channel velocity of −0.5 km s⁻¹, but not detected in NaCl J = 27 → 26 where its flux it might be split between the −1.0 km s⁻¹ and 0.6 km s⁻¹ channels.

Overall, the NaCl J = 26 → 25 channel maps have clumps detected for local standard of rest (LSR) velocities between υ = −16.1 km s⁻¹ and υ = 18.5 km s⁻¹, while the NaCl J = 27 → 26 channel maps have clumps detected between υ = −14.4 km s⁻¹ and υ = 14.0 km s⁻¹. Table 2 shows at which velocities each labelled clump is detected. There is mostly agreement between the two transitions, although it should be noted that sometimes clumps are not detected above the noise in the NaCl J = 27 → 26 channel maps at velocities for which they are detected in the NaCl J = 26 → 25 channel maps, which, again, can be explained by the difference in signal-to-noise at different frequencies and the slightly weaker apparent flux of NaCl J = 27 → 26.

While the focus of the present work is on NaCl, some other molecular lines were also considered to aid our analysis. These are the SiS J = 19→ 18 (v = 0) and CS J = 7 → 6 (v = 0) lines (both analysed in detail by Danilovich et al. 2019) and the CO J = 3 → 2 (v = 0) and HCN J = 4 → 3 (v = 0) lines. The spectra of these lines are plotted in Fig. B.1 and are used for the velocity determination discussed in Sect. 3.1.

3. 3D modelling

In this section we describe how the observed channel maps can be turned into 3D models of the NaCl abundance by de-projection. This de-projection is then used with a radiative transfer modelling code to determine the abundance of NaCl in the clumps by creating synthetic channel maps and comparing these to the observed ones. To do this, first a new maximum expansion velocity is determined, as ALMA’s increased sensitivity allows for the detection of high velocity wings that were not taken into account when the value from the literature was determined. Table 3 summarises the stellar parameters of IK Tau that are used in this section.

3.1. Expansion velocity determination

Prior to the higher sensitivity of ALMA observations, the terminal expansion velocity of IK Tau was calculated to be υ_exp = 17.5 km s⁻¹ (e.g., Maercker et al. 2016, from APEX and Herschel/HIFI observations). Decin et al. (2018) found, by fitting a velocity profile to the measured widths of various molecular lines at zero intensity, a maximal expansion velocity of ∼25 km s⁻¹ for IK Tau, owing to high-velocity wings that were not detected above the noise in earlier single-antenna observations (Kim et al. 2010; Decin et al. 2010a,b; Maercker et al. 2016; Velilla Prieto et al. 2017). This leaves some ambiguity as to the true terminal expansion velocity of IK Tau¹; the bulk of the circumstellar material appears to have a maximum expansion velocity of 17.5 km s⁻¹, and the line profiles without wings are reproduced by radiative transfer models with this expansion velocity (Maercker et al. 2016; Danilovich et al. 2019), but some material expands at a higher velocity. Furthermore, the underlying gas velocity field should be independent of molecular probes; in other words, although the NaCl emission is contained within 17.5 km s⁻¹ of the LSR velocity, we should not ignore the fact that emission from some other lines lies outside this velocity range. For example, in Fig. B.1 we plot the SiS J = 19 → 18 (v = 0), CO J = 3 → 2 (v = 0), CS J = 7 → 6 (v = 0), and HCN J = 4 → 3 (v = 0) spectral lines, where the wings of the lines extend beyond υ_exp = 17.5 km s⁻¹.

Our de-projection technique (see Sect. 3.2) assumes a monotonically increasing radial velocity profile. For consistency with the observed high velocity wings, we derive a new maximum expansion velocity by carefully examining the SiS J = 19 → 18 and CS J = 7 → 6 lines, which do not show absorption in their wings as the CO J = 3 → 2 and HCN J = 4 → 3 lines do (Fig. B.1). These lines were also chosen because their emission is thought to be approximately spherical, based on channel maps (Decin et al. 2018; Danilovich et al. 2019), and because they are relatively abundant species that are not expected to have significant maser activity. The de-projection was done by finding the velocity corresponding to the red and blue line wings and taking the one with the largest magnitude (method outlined in Appendix A of Gottlieb et al. 2022), as illustrated in Fig. 3 for the SiS line. There, the horizontal dotted orange line at 0.047 Jy beam⁻¹ corresponds to 3σ peak rms noise. The velocities of the red and blue line wings are then determined by finding where the spectral line first goes below this criterion of 3σ_rms. The red wing velocity corresponds to 26.0 km s⁻¹, while the blue wing velocity corresponds to −28.4 km s⁻¹. The maximum expansion velocity is then the wing velocity with the greatest magnitude, which is the blue wing velocity, υ_exp = 28.4 ± 1.7 km s⁻¹. We used the velocity width of a channel map as the uncertainty. Doing the same calculation for CS, we find a lower expansion velocity, υ_exp = 23.6 ± 1.7, probably because CS has a lower signal-to-noise ratio than SiS. Hence, we used the value obtained from SiS in our calculations.

Fig. 3. Line profiles of SiS. Left: maximum expansion velocity determination using the SiS 19 → 18 (v = 0). The horizontal dotted orange line is the used criterion, 3σrms, and the vertical dashed lines indicate the red and blue wing velocities. Right: same as on the left, but zoomed in to better discern the red and blue wings. The method used is from Gottlieb et al. (2022).

3.2. De-projection

To model the observed NaCl clumps towards IK Tau, we de-project the channel maps, shown in Fig. 2, to transform the observed (RA, Dec, υ) coordinates to spatial (x, y, z) coordinates relative to the stellar position. Here, x and y are related to the coordinates in the plane of the sky, the right ascension and declination, respectively, while the z coordinate is the spatial coordinate along the line of sight. These new coordinates (x, y, z) all share the same units of distance, while the old coordinates have differing units (RA and Dec are angles and υ has velocity units). The spatial (x, y) coordinates of a point in a channel map can be found using the right ascension and declination offsets (Δα, Δδ) and the distance, D, to the star (see Table 3):

(1)

(2)

Here, the right ascension and declination offsets (Δα, Δδ) are defined relative to the position of IK Tau, such that the star’s coordinates on the channel maps are (Δα, Δδ) = (0, 0) mas and (x, y) = (0, 0) AU. An equation for the z coordinate can be calculated based on the Doppler shift of the emission, assuming a monotonic velocity field along the line of sight (see e.g., Guélin et al. 2018; Montargès et al. 2019), so that

(3)

A full derivation of this equation can be found in Appendix C.1. Here, υ_channel is the velocity of a given channel map, and υ_exp is the constant, non-zero radial velocity of the stellar wind, taken to be the maximum expansion velocity. However, previous studies of IK Tau have used a β-type velocity law, which is a common 1D model for the accelerating gas velocity of the bulk material of the stellar wind of an AGB star (see for example Maercker et al. 2016):

(4)

Here υ₀ = 3 km s⁻¹ is the inner wind velocity (Maercker et al. 2016), generally assumed to be the sound speed; υ_exp = 28.4 ± 1.7 km s⁻¹ is the maximum expansion velocity that was derived in Sect. 3.1; r_dust = 2.38 × 10¹⁴ cm is the dust condensation radius (Decin et al. 2018); and β = 1.5 is a parameter determining how fast the wind accelerates (Maercker et al. 2016). The values of these parameters are also given in Table 3.

Since the β-type velocity law is monotonic along the line of sight, we can replace the constant wind velocity υ_exp in Eq. (3) with the velocity profile Eq. (4), and get the following equation:

(5)

which we used to de-project the channel maps. Here, υ_channel is the central velocity of a channel map. It should be noted that since the velocity profile given in Eq. (4) is dependent on the radial distance , the equation is implicit and needs to be solved numerically. For this, the scipy (Virtanen et al. 2020) implementation of Brent’s method (Brent 1973) was used (scipy.optimize.brentq). On a single channel map, the points with a small impact parameter, (i.e., those that lie close to the star along the line of sight), will be de-projected to smaller z values than the points with a larger impact parameter. This effect is enhanced for channel maps with velocities that are much greater than zero. The z-coordinate thus increases in size when moving away from (x,y) = (0,0) AU and/or υ_channel = 0 km s⁻¹. This is illustrated in Fig. C.2.

Before the NaCl J = 26 → 25 channel maps are de-projected, a few practical adjustments need to be made. First, only the channel maps with an absolute velocity lower than the maximum expansion velocity are de-projected. No NaCl emission is detected outside of these channel maps. Then the pixels that lie within the 3σ contours that are either larger than the size of the synthetic beam or are in a channel map and spatial position adjacent to such a 3σ contour are identified. The intensities of the pixels that do not lie within these contours are put to zero to avoid complicating the calculation with noise-dominated regions. After this, the coordinate transformation from (Δα, Δδ) to (x, y) (Eqs. (1) and (2)) is done, and the channel maps were cut down to a smaller size in x, y, encompassing all the clumps.

Next, the intensity is interpolated to a grid that is six times as dense in the x and y directions, and 12 times as dense in the υ direction, yielding a grid of 132 × 156 × 300 points. The higher grid-density in the υ coordinate is chosen because there are large spacings in the z direction of the (x, y, z) grid corresponding to the given (x, y, υ) grid. As mentioned before, the z coordinate will have a larger (smaller) absolute value when the x, y, and/or υ coordinates have absolute values far from (close to) 0. This leads to a non-constant spacing in the z direction, which can be large for grid points far from the origin. The higher grid-density reduces the distance between grid points in the z-direction. This grid is then de-projected using Eq. (5). It is later used as an input for the radiative transfer modelling code Magritte (see Sect. 3.3), and the higher grid density allows for more granular mapping of quantities such as molecular number density and temperature.

Figure 4 shows the de-projected clumps. The colours relate to the observed intensity at those points, and the little yellow dot shows the central star, IK Tau, to scale (Decin et al. 2018). The observed structures extend about ∼200 AU in the x-direction, ∼240 AU in the y direction and ∼140 AU in the z-direction, meaning that the structures are flatter in the z-direction, given the assumed υ_β law.

Fig. 4. De-projection of the NaCl J = 26 → 25 channel maps shown from different angles. Only points detected at ≥3σ are shown, and the colours indicate intensity, using the same colour bar as for the original channel maps (Fig. 2). The yellow arrow points in the direction of the line of sight.

3.3. Radiative transfer modelling

The channel maps (Figs. 2 and A.2) and our de-projection show that the distribution of NaCl is very asymmetric, so a 1D radiative transfer model cannot be used to study the emission. Instead, the 3D radiative transfer code Magritte² (De Ceuster et al. 2019, 2020) is used. Magritte computes the radiation field by solving the radiative transfer equation along rays (straight lines) through the medium. For numerical stability, Magritte solves the radiative transfer equation in a second-order, Schuster-Feaurier form (Schuster 1905); feautrier. More technical details are available in De Ceuster et al. (2019).

The input for Magritte consists of (1) a 3D point cloud with (a) spatial coordinates, (b) a velocity field, (c) a temperature profile, (d) the number densities of all considered molecular species, and (e) a turbulent velocity field; (2) the boundary of the grid; (3) the transition rates of the considered molecular species ; and (4) the frequencies at which it needs to compute the output. The de-projected grid described in Sect 3.2 is used as the 3D grid with spatial coordinates, after it has been interpolated a second time to a new grid of 132 × 156 × 115 points. This grid gives the z direction the same spatial scale as the x and y directions, such that the grid is evenly spaced, and crops points representing empty space to make the model more efficient. To reduce computation times further, the Magritte model is re-meshed (De Ceuster et al. 2020), meaning that the grid will be resampled to have more points in the high density regions and have fewer points in the low density regions. The reduced Magritte model has a spherical inner boundary centred on the star, with a radius equal to the stellar radius (∼2.54 AU), and a cuboidal outer boundary.

The radiative transfer code requires the three vector components of the velocity profile rather than the total radial velocity (calculated using Eq. (4)). Assuming a radially outward spherically symmetric outflow, the components are given by

(6)

with υ_i the component of the velocity along the x, y or z direction, υ the radial velocity at point (x, y, z), and .

We used a temperature profile based on the results of Maercker et al. (2016). The inner boundary of this temperature profile is the dust condensation radius R = 18 × 10¹³ cm, which is not small enough for our purposes. Therefore, we linearly extrapolate it to match the effective temperature of the star at r = R_⋆. The extrapolated temperature profile can be seen in the top panel of Fig. 5.

Fig. 5. Physical conditions in the CSE of IK Tau. Top: temperature profile derived by Maercker et al. (2016), extrapolated to the stellar radius and effective temperature of IK Tau. Bottom: number density profile of H2 calculated from the mass-loss rate using Eq. (7).

In the model, we consider H₂ and NaCl, where H₂ serves as an approximation of the total circumstellar gas. Since the H₂ abundance is typically much greater than that of any other molecular species, its number density can be approximated from the mass-loss rate of the star as follows, assuming a spherically symmetric outflow and that hydrogen is mainly in its molecular form H₂:

(7)

where υ_β(r) is the gas expansion rate at radial distance r, which is calculated using the beta-type velocity law Eq. (4). The radial dependence of this formula is plotted in the bottom panel of Fig. 5.

We assume LTE in our radiative transfer modelling and we include only the ground vibrational levels. The modelling time would increase by about four orders of magnitude without the assumption of LTE, with an additional multiplier for higher vibrational states. We expect that models using non-LTE would provide more accurate results, especially farther from the star, where the density is low. Infrared radiative pumping might have a significant contribution, but this cannot be taken into account for the model without including more vibrationally excited levels. The radiative transition rates and energy levels for NaCl are taken from the ExoMol line list database (Barton et al. 2014), and the collisional rates for NaCl were calculated using the formulation provided by Quintana-Lacaci et al. (2016). We consider 40 radiative transitions between 41 energy levels of NaCl.

By construction, the observations are binned in frequency channels yielding one channel map per frequency bin. As a result, any variation of the intensity within a frequency bin is integrated out, and we only observe the integrated intensity over the frequency bin. The radiative transfer solver, however, constructs the synthetic channel maps by computing the intensity at specific frequencies. Therefore, any variation of the intensity in between those frequencies is not necessarily integrated out, and we obtain the intensity at a certain frequency rather than the integrated one over a frequency bin. To simulate the integration over the frequency bins, we computed intensity maps for three equally spaced frequencies within each bin, and summed the intensity maps. An example of these channel maps can be found in Fig. A.3.

The resulting synthetic channel maps do not yet take the ALMA beam size into account. This can be resolved by using the Common Astronomy Software Applications (CASA, CASA Team 2022), which is a data processing software for radio interferometers, including ALMA. By using the simalma tool, it is possible to simulate what the synthetic channel maps from Magritte would look like if they were actually observed by ALMA. IK Tau was observed when the precipitable water vapour (pwv) was in the range [0.2, 1.7] mm. To account for this, we ran simalma for both values at the edges of the range. An example of these channel maps can be found in Fig. A.4.

Due to the trimming that was described earlier in this section, the angular extent of the resulting synthetic channel maps is too small for CASA to run efficiently. To solve this, a border that has the exact same value as the background was added, extending the channel maps.

3.4. NaCl abundance determination

The number density of NaCl within each clump is defined relative to H₂: n(NaCl) = f_i × n(H₂), with f_i the NaCl abundance of clump i. As a first approximation, we used the NaCl abundances from the literature, which are listed in Table 4. Within a clump, the number density is assumed to be constant, and we define n(H₂) at the centre, r_centre, of each clump, using Eq. (7). For each clump, Table 5 lists the distance to IK Tau from both the centre of the clump r_centre (calculated as the median of the distances of all the points making up the clump) and the closest point to the star r_closest. At the grid points outside of the clumps, the NaCl number density has an arbitrarily small value of 10⁻⁴⁰ × n(H₂), as Magritte requires a number density greater than zero.

A grid of Magritte simulations was run for different clump abundances. For each model, the intensities of the resulting synthetic channel maps were compared to the observed ones. To make an unambiguous comparison, spectra were extracted for regions covering each of the clumps, and for a region containing all the clumps of both the observed and synthetic channel maps. These regions are plotted in Fig. 6. The following χ² statistic was then calculated to find the best model:

(8)

Fig. 6. Apertures used for extracting spectra from the observed and synthetic channel maps. The ellipses labelled A through E represent each of the clumps, and the big circle is the general emission area. Areas A, B, C, and general are also used for the rotational temperature diagram (see Appendix D). The background is the zeroth moment map of the NaCl J = 26 → 25 transition, where yellow (blue) is the highest (lowest) intensity, and the white ellipse is the beam.

Here, I is the integrated line intensity, N is the number of spectra (N = 6 for five clumps and one total region), and σ is the uncertainty on the observed spectra, which is σ = 7% (Decin et al. 2018).

We modelled transitions NaCl J = 26 → 25 and NaCl J = 27 → 26 at the same time, calculating the χ²-function across all clumps and transitions. To explore the parameter space, we first set the abundance of all clumps to f = 3.1  ×  10⁻⁷ (Velilla Prieto et al. 2017), which we find to be too large. Then, we set the abundance to f = 4  ×  10⁻⁹ (Milam et al. 2007), which we find is too small. Thereafter, to determine the best fitting abundance for each clump, we varied the abundances between these values. Figure 7 shows the spectra of the final model for both lines, and Table 6 gives the retrieved abundances and the corresponding number densities. The uncertainties in this table are for a confidence interval of 50%.

Fig. 7. Synthetic spectra of transitions NaCl J = 26 → 25 (left) and NaCl J = 27 → 26 (right), extracted from the emission region of each of the clumps, and of the general emission region. Figure 6 shows the extraction apertures used. For clump E, the only relevant peak is indicated in red as the rest of the emission is caused by the overlap with the emission regions of other clumps.

4. Discussion

4.1. De-projection

4.1.1. Choice of velocity profile

In this work we have assumed, as have previous studies (Decin et al. 2010a, 2018; Maercker et al. 2016), a beta-law expansion velocity profile (Eq. (4)) for the CSE of IK Tau. This ensures a monotonically increasing radial velocity, yielding an invertible projected velocity along the line of sight, which is necessary for our de-projection method. However, it is possible that the CSE of IK Tau has a more complex velocity field. For example, there may be shocks propagating through the CSE – especially in the inner regions close to the star where pulsations could drive shocks – or the wind may be partially shaped by a stellar or planetary companion, affecting the velocity field gravitationally. In principle, one could model these effects and as long as the resulting velocity along the line of sight is invertible, one could still employ our de-projection method. Furthermore, from the synthetic observations of the de-projected model, one could then fit the additional model parameters. However, this would greatly extend the parameter space, and it would be very difficult to impose the invertibility of the velocity along the line of sight. In future work, we will explore these possibilities with another (probabilistic) de-projection method that relies on fewer and less strict assumptions.

We tested a few different values of υ_exp in our de-projection, including υ_exp = 17.5 km s⁻¹, the value derived from single-dish observations, and values between this and our final value of υ_exp = 28.4 km s⁻¹. In general, assuming larger values of υ_exp in Eq. (5) results in a more compressed 3D distribution in the z-direction. If υ_exp is smaller than the largest absolute channel velocity considered, Eq. (5) will have an imaginary denominator and the de-projection will fail. The physical implication of this is that if the high velocity wings are created by processes (such as shocks) close to the star, which we consider to be more likely than the wings resulting from some gas being accelerated at the edge of the wind, then the emission corresponding to the high velocity wings will be placed on an outer edge in the de-projected model and will never be placed close to the star (in the z-direction). This has some consequences for the 3D structure of the NaCl distribution and, in particular, for the location of clump E, which we discuss in more detail in the following section.

Finally, we reiterate that whatever the true velocity field may be, it should be equally applicable to all the molecular emission in the CSE of IK Tau, although different molecular lines may trace different parts of the velocity field. The main driver of the CSE expansion is radiation pressure on the dust that then drags the gas with it (Höfner & Olofsson 2018). With this in mind, there is no apparent reason why NaCl should have different kinematic behaviour to other molecules. The asymmetric distribution points to an unusual formation mechanism, but the distribution of the clumps suggests that they are expanding with the rest of the CSE, especially if we consider clump C to have formed most recently, with clumps A, B, and D getting progressively older (see the more detailed discussion on clump ages below).

4.1.2. 3D structure of NaCl

In the 3D model of the NaCl distribution around IK Tau (Fig. 4), a tentative spiral-like shape can be discerned from the clumps, starting from clump C, moving outwards through clumps A and B to end in clump D. This spiral-like shape lies more or less in the xy-plane and could be caused by a companion, although none has been directly observed as of yet. The structures seen in the ALMA CO data of IK Tau are also indicative of a spiral-like structure (Decin et al. 2018). Another possible cause of the spiral-like pattern could be rotation of the star, much like the sprinkler-like scenario described by Quintana-Lacaci et al. (2023). Other possible scenarios for the formation of the observed distribution of NaCl are discussed in Sect. 4.2.3.

Although we have generally treated the clumps of NaCl as discrete in this work, their apparent separation may be an observational limitation. That is, lower abundances (corresponding to lower intensities) of NaCl may be present between the clumps we have defined here. In Fig. 4 it can already be seen that clumps C and A, and A and B are almost touching. A less abundant trail of NaCl may in fact be linking them (and also clumps B and D) while not being detected above the noise of our observations. The clumps we do detect may hence be enriched regions of NaCl along a putative spiral.

To better understand the putative spiral-like structure, we examined the maximum intensity plot of NaCl J = 26 → 25. The maximum intensity plot was made (using the CASA software, CASA Team 2022) by selecting the maximum flux for each pixel in any channel across the NaCl line (i.e., the range of channels in Fig. 2). In this way, the spatial locations of all four main clumps (A, B, C, and D) can be seen in a single plot, whereas clumps B and D are difficult to discern in an integrated intensity map such as the zeroth moment map shown in Fig. 6. In Fig. 8 we plot the radially distributed emission of the maximum intensity map as a function of angle for a circle centred on the star. We show the plot for three full revolutions and label the four visible clumps (A to D). We fit a straight line to the brightest points in each clump and then plot this fit over the maximum intensity plot in the usual RA and Dec coordinates, with the result of a spiral fit to the clumps. The fact that clump D lies close to the LSR velocity (see Fig. 2 and discussion in Sect. 4.1.2) indicates that the spiral arrangement of the clumps lies close to the plane of the sky, though the offset velocities of clumps A and C in opposite directions suggest that the spiral might be tilted along the axis between clump D and the central star, such that clump A is angled towards us and clump C is angled away from us by 31 ± 11° (with the uncertainty based on the beam size).

Fig. 8. Maximum intensity plot of NaCl J = 26 → 25 towards IK Tau. Left: radial emission distribution against the angle for the maximum intensity plot, with a full revolution repeated twice. The clumps are labelled, and an orange line is fit to the brightest point in each clump. Right: maximum intensity plot with the same line, now appearing as a spiral passing through clumps C, A, B, and D.

The relation between the size of a clump and the distance to the star can tell us something about the evolution of the clumps in time and how the NaCl abundance varies between them. The distances to each of the clumps are given in Table 5. Clump D is one of the smallest clumps and the farthest from the star, while clump C is both the largest and the closest to the star. Similarly, clump A is slightly closer to the star than clump B, and is larger as well. Given that the CSE in which these clumps lie is radially expanding, it can be assumed that the closer a clump is located to the star, the more recently it was formed, and vice versa. This implies that clump D is the oldest, and clump C is the youngest. Two possible explanations for the decreasing size with distance can be invoked: (1) the clumps expand as they travel outwards, making the outer regions of the clumps too diluted to be detected, or (2) NaCl is destroyed more readily farther from the star, either through condensation onto dust or via photodissociation. In both cases, the NaCl number density of a clump would decrease as it travels outwards. We thus expect clump D to have the lowest NaCl number density and abundance, while clump C would have the highest. We find that this holds for the NaCl number density, but not for the fractional abundance: the clumps farther from the star, such as clump D, have a higher abundance than the clumps closer to the star, such as clump C. This is discussed in more detail in Sect. 4.3.

The one outlier is clump E, which is similar in size to clump D, but which lies at a similar distance to clumps C and A. Two possible explanations for this are: (1) because the contour defining clump E is only slightly bigger than the beam size (see Sect. 2), it could possibly be misinterpreted noise; and (2) clump E could also be a part of clump C, similarly to how clumps A and B are connected. In the top-right panel of Fig. 4, it can be seen that clump C already has a similar protrusion in the negative y-direction, so clump E could be connected to clump C in the same way. To expand on this idea, it is possible that clump E and the negative-y protrusion of clump C are not located where our de-projection model (with the monotonic radial velocity profile) places them, but are instead high velocity components located close to the star. The observations presented here were taken close to the peak of IK Tau’s pulsation phase (Decin et al. 2017). It could be that these two features are part of a shock that has accelerated the gas in a relatively small region close to the star. (Such a shock could also account for the high velocity wings seen in other molecular line profiles, but a detailed analysis of other molecular emission is beyond the scope of the present work.)

4.2. Origin of the NaCl clumps

No other molecules detected towards IK Tau with ALMA show the same kind of irregular emission as NaCl. If one of the more abundant species were to display such clumpy emission, it is possible that the clumps could be obscured by the rest of the more extended and symmetric emission. However, an examination of the data, particularly of CO, SiO and HCN, does not reveal any evidence of clumps or overdensities corresponding to the NaCl clumps. Other metal-bearing molecules that have been (tentatively) detected towards IK Tau, such as TiO (Danilovich et al. 2020a), AlCl and AlOH (Decin et al. 2017), exhibit emission spatially centred on the star. This includes TiO, which has a blue-shifted spectral line profile, similar to what we see for NaCl owing to the dominance of clump A, but spatially centred emission (Danilovich et al. 2020a), unlike NaCl.

Our results can be compared to the carbon-rich AGB star CW Leo (Quintana-Lacaci et al. 2016) and the RSG VY CMa (Decin et al. 2016; Quintana-Lacaci et al. 2023), which both exhibit structures in their NaCl emission. Quintana-Lacaci et al. (2016) found that the emission of NaCl and KCl around CW Leo could be explained by either a spiral or a torus. However, this emission is predominantly centred on the star and most of the emission is part of a continuous structure in each channel. This is not the case for IK Tau and could suggest a different formation mechanism in IK Tau than for CW Leo. For VY CMa, Decin et al. (2016) found multiple clumps, whose total emission is not centred on the star, with two clumps lying at least 1″ farther from the star. They suggest that localised mass ejections that drive sputtering off dust could be the origin of these clumps. Quintana-Lacaci et al. (2023) obtained higher resolution observations of molecular emission towards VY CMa and found that the emission from NaCl and H₂S trace jets with a Hubble-like velocity field. We see no evidence of jets around IK Tau. Although there are no spatially resolved observations of H₂S towards IK Tau, a previous study based on Atacama Pathfinder Experiment (APEX) observations found much wider H₂S line profiles than those of NaCl and modelling indicated a relatively large emitting envelope centred on the star (Danilovich et al. 2017b). This is a strong indication that NaCl and H₂S do not trace the same regions around IK Tau and that one or both likely do not have the same formation mechanisms around IK Tau as around VY CMa.

For IK Tau, the irregular distribution of the NaCl clumps and their apparent spiral distribution suggests a temporally variable process is behind their formation. We considered several possible contributing scenarios for the formation of the NaCl clumps, which we outline in the following subsections. These include the impact of pulsations, correlations with the continuum and hence dust, and the chemistry of NaCl.

4.2.1. Pulsations

We tested whether the formation of the clumps could be partly driven by pulsations by estimating the time taken for each clump to reach its current location from the surface of the star. Even though the clump travel time from the stellar surface to a few stellar radii is uncertain, we are only considering the differences between clump travel times, so this region is cancelled out. We find the differences in travel times between clumps A, B and C to be close to integer multiples of the pulsation period, within the uncertainties, which could suggest they were formed during similar pulsation phases. This does not hold for clumps D and E, however. Also, the spacing is uneven, such that we do not see a clump forming from each pulsation, making this result inconclusive. If clump E and the negative-y protrusion of clump C are presently being shocked and only formed recently, that would suggest that a clump could be formed over several pulsations, since we estimate that the centre of clump C is 41 AU from the star (and the nearest edge, excluding the aforementioned protrusions, is only 3 AU from the star). Overall, we estimate that the difference in expansion time between the centres of clump C and clump D (the farthest clump) is ∼20–35 years, depending on the choice of expansion velocity. This is a remarkably short time period if we consider that the approximate amount of time it would take for the envelope to expand out to 1″ (the approximate radius of the HCN emission, Decin et al. 2018) is ∼70 years and the time taken to expand to the radius of the dust observed with the Herschel Photoconductor Array Camera and Spectrometer (PACS, 85″, Cox et al. 2012) is around 6000 years. We do not see any NaCl beyond 0.6″ or 150 AU from the star (Figs. 6 and 10), from which we can deduce that beyond this radius either the excitation conditions are not favourable for NaCl to be detected or the NaCl is removed from the gaseous phase, either because the molecule is destroyed or because it is condensed onto dust grains.

4.2.2. Continuum and dust

In Fig. 9 we plot the continuum emission of IK Tau as observed with ALMA (for details of the continuum data reduction, see Decin et al. 2018). Although the emission close to the star is almost circular (as shown by the 100σ contour), the contour at the 15σ level and the lower level contours show significant asymmetry. In the right panel of Fig. 9 we plot the same continuum contours over the maximum intensity plot of NaCl J = 26 → 25. Hence, it becomes clear that the asymmetric continuum contours partly follow the spatial locations of the NaCl clumps. The 15σ continuum contour protrudes in the direction of clump A, the 10σ contour approximately encompasses clump A and protrudes in the direction of clump B, the 5σ contour encompasses clump B and the 3σ contour protrudes past clump D. There are also two small 5σ regions aligned with the clump D protrusion but located farther from the star. This indicates that clump D is real but suggests that the excitation conditions for NaCl are such that clump D is too cool or distant from the star for NaCl J = 27 → 26 to be detected. We note also that the centre of clump D corresponds to a kinetic temperature of ∼200 K (Fig. 5), which is already slightly lower than the J = 26 level energy of 220 K. It is possible that lower-J NaCl lines might result in a clearer detection of clump D and possibly even other cooler clumps that are not visible in either of the transitions covered by the present observations.

Fig. 9. ALMA continuum observed towards IK Tau. Left: continuum flux for IK Tau as observed by ALMA plotted on a log scale and with contours drawn at levels of 3, 5, 10, 15, and 100σ. Right: maximum intensity plot of NaCl J = 26 → 25, overlaid with the same red continuum contours. The white ellipse in the bottom-right corner indicates the size of the synthetic beam for the NaCl emission. The red star indicates the location of the continuum peak, and the red ellipse in the bottom-right corner indicates the size of the synthetic beam. North is up and east is left.

4.2.3. Chemical considerations

NaCl forms through the reaction of HCl with free sodium (Na + HCl → NaCl + H), which has been measured to be fast and with a high energy barrier (∼4100 K, Husain & Marshall 1986; Cherchneff 2012), meaning that it will form more readily at high temperatures. In the absence of external forces on a spherically symmetric AGB wind, NaCl is expected to form uniformly and at relatively low abundances (≲10⁻⁸ relative to H₂) in the stellar atmosphere and inner wind at chemical equilibrium (Cherchneff 2012; Agúndez et al. 2020, and see the further discussion in Sect. 4.3), but such a uniform scenario does not explain our observations. The chemical model of Gobrecht et al. (2016), which includes shock-induced chemistry as a result of AGB pulsations, does not result in significantly higher abundances of NaCl (less than an order of magnitude increase) than the equilibrium chemical model of Agúndez et al. (2020).

Recent studies have shown that dust and gas do not form isotropically in AGB stellar atmospheres, with large convective cells likely playing a role in gas and dust anisotropies close to the star (Velilla-Prieto et al. 2023; but see also Ohnaka et al. 2016; Takigawa et al. 2017; Khouri et al. 2018). In principle, the anisotropic formation of NaCl could be explained by some dependence on convective cells; however, it is unclear why this would result in NaCl forming in only one region at a time, as we see in our data, or why there would be relatively long gaps between the enhancements in NaCl seen as clumps, as can be deduced from the distances given in Table 5 and briefly discussed in Sects. 4.1.2 and 4.2.1. Although the spiral structure discussed in Sect. 4.1.2 is only tentative, the emission is largely grouped in a plane in the CSE. Randomised ejections could form such a structure, but so could something moving in a roughly circular trajectory around IK Tau.

Decin et al. (2016) proposed that the clumpy NaCl distribution towards VY CMa could be a result of sputtering off dust, a process that could be driven by shocks. For IK Tau, the scenario of NaCl being formed on a single, rotating, location on the surface of the AGB star seems unlikely, so we suggest that the formation of NaCl could be assisted by a body orbiting the star, which would also be a source of shocked gas if it were moving faster than the local sound speed, υ_sound ≈ 3 km s⁻¹. Such a supersonic velocity is likely since the orbital velocity for a planet or brown dwarf at 50 AU from IK Tau would be 4–5 km s⁻¹ (based on the mass estimate for IK Tau of Danilovich et al. 2017a), higher than the sound speed in the innermost regions of the wind. For a circular orbit at 5 AU, this goes up to 14–16 km s⁻¹, well above the local sound speed. The shock conditions (i.e., higher temperatures and densities) generated by the orbit could then drive the NaCl formation, either through the bimolecular chemical reaction of Na and HCl or through sputtering off dust, as suggested by Decin et al. (2016), which would also explain why the NaCl clumps correlate well with the dust seen in the ALMA continuum (Fig. 9). Peaks in the AGB stellar pulsation might also contribute to heightened shock conditions and hence increased NaCl formation (see Sect. 4.2.1). The atmospheric models of Bowen (1988) predict that shocks from AGB pulsations will cause increases of a factor of a few in temperature and density in the inner wind, whereas the binary hydrodynamic models of Maes et al. (2021) and Malfait et al. (2021) predict increases up to an order of magnitude in density and temperature, depending on the model parameters. The hydrodynamic models of Aydi & Mohamed (2022), which consider both stellar pulsations and substellar companions in the inner wind, predict enhanced shock conditions when the shocked gas in the wake of a closely orbiting companion (orbital periods of the order of ∼1000 days) interacts with the shock from the pulsations of the AGB star. Such shocks also result in a shock velocity higher than the velocity of the ambient gas, which could explain the high velocities observed for clump E and the protrusion from clump C in the negative y-direction (see the further discussion in Sect. 4.1.2), whereas (portions of) clumps that formed less recently have decelerated and cooled as they moved outwards in the CSEs. The models discussed here suggest that both stellar pulsations and shocks from orbiting companions could contribute to conditions that facilitate the production of NaCl. Collisions between the AGB wind and the shock generated by the supersonic orbit could also contribute to enhanced dust formation, potentially further enhanced during pulsations, (i.e., similar to, but less intense than, the process described by Usov 1991, for Wolf-Rayet stars), which could also explain the continuum emission and why it correlates with the NaCl emission (see Fig. 9 and Sect. 4.2.2).

If the orbiting body is a star, it could also contribute to the heating of the gas through its radiation, further enabling NaCl formation by overcoming the formation energy barrier, but such a star should not be hot enough to destroy other molecules, since this is not seen in the observations. So a stellar companion to IK Tau might be, for example, a cooler main sequence star, but not a white dwarf (see a more detailed discussion of the contribution of stellar companions to the circumstellar chemistry of AGB stars in Van de Sande & Millar 2022).

A less likely scenario is a rocky planetary body in a close orbit around the AGB star that is in the process of being destroyed. However, in such a case we would expect to see water and other metal-bearing species, composed of elements common to rocky planets, in clumps similar to NaCl, which we do not (e.g., Decin et al. 2017, 2018; Danilovich et al. 2020a, covering AlOH, AlCl, AlO, H₂O, and TiO, and note that no spatially resolved observations of cool water lines presently exist, but the observed hot water lines exhibit emission centred on the star, not in clumps). So we conclude that, an orbiting body (planetary or substellar) would mainly contribute to the aforementioned shock conditions in the wind, which then drive the formation of NaCl.

4.3. NaCl abundance

The abundances and number densities we obtain for NaCl (see Table 6) depend on the pwv we used: for all clumps other than clump D, the obtained abundance is smaller at pwv = 1.7 mm than at pwv = 0.2 mm. Since a larger pwv value means that the simulated noise has a larger amplitude, the resulting lower abundance can be thought of as a correction for the higher noise. The clump abundances for the different pwv values all agree within the uncertainties, except for clump D. Figure 7 shows that the synthetic NaCl J = 26 → 25 spectra of clump D do not feature the peak that is visible in the observed spectrum, meaning that this clump is not visible in the output. This is true whether or not the CASA ALMA simulator is applied to the synthetic channel maps (see Figs. A.3 and A.4). This discrepancy can be explained by the fact that clump D was only detected in a single channel map in only one of the two observed transitions, meaning that there is only one data point in frequency or velocity space that is not noise.

In Fig. 10 we plot the fractional abundances and number densities of each clump against the distance of each clump to the star, and also the fractional abundances against the number densities. The error bars on the radial distance r indicate the nearest and farthest point of each clump to the star. This shows that while the clumps farther from the star have a higher fractional NaCl abundance (with the exception of clump D at pwv = 0.2 mm), the absolute NaCl number density does decrease with distance. This could mean that the NaCl number density decreases more slowly than the H₂ number density. It could also be an indication that NaCl formation has not run to completion in the clumps closer to the star.

Fig. 10. NaCl clump properties. The relative NaCl abundances (left) and NaCl number densities (middle) of each of the clumps are plotted as a function of the distance of said clumps to the star. The right figure shows the NaCl number densities as a function of the relative NaCl abundances. The error bars on f and n represent the 50% confidence interval, while the error bars on r show the closest and farthest points of each clump.

Comparing our abundances with the literature values in Table 4, we can see that none of the clump abundances agree with the observation-based literature values that were found by Milam et al. (2007, f ∼ 4  ×  10⁻⁹) and Velilla Prieto et al. (2017, f ∼ 3.1  ×  10⁻⁷), which are respectively smaller and larger than our values (f ∼ 9  ×  10⁻⁹ − 5  ×  10⁻⁸). Milam et al. (2007) used single-dish data from the IRAM SMT and 12 m telescopes, while Velilla Prieto et al. (2017) used single-dish data from the IRAM 30 m telescope, and they both assumed a spherical distribution of NaCl with a source size of 0.3″. In both datasets, the NaCl distribution is not spatially resolved, and the assumed source size is too small, covering only a part of clump C rather than the entire emitting region. Velilla Prieto et al. (2017) predicts that an underestimation of the emission region in their calculations would lead to an overestimation of their abundance, which we can confirm from our results.

Gobrecht et al. (2016) modelled non-equilibrium chemistry for IK Tau in shocked gas layers close to the star. They found a range of NaCl fractional abundances from f ∼ 3.7 × 10⁻¹⁰ − 2 × 10⁻⁸ in the range 1R_⋆–9R_⋆. These values only agree with our abundance for clump C, while we find the other clumps to have larger abundances. This is consistent with clump C lying partly in the region modelled by Gobrecht et al. (2016, r ∈ [1, 6]R_⋆), including the shock where part of clump C is forming, while the other clumps we observed lie farther from the star. The difference between our result and the Gobrecht et al. (2016) model can also be explained by a companion orbiting IK Tau, as suggested in Sect. 4.2.3, because the shock generated by the companion is not included in the chemical model. Agúndez et al. (2020) calculated the thermochemical equilibrium abundances of molecules in the atmospheres of M-type AGB stars, between radial distances of 1R_⋆–10R_⋆. At the outer part of their model, (r = 10R_⋆), they find an NaCl abundance of f ∼ 1  ×  10⁻⁹ relative to H₂. A part of clump C lies within 10R_⋆, but the fractional abundance we find for clump C is almost one order of magnitude larger. This is further evidence that equilibrium conditions are not responsible for the formation of NaCl and underscores the importance of shock chemistry.

5. Conclusions

In this paper we study the clumpy distribution of NaCl in the CSE of the AGB star IK Tau. First, a 3D model of the spatial distribution of NaCl around the star was obtained by de-projecting the channel maps of the NaCl J = 26 → 25 transition observed with the ALMA telescope. To do this, we derived a new value for the maximum expansion velocity, υ_exp = 28.4 ± 1.7 km s⁻¹, using the SiS J = 19 → 18, v = 0 emission line in the spectrum of IK Tau. The 3D model indicates that each of the clumps has a different shape and size and is located at a different distance from the star. The overall distribution of these clumps suggests a tentative spiral shape, and the relative size of these clumps and their distances from the central star suggest that the farthest clumps were formed earliest and are now smaller either because of dust condensation or because they expand while moving outwards, making their outer regions too diluted to be detected. We speculate that the formation of NaCl is driven by shocks originating from stellar pulsations and/or an orbiting body, explaining the clumpy spiral distribution.

We derived the NaCl abundances of each of the clumps using the 3D radiative transfer modelling code Magritte, with our de-projected model as input. This modelling was performed under the assumption of LTE to keep the computational times reasonable. The clump abundances lie in the range 9  ×  10⁻⁹–5  ×  10⁻⁸ relative to H₂. The abundance of the clump closest to the star falls within the range predicted by shock-induced chemistry models. The higher relative abundances of the clumps farther from the star indicate that NaCl formation may continue at larger distances.

Acknowledgments

The authors grateful to the anonymous referee for providing many constructive comments. We also thank Daniel Price for fruitful discussions regarding shocks. This paper makes use of the following ALMA data: ADS/JAO.ALMA2013.1.00166.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Credit CASA: International consortium of scientists based at the National Radio Astronomical Observatory (NRAO), the European Southern Observatory (ESO), the National Astronomical Observatory of Japan (NAOJ), the CSIRO Australia Telescope National Facility (CSIRO/ATNF), and the Netherlands Institute for Radio Astronomy (ASTRON) under the guidance of NRAO. TD and FDC acknowledge support from the Research Foundation Flanders (FWO) through grants 12N9920N, and 1253223N, respectively. TD is supported in part by the Australian Research Council through a Discovery Early Career Researcher Award (DE230100183). LD thanks the support of the Fund of Scientific Research via the Senior Research Project G099720N.

References

Appendix A: Additional plots of NaCl

A.1. Spectral lines

In Fig. A.1 we plot all four lines of NaCl covered by the ALMA spectral scan of IK Tau Decin et al. (2018) using four different extraction apertures: circular apertures with radii of 800 mas and 320 mas centred on the star, and two irregular elliptical apertures, chosen in CASA to capture the majority of the emission from clumps A and C, the two brightest clumps. Further discussion can be found in Sect. 2.

Fig. A.1. Four spectral lines of NaCl covered by ALMA. The transitions are given at the top of each column and the extraction apertures to the right of each row. The spectra in the top two (black) rows were extracted for circular apertures (with radii of 800 and 320 mas) centred on the star, while the bottom two rows show spectra extracted for irregular apertures encompassing clumps A and C, respectively. The vertical dashed grey lines indicate υLSR = 34 km s−1.

A.2. Observed channel maps

Figure A.2 shows the observed channel maps of the NaCl J = 27 → 26 transition. This figure is discussed in Sect. 2.

Fig. A.2. Channel maps of the NaCl J = 27 → 26 line. The red star indicates IK Tau’s position, and the white ellipse shows the beam size. The white contours are at 3σrms, 5σrms, and 10σrms (σrms = 2.6 mJy). The υLSR = 34 km s−1 has been subtracted from the velocities.

A.3. Synthetic channel maps

Figures A.3 and A.4 both show the synthetic channel maps of the same model, before and after the CASA simulator has been applied to them, respectively. These figures are mentioned in Sects. 3.3 and 4.3.

Fig. A.3. Synthetic Magritte channel maps of the NaCl J = 26 → 25 line, after the observed velocity bins are recreated by combining multiple maps. The red star indicates IK Tau’s position.

Fig. A.4. Synthetic Magritte channel maps of the NaCl J = 26 → 25 line from Fig. A.3, after the CASA simulator has been applied to them, using pwv = 0.2 mm. The red star indicates IK Tau’s position, and the white ellipse shows the (simulated) beam size.

Appendix B: Spectral lines of CO and SiS

Figure B.1 shows the observed spectral lines of the CO 3 → 2 and SiS 19 → 18 (v = 0) transitions. They are used in Sect. 3.1 to determine the maximum expansion velocity of IK Tau.

Fig. B.1. Spectral lines corresponding to the CO 3 → 2 (top left), SiS 19 → 18 (v = 0) (top right), HCN 4 → 3 (bottom left), and CS 7 → 6 (v = 0) (bottom right) transitions towards IK Tau. The dashed red lines indicate the υexp = ±17.5 km s−1. The used aperture has a radius of 320 mas.

Appendix C: De-projection

C.1. Derivation of the de-projection equation

An equation for the z coordinate can be derived from the fact that the observed channel map velocity is the component of the expansion velocity along our line of sight. If we assume a spherically symmetric stellar wind with a constant expansion velocity υ_exp, we can then define θ as the angle between the channel map velocity υ_channel, from which the υ_LSR has been subtracted, and the expansion velocity, υ_exp, such that υ_channel = υ_expcos(θ). This is illustrated in Fig. C.1. The expression can then be rewritten as

(C.1)

Fig. C.1. Principles of de-projection from radial velocity to spacial coordinates, using a constant radial expansion velocity, υexp. The observer is located on the right. The impact parameter is .

For any point on a particular channel map, this angle θ also connects the coordinate along our line of sight z with the radial distance r to said point, such that z = rcos(θ), as shown in Fig. C.1. Using the impact parameter , we can write the following expression:

(C.2)

Combining Eqs. (C.1) and (C.2) yields the expression for the coordinate along our line of sight:

(C.3)

We also note that, from the definitions of υ_channel and z, that this equation can also be written as

(C.4)

C.2. De-projection of planes of constant velocity along the line of sight

Figure C.2 shows the results of a general de-projection. In this figure, five planes of constant velocity, evenly spaced along the line of sight, are de-projected. Only the positive x side of each plane is shown. On each plane, the points with a small impact parameter are mapped to smaller z values, and thus closer to the star, while points with a larger impact parameter p are mapped to much greater z values. The larger the velocity associated with the plane, the larger the spread in z values will be. Since every channel map has a constant velocity everywhere, the de-projections of these planes are a good representation of how the contents of the channel maps are de-projected: on a single channel map, the points that lie close to the star will be de-projected to smaller z values than the points that lie farther from the map. This effect is enhanced for channel maps with velocities that are much greater than zero.

Fig. C.2. Two views on 3D images of de-projected planes. We show five de-projected planes where every point within a plane has the same velocity. Planes with the same colour have the same velocity with opposite signs. For clarity, the planes are only plotted for positive x values.

Appendix D: Excitation analysis

D.1. Calculation

Rotational temperature diagrams or population diagrams can be used to deduce rotational temperatures and column densities (Goldsmith & Langer 1999). This method assumes that the emitting region is in LTE and is optically thin. We expect NaCl to be optically thin, because we expect relatively low abundances (≲10⁻⁷). According to the Magritte modelling, the maximum optical depth is τ = 6  ×  10⁻⁵. The assumption of LTE might not hold for the clumps, due to the low density of the CSE (∼10¹¹ − 10¹² cm⁻³) and the high Einstein A-values. The rotational temperature T_rot is the excitation temperature of the purely rotational transitions.

In LTE, the column density of each level is proportional to its population, and can be expressed as

(D.1)

Here, N_u is the column density of the upper level u in cm^-2, N is the total column density of this particular species in cm^-2, Z = ∑_iN_i is the partition function, g_u is the statistical weight of the level u, E_u is the energy of the level u in eV, k is the Boltzmann constant in eV K⁻¹ and T_rot is the excitation temperature in K.

This equation can then be rewritten as

(D.2)

The rotational temperature diagram or population diagram is then a scatter plot with (E_u/k) on the x-axis and ln(N_u/g_u) on the y-axis. Under the assumptions of LTE and optically thin lines, the data points should form more or less a straight line. The rotational temperature T_rot can then be found from the slope, and the column density N can be found from the intercept (as long as the partition function Z is known). We used this method to obtain the rotational temperatures of the NaCl clumps.

To calculate N_u for each transition, the method of Saberi et al. (2022) was used. Here, N_u is defined as

(D.3)

Here, D is the distance to the star, r_e is the radius of the emitting region, W is the line flux in units (W/m²), A_ul is the Einstein coefficient of spontaneous emission for a transition linking levels u and l, h is Planck’s constant and ν is the frequency of the transition. Figure 6 shows the locations of the emission regions for clumps A, B and C, and for a larger region encompassing all the emission. We exclude clumps D and E from the rotational diagram analysis because they are not seen above the noise in the NaCl J = 27 → 26 channel maps. The emission regions are ellipses with major and minor axes given in Table D.1. W is found by integrating over the spectral line extracted for each emission region. An error of 7% on the flux density (Decin et al. 2018) and an error of 0.02″ on the radii of the emission regions are taken into account to find the uncertainties on the rotational temperature. Figure D.1 shows the rotational temperature diagram calculated using the regions shown in Fig. 6. We find rotational temperatures ranging from 14 – 49 K, as listed in Table D.1. The uncertainties given in this table are underestimates, however, because there are only two data points, and thus the slopes were calculated directly rather than fitted. If more NaCl lines were observed, more accurate rotational temperatures could be derived.

Fig. D.1. Rotational temperature diagram for four different emitting regions of NaCl around IK Tau.

D.2. Discussion

Comparing the T_rot of the three different clumps, it becomes clear that the T_rot is larger for clumps nearer to the star than for those farther from the star. This makes sense, since the physical temperature is highest near the star. That the temperature of the entire emission region lies near the average of the three clumps makes sense as well, since it takes the entire region into account. Our results are very sensitive to noise, because they were derived from only two data points, because there were only two NaCl transitions in the observed range. It is thus necessary to observe more transitions before conclusive results can be obtained.

Velilla Prieto et al. (2017) determined a rotational temperature of T_rot = 67 ± 7 K for NaCl, using single-dish data from the IRAM 30m telescope, and assumed a circular NaCl emitting region with a diameter of 0.3″. The disagreement between their and our results can be explained by two factors, the first being that, as mentioned before, in their observations, the NaCl distribution is not spatially resolved and their assumed emission region only covers a part of clump C. A second factor is that we only had two data points available, while they used 13, but the first factor overshadows this benefit. The diameter of their assumed emission region is four times smaller than the diameter of our general emission region, making its area too small by a factor of 16. They were thus unable to take most of the emission into account for their calculations. Interestingly, the T_rot we derived for the separate clumps lie in a similar range as the T_rot found by Velilla Prieto et al. (2017) for most other molecules in the CSE of IK Tau. They found that most molecules, other than NaCl and SiS, display a rotational temperature between 15 K and 40 K in the CSE of IK Tau. Whether this means that the T_rot of NaCl should actually be closer to this range can only be determined with more spatially resolved observations of more NaCl lines.

Finally, the derived rotational temperatures are much lower than those of the kinetic temperature profile of IK Tau derived by Maercker et al. (2016), at similar radial distances (this temperature profile is shown in Fig. 5). This is especially so for clump C, which lies very close to the star, yet has a rather cool rotational temperature of T_rot = 48.5 ± 0.4 K. This large discrepancy indicates that the assumption of LTE does not hold, as in LTE, the rotational temperature is equal to the gas kinetic temperature (Saberi et al. 2022). These results should be further verified by observing more spatially resolved transitions of NaCl.

All Tables

All Figures

	Fig. 1. Four NaCl lines covered in the spectral scan of IK Tau by Decin et al. (2018). The lines were extracted using a circular aperture with a radius of 320 mas, centred on the star. Spectra for additional extraction apertures are shown in Fig. A.1.
In the text

	Fig. 2. Channel maps of the NaCl J = 26 → 25 line. The red star indicates IK Tau’s position, and the white ellipse shows the beam size. The white contours are at 3σrms, 5σrms, and 10σrms (σrms = 2.4 mJy). The υLSR = 34 km s−1 has been subtracted from the velocities.
In the text

	Fig. 3. Line profiles of SiS. Left: maximum expansion velocity determination using the SiS 19 → 18 (v = 0). The horizontal dotted orange line is the used criterion, 3σrms, and the vertical dashed lines indicate the red and blue wing velocities. Right: same as on the left, but zoomed in to better discern the red and blue wings. The method used is from Gottlieb et al. (2022).
In the text

	Fig. 4. De-projection of the NaCl J = 26 → 25 channel maps shown from different angles. Only points detected at ≥3σ are shown, and the colours indicate intensity, using the same colour bar as for the original channel maps (Fig. 2). The yellow arrow points in the direction of the line of sight.
In the text

	Fig. 5. Physical conditions in the CSE of IK Tau. Top: temperature profile derived by Maercker et al. (2016), extrapolated to the stellar radius and effective temperature of IK Tau. Bottom: number density profile of H2 calculated from the mass-loss rate using Eq. (7).
In the text

Fig. 6. Apertures used for extracting spectra from the observed and synthetic channel maps. The ellipses labelled A through E represent each of the clumps, and the big circle is the general emission area. Areas A, B, C, and general are also used for the rotational temperature diagram (see Appendix D). The background is the zeroth moment map of the NaCl J = 26 → 25 transition, where yellow (blue) is the highest (lowest) intensity, and the white ellipse is the beam.

In the text

	Fig. 7. Synthetic spectra of transitions NaCl J = 26 → 25 (left) and NaCl J = 27 → 26 (right), extracted from the emission region of each of the clumps, and of the general emission region. Figure 6 shows the extraction apertures used. For clump E, the only relevant peak is indicated in red as the rest of the emission is caused by the overlap with the emission regions of other clumps.
In the text

	Fig. 8. Maximum intensity plot of NaCl J = 26 → 25 towards IK Tau. Left: radial emission distribution against the angle for the maximum intensity plot, with a full revolution repeated twice. The clumps are labelled, and an orange line is fit to the brightest point in each clump. Right: maximum intensity plot with the same line, now appearing as a spiral passing through clumps C, A, B, and D.
In the text

Fig. 9. ALMA continuum observed towards IK Tau. Left: continuum flux for IK Tau as observed by ALMA plotted on a log scale and with contours drawn at levels of 3, 5, 10, 15, and 100σ. Right: maximum intensity plot of NaCl J = 26 → 25, overlaid with the same red continuum contours. The white ellipse in the bottom-right corner indicates the size of the synthetic beam for the NaCl emission. The red star indicates the location of the continuum peak, and the red ellipse in the bottom-right corner indicates the size of the synthetic beam. North is up and east is left.

In the text

Fig. 10. NaCl clump properties. The relative NaCl abundances (left) and NaCl number densities (middle) of each of the clumps are plotted as a function of the distance of said clumps to the star. The right figure shows the NaCl number densities as a function of the relative NaCl abundances. The error bars on f and n represent the 50% confidence interval, while the error bars on r show the closest and farthest points of each clump.

In the text

Fig. A.1. Four spectral lines of NaCl covered by ALMA. The transitions are given at the top of each column and the extraction apertures to the right of each row. The spectra in the top two (black) rows were extracted for circular apertures (with radii of 800 and 320 mas) centred on the star, while the bottom two rows show spectra extracted for irregular apertures encompassing clumps A and C, respectively. The vertical dashed grey lines indicate υLSR = 34 km s−1.

In the text

	Fig. A.2. Channel maps of the NaCl J = 27 → 26 line. The red star indicates IK Tau’s position, and the white ellipse shows the beam size. The white contours are at 3σrms, 5σrms, and 10σrms (σrms = 2.6 mJy). The υLSR = 34 km s−1 has been subtracted from the velocities.
In the text

	Fig. A.3. Synthetic Magritte channel maps of the NaCl J = 26 → 25 line, after the observed velocity bins are recreated by combining multiple maps. The red star indicates IK Tau’s position.
In the text

	Fig. A.4. Synthetic Magritte channel maps of the NaCl J = 26 → 25 line from Fig. A.3, after the CASA simulator has been applied to them, using pwv = 0.2 mm. The red star indicates IK Tau’s position, and the white ellipse shows the (simulated) beam size.
In the text

	Fig. B.1. Spectral lines corresponding to the CO 3 → 2 (top left), SiS 19 → 18 (v = 0) (top right), HCN 4 → 3 (bottom left), and CS 7 → 6 (v = 0) (bottom right) transitions towards IK Tau. The dashed red lines indicate the υexp = ±17.5 km s−1. The used aperture has a radius of 320 mas.
In the text

	Fig. C.1. Principles of de-projection from radial velocity to spacial coordinates, using a constant radial expansion velocity, υexp. The observer is located on the right. The impact parameter is .
In the text

	Fig. C.2. Two views on 3D images of de-projected planes. We show five de-projected planes where every point within a plane has the same velocity. Planes with the same colour have the same velocity with opposite signs. For clarity, the planes are only plotted for positive x values.
In the text

	Fig. D.1. Rotational temperature diagram for four different emitting regions of NaCl around IK Tau.
In the text