The morpho-kinematics of the circumstellar envelope around the AGB star EP Aqr

TL;DR

This study uses ALMA observations to analyze the complex morpho-kinematics of the circumstellar envelope around the AGB star EP Aqr, revealing detailed structures and dynamics that inform models of stellar mass loss.

Contribution

It introduces an original method for separating data-cube components, providing new insights into the envelope's morphology, kinematics, and inhomogeneities around EP Aqr.

Findings

Mass loss rate of (1.6±0.4)×10⁻⁷ solar masses per year

Detection of a spiral enhancement at low velocities

Identification of a circular enhancement in expansion velocity

Abstract

ALMA observations of CO(1-0) and CO(2-1) emissions of the circumstellar envelope of EP Aqr, an oxygen-rich AGB star, are reported. A thorough analysis of their properties is presented using an original method based on the separation of the data-cube into a low velocity component associated with an equatorial outflow and a faster component associated with a bipolar outflow. A number of important and new results are obtained concerning the distribution in space of the effective emissivity, the temperature, the density and the flux of matter. A mass loss rate of (1.6$\pm$0.4)10$^{-7}$ solar masses per year is measured. The main parameters defining the morphology and kinematics of the envelope are evaluated and uncertainties inherent to de-projection are critically discussed. Detailed properties of the equatorial region of the envelope are presented including a measurement of the line width and a precise description of the observed inhomogeneity of both morphology and kinematics. In particular, in addition to the presence of a previously observed spiral enhancement of the morphology at very small Doppler velocities, a similarly significant but uncorrelated circular enhancement of the expansion velocity is revealed, both close to the limit of sensitivity. The results of the analysis place significant constraints on the parameters of models proposing descriptions of the mass loss mechanism, but cannot choose among them with confidence.

Tables

Table 1. Table 1: Main observation parameters. R 𝑅 R is the distance from the star projected on the sky plane.

	CO(1-0)	CO(2-1)
Beam FWHM (arcsec2)	0.78$\times$0.70	0.33$\times$0.30
Beam PA (∘)	$-$56	$-$80
Maximum baseline (m)	1124	1400
Minimum baseline (m)	19	15
Time on source (min)	98.3	52.4
Zero spacing	Yes	No
Area covered by mosaic (arcsec2)	$\sim 60\times 45$	$\sim 60\times 45$
Noise (mJy beam-1/(0.2 km s-1))	7	6
Flux, 13 arcsec (Jy km s-1)	110	549

Table 2. Table 2: Result of the fits of the mean Doppler velocity of the narrow, broad, blue-shifted broad and red-shifted broad components of the form < V z > = a + b ​ cos ⁡ ( ψ − ψ 0 ) expectation subscript 𝑉 𝑧 𝑎 𝑏 𝜓 subscript 𝜓 0 <\!\!V_{z}\!\!>=a+b\cos(\psi-\psi_{0}) . For the blue-shifted broad component we use < − V z > expectation subscript 𝑉 𝑧 <\!\!-V_{z}\!\!> rather than < V z > expectation subscript 𝑉 𝑧 <\!\!V_{z}\!\!> .

		$a$ (km s-1)	$b$ (km s-1)	${\psi }_0$ (degree)
CO(1-0)	Narrow component	0	0.30	$-$4.1
	Broad component	0.44	$-$0.61	3.8
	Blue-shifted broad component	5.2	0.32	16
	Red-shifted broad component	5.4	$-$0.36	9
CO(2-1)	Narrow component	0.02	0.36	$-$5.0
	Broad component	0.50	$-$0.68	1.8
	Blue-shifted broad component	4.5	0.39	16
	Red-shifted broad component	4.7	$-$0.34	$-$4

Table 3. Table 3: Parameters of the fits I 0 × 10 − R / R 0 subscript 𝐼 0 superscript 10 𝑅 subscript 𝑅 0 I_{0}\times 10^{-R/R_{0}} to the R 𝑅 R -dependence of the intensity ( R > 2.5 𝑅 2.5 R>2.5 arcsec).

	Line	$I_0$ (Jy beam-1 km s-1)	$R_0$ (arcsec)
Narrow component	CO(1-0)	0.076	16.2
	CO(2-1)	0.177	9.0
Broad component	CO(1-0)	0.567	11.7
	CO(2-1)	1.18	6.1

Table 4. Table 4: Results of Gaussian fits to the Φ ± subscript Φ plus-or-minus \mit{\Phi}_{\pm} and δ ​ V ± 𝛿 subscript 𝑉 plus-or-minus \delta{V_{\pm}} distributions.

${\Phi }_{+}$		${\Phi }_{-}$		$\delta V_{+}$ (km s-1)		$\delta V_{-}$ (km s-1)
Mean	$\sigma$	Mean	$\sigma$	Mean	$\sigma$	Mean	$\sigma$
1.01	0.32	0.01	0.15	0.03	0.26	$-$0.02	0.14

Table 5. Table 5: Parameters of the fits of the form I 1 ​ exp ⁡ ( − R / R 1 − k 1 ​ R 2 ) + I 2 ​ exp ⁡ ( − R / R 2 − k 2 ​ R 2 ) subscript 𝐼 1 𝑅 subscript 𝑅 1 subscript 𝑘 1 superscript 𝑅 2 subscript 𝐼 2 𝑅 subscript 𝑅 2 subscript 𝑘 2 superscript 𝑅 2 I_{1}\exp(-R/R_{1}-k_{1}R^{2})+I_{2}\exp(-R/R_{2}-k_{2}R^{2}) to the R 𝑅 R -dependence of the integrated intensity for R < 8 𝑅 8 R<8 arcsec.

	Line	$I_1$ (Jy beam-1 km s-1)	$I_2$ (Jy beam-1 km s-1)	$R_1$ (arcsec)	$R_2$ (arcsec)	$k_1$ (arcsec-2)	$k_2$ (arcsec-2)
Narrow component	CO(1-0)	0.170	0.0593	2.69	24.4	0.458	0.00945
	CO(2-1)	0.334	0.141	2.38	6.20	0.687	0.00895
Broad component	CO(1-0)	0.781	0.506	2.93	6.58	1.65	0.00400
	CO(2-1)	4.00	0.875	0.347	3.99	0.138	0.0122

Table 6. Table 6: Narrow component line widths (Gaussian σ 𝜎 \sigma ’s, km s -1 ) measured after correction. Distances on the sky plane are measured in arcsec.

Selection		CO(1-0)	CO(2-1)
All ($R$ & $\psi$ corrected)		0.51	0.55
All (each px recentred)		0.44	0.47
$\psi \pm {22.5}^{\circ }$	0∘	0.42	0.45
	45∘	0.46	0.48
	90∘	0.43	0.49
	135∘	0.41	0.44
	180∘	0.39	0.42
	225∘	0.45	0.45
	270∘	0.43	0.46
	315∘	0.42	0.44
, $R>1.5$		0.51	0.54
$\left|R-R_{sp}\right|>h/4$, $R>1.5$		0.56	0.61
, $R>1.5$		0.56	0.61
$\left|R-R_{circle}\right|>0.95$, $R>1.5$		0.50	0.54
		0.40	0.42

Table 7. Table 7: Parameters of the reference distributions used for the broad component. They are defined as F 1 ​ ( r ) subscript 𝐹 1 𝑟 F_{1}(r) being a fit to ρ ​ r 2 𝜌 superscript 𝑟 2 \rho r^{2} of the form exp ⁡ ( a + b ​ r + c ​ r 2 + d ​ r 3 + e ​ r 4 ) 𝑎 𝑏 𝑟 𝑐 superscript 𝑟 2 𝑑 superscript 𝑟 3 𝑒 superscript 𝑟 4 \exp(a+b\,r+c\,r^{2}+d\,r^{3}+e\,r^{4}) and F 2 ​ ( sin ⁡ α ) subscript 𝐹 2 𝛼 F_{2}(\sin\!\alpha) being a fit to ρ ​ r 2 / F 1 ​ ( r ) 𝜌 superscript 𝑟 2 subscript 𝐹 1 𝑟 \rho r^{2}/F_{1}(r) of the form a + b ​ sin ⁡ α + c ​ sin 2 ⁡ α + d ​ sin 4 ⁡ α 𝑎 𝑏 𝛼 𝑐 superscript 2 𝛼 𝑑 superscript 4 𝛼 a+b\sin\!\alpha+c\sin^{2}\!\alpha+d\sin^{4}\!\alpha .

	Parameter	CO(1-0)		CO(2-1)
		arcsec	$r>4$ arcsec
$F_1\left(r\right)$	$a$	$-$1.918	$-$1.632	0.514
	$b$	$-$0.0405	0.206	0.00398
	$c$	0.149	$-$0.0118	0.105
	$d$	$-$0.0206	0	$-$0.0229
	$e$	0	0	0.00125
$F_2\left(\sin \alpha \right)$	a	0.370		0.304
	b	0.164		0.278
	c	4.755		5.117
	d	$-$3.882		$-$3.919

Equations

$$\begin{array}{rl}\delta V_{+}& =\frac{1}{2}[\delta V(\text{mbox}CO(1-0))+\delta V(\text{mbox}CO(2-1))],\\ \delta V_{-}& =\frac{1}{2}[\delta V(\text{mbox}CO(1-0))-\delta V(\text{mbox}CO(2-1))],\\ \text{mit}{\Phi }_{+}& =\frac{1}{2}[\text{mit}\Phi (\text{mbox}CO(1-0))+\text{mit}\Phi (\text{mbox}CO(2-1))]\text{mbox}and\\ \text{mit}{\Phi }_{-}& =\frac{1}{2}[\text{mit}\Phi (\text{mbox}CO(1-0))-\text{mit}\Phi (\text{mbox}CO(2-1))].\end{array}$$

Full text

††pubyear: 2018††pagerange: The morpho-kinematics of the circumstellar envelope around the AGB star EP Aqr–7

1 Introduction

EP Aqr is an oxygen-rich M type AGB star at a distance of only 119$\pm$6 pc from the solar system (van Leeuwen, 2007; Gaia Collaboration, 2018) with a luminosity of 3450 solar luminosities. There is a silicate feature at 10 $\mu$m (Speck et al., 2000), which shows that the star is presently undergoing mass loss. Several observations suggest a mass loss episode at the scale of 104 to 105 years: the Herschel 70 $\mu$m observation (Cox et al., 2012) of a large trailing circumstellar shell and the HI observation (Le Bertre & Gérard, 2004) of its interaction with the interstellar medium. In CO rotational lines, EP Aqr is characterized by a composite profile with narrow ($\sim 2$ km s*-1*) and broad ($\sim 10$ km s*-1*) components (Knapp et al., 1998; Winters et al., 2003). Few oxygen-rich AGB stars have been observed with sufficient sensitivity and spectral resolution to reveal such composite profiles in CO (or SiO) rotational lines, making it difficult to reliably evaluate which fraction of the total population they represent and to state whether it corresponds to an intrinsic property of some stars, or to a common stage in the evolution of oxygen-rich AGB stars. Detailed studies of the kinematics and 3D-structure of the inner circumstellar shells ($<$ 100 AU) are needed. We are making efforts to obtain and model high spectral-resolution interferometric data on such objects Le Bertre et al. (2016). EP Aqr is particularly interesting because, from the absence of technetium in the spectrum (Lebzelter & Hron, 1999), and from its low 12C/13C ratio (Cami et al., 2000), it appears to be in an early stage of evolution on the AGB.

Spatially resolved observations obtained at IRAM (Winters et al., 2003; Nhung et al., 2015a) on EP Aqr, and on other sources such as RS Cnc (Libert et al., 2010; Hoai et al., 2014; Nhung et al., 2015b) and X Her (Kahane & Jura, 1996) are generally interpreted as the superposition of a slowly expanding wind with a bipolar outflow. The presence of irregularities in the morpho-kinematics has been interpreted as suggesting an episode of increased mass loss rate (Knapp et al., 1998; Nakashima, 2006; Winters et al., 2007). De-projection (Diep et al., 2016; Nhung et al., 2015a) has shown that morphology and kinematics of EP Aqr display approximate axi-symmetry about an axis making a small angle with the line of sight (typically smaller than $\sim 20$∘) and projecting on the sky plane some 20∘ to 30∘ west of north. This was corroborated using recent ALMA observations having four times better spatial resolution than the preceding IRAM observations by Nhung et al. (2018a, b) and Homan et al. (2018). The latter authors have also observed SO2 and SiO line emissions exploring shorter distances to the star.

The present work presents a detailed study of the morphology and kinematics of the circumstellar envelope of EP Aqr, taking due account of what has been learned from previous studies and using observations of CO(1-0) and CO(2-1) emissions obtained at ALMA. The two components of the wind, bipolar and equatorial, are studied separately, an approach that is used for the first time and that helps greatly with identifying and describing properties of the morpho-kinematics more thoroughly than could be done previously. The method is described in Section 3; Section 4 studies the equatorial component and Section 5 studies the bipolar component. The results are then used in Section 6 to present a joint description independent of the separation into two components and accounting for absorption and temperature effects. The under-determination inherent to the interpretation of radio observations and their de-projection in space plays a particularly important role in the case of EP Aqr as has been amply demonstrated by earlier studies. The emphasis of the present work is accordingly to attempt a critical and thorough study of the morpho-kinematics with the aim to identify those features that can be reliably ascertained. The results provide important new information on the properties of the circumstellar envelope at distances from the star exceeding $\sim 250$ au, where CO molecules are particularly efficient tracers. Understanding the morpho-kinematics at shorter distances requires in addition different tracers, such as SiO and SO2, a study that is postponed to a forthcoming publication (Tuan-Anh et al. , in preparation). The results are summarized and discussed in Section 7.

2 Observations and data reduction

The observations used in the present article were made in cycle 4 of ALMA operation (2016.1.00026.S) between October 30th 2016 and April 5th 2017 and were presented in detail in Nhung et al. (2018a). CO(1-0) emission was observed in 3 execution blocks in mosaic mode (3 pointings) with the number of antennas varying between 41 and 45. CO(2-1) emission was observed in 2 execution blocks in mosaic mode (10 pointings) with the number of antennas varying between 38 and 40.

Both lines were observed in two different configurations, C40-2 and C40-5 and the data were then merged in the uv plane. The main parameters attached to these observations are listed in Table 1. The CO(1-0) and CO(2-1) data have been calibrated using CASA111http//casa.nrao.edu and the mapping was done with GILDAS222https://www.iram.fr/IRAMFR/GILDAS.

We also included IRAM 30-m observations of CO(1-0) emission that have been reported in Winters et al. (2007). We used GILDAS to re-centre the map and re-sample velocities in the same way as for ALMA CO(1-0) data. We then produced a list of pseudo visibilities containing short-spacing information and merged these with the 12-m ALMA data.

The origin of coordinates at RA=21h46m31.848s and DEC=$-$02∘12’45.93” corresponds to year 2000. Between 2000 and the time of observation, the source has moved by 0.40 arcsec east and 0.31 arcsec north (proper motion of (24.98, 19.54) mas/year (van Leeuwen, 2007)); the data have been corrected accordingly. The spectral resolution (channel) has been smoothed to 0.2 km s*-1* and the Doppler velocity covers between $-20$ and 20 km s*-1*. We use as origin the Doppler velocity of $-33.6$ km s*-1* (LSR) about which the profile is well symmetric.

In addition to line emission, we have also observed the continuum emission in the same frequency bands (2.6 mm and 1.3 mm). As illustrated in the left panels of Figure 1, the source is unresolved, with no hint at the presence of a possible companion. It is centred at the same location as CO line emission to within 30 mas. The noise is at the level of 23 $\mu$Jy for band 3 and 75 $\mu$Jy for band 6. The continuum fluxes are measured as 4.9$\pm$0.1 mJy at 2.6 mm and 17.8$\pm$0.2 mJy at 1.3 mm, consistent with the result of Winters et al. (2007) and with expectation from the black-body emission of the central star. We note that Homan et al. (2018) erroneously quote a 75 times larger flux, implying surprisingly slow winds and unusual amounts of central dust.

3 Separation of the data cubes into broad and narrow components

3.1 General properties

In order to ease the analysis and interpretation of the observations, we use a system of coordinates rotated clockwise by 20∘ about the line of sight, meaning that the $y$ axis points 20∘ west of north and the $x$ axis 20∘ north of east, the $z$ axis pointing away from us along the line of sight. Indeed, earlier analyses (see in particular Nhung et al., 2018b; Homan et al., 2018) have given evidence for an approximate symmetry of the data cubes with respect to the plane containing the line of sight and projecting on the sky plane some 20∘ to 30∘ west of north. A summary of the geometry and of the variables used in the study is presented in Appendix A. Under the hypothesis that the morphology and kinematics of the star obey axi-symmetry, the star axis would project on the sky plane along the $y$ axis. The arguments that are developed below confirm this result. Note that strictly speaking, what we call “star axis” for brevity should be called “axis of the circumstellar envelope” as no evidence is presented in the present work for these to be identical. Moreover, we apply some very small shifts to the Doppler velocity scales, by adding 0.12 km s*-1* to the CO(1-0) velocities and subtracting 0.06 km s*-1* from the CO(2-1) velocities; these are smaller than the spectral resolution (the velocity is given in bins of 0.2 km s*-1* for both CO(1-0) and CO(2-1) emissions) but fine tune the mean position of the narrow central peak of the Doppler velocity spectrum to be at the origin ($V_{z}=0$). Finally, in order to ease the comparison between the CO(1-0) and CO(2-1) observations, we normally use a common pixel size, $0.3\times 0.3$ arcsec2 in spite of the difference between the associated beam widths. The grouping is done in the image plane after cleaning, the latter operation being made as usual with optimal pixel sizes, at the level of $\sim 20$ to 25% of the beam size. The present analysis requires a good spectral resolution but is less demanding on spatial resolution, most resolved features being significantly larger than the beam. The new data cubes are obtained from the original data cubes by sharing the content of an original pixel between the new pixels that overlap it, in proportion to the overlap fraction of the original pixel area.

Figure 1 (right panels) displays the distributions of the brightness measured in each pixel over the circle of projected distances to the star not exceeding 12 arcsec. It gives evidence for an effective noise level ($1\sigma$) of 6.7 mJy beam*-1* for CO(1-0) and 5.8 mJy beam*-1* for CO(2-1), consistent with thermal noise levels between 1 and 1.5 mJy per 0.2 km s*-1* bin. As expected from the pixel/beam ratios, the grouping of pixels has only little effect on the effective noise. Before grouping, these values were respectively 7.2 and 6.7 mJy beam*-1*, namely respectively $\sim 7$% and $\sim 16$% larger than after grouping.

Figure 2 displays the distributions of Doppler velocity, $V_{z}$. In this figure, like in the remainder of the present work, we limit our analysis to pixels distant by less than 8 arcsec from the central star. The rightmost panel of Figure 2 compares the CO(1-0) and CO(2-1) spectra normalised to a same area, giving evidence for a significant difference of shape near the end-points. This feature has been discussed in Nhung et al. (2018a) and associated with a depression of the effective emissivity near the poles of the star. We discuss it in some detail in Section 5.2.

3.2 Splitting the data cubes

Earlier analyses have associated, at least qualitatively, the narrow central peak of the Doppler velocity spectra with an equatorial outflow being seen nearly face-on, and the wings of the spectra with a bipolar outflow being seen nearly pole-on. We exploit this remark by separating the data cubes in two parts, referred to as narrow and broad components and associated with equatorial and bipolar outflows respectively. To do so, we use in each pixel two control-regions, each 2 km s*-1* wide, located symmetrically with respect to the position of the central peak in the pixel and covering between 2 km s*-1* and 4 km s*-1* away from it on either side. We then interpolate linearly between the two control-regions what we attribute to the broad component, the narrow component being obtained by subtraction. We have checked thoroughly that the results are very robust and that the algorithm being used does not produce any significant bias. Indeed most of these results would not be strongly affected if the narrow component were simply defined as obeying the inequality $|V_{z}|<2$ km s*-1*, the broad component being defined as its complement, $|V_{z}|>2$ km s*-1*.

We show in Figure 3 two distributions that are meant to illustrate the quality of the results: one is of the relative slope of the interpolated broad component below the peak, which is of the form $f=f_{0}+\varepsilon V_{z}$; the quantity plotted, $\varepsilon/f_{0}$, is well-behaved with an approximate Gaussian shape centred at 0.02 (0.04) km*-1* s and having a $\sigma$ of 0.06 (0.09) km*-1* s for CO(1-0) (CO(2-1)). The other distribution is of the difference between the Doppler velocity of the central peak estimated from the total spectrum in a first step and that obtained for the narrow component in a second iteration. Again, it is very well-behaved with an approximate Gaussian shape centred at 0.01 (0.04) km s*-1* with a $\sigma$ of 0.24 (0.28) km s*-1* for CO(1-0) (CO(2-1)). The separation between the narrow and broad components is indicated in Figure 2 on each of the CO(1-0) and CO(2-1) Doppler velocity spectra.

3.3 Dependence of the Doppler velocity on position angle: evidence for a tilt of the star axis

Figure 4 displays the sky maps of the integrated intensity of the broad and narrow components, the red-shifted and blue-shifted halves of the broad component being displayed separately. Figure 5 displays the dependence on position angle $\psi$ (measured counter-clockwise from the $y$ direction) of the mean Doppler velocity for each of the narrow, broad, blue-shifted broad and red-shifted broad components separately. In the case of the blue-shifted broad component we plot $<\!\!-V_{z}\!\!>$ instead of $<\!\!V_{z}\!\!>$. Fits of the form $<\!\!V_{z}\!\!>=a+b\cos(\psi-\psi_{0})$ are illustrated on the figure and their results are listed in Table 2. The values of the offset $a$ and of the phase-shift $\psi_{0}$ are negligible for the narrow component; indeed excellent fits are obtained by simply writing $<\!\!V_{z}\!\!>=b\cos\!\psi$ with $b=0.30$ (0.36) km s*-1* for CO(1-0) (CO(2-1)).

A radial expansion velocity $V$ produces a Doppler velocity $V_{z}=V\cos\!\psi\sin\!\varphi$ for the star equator and $V_{z}=V\cos\!\psi\cos\!\varphi$ for the star poles; $b$ is therefore a measure of some effective value of $V\sin\!\varphi$ for the narrow component, effective value meaning averaged over the relevant range of stellar latitudes: in particular, for $\varphi=10$∘ ($\sin\!\varphi=0.17$) $b=0.33\pm 0.03$ km s*-1* means an average expansion velocity of $\sim 0.33/0.17=1.9\pm 0.2$ km s*-1*.

3.4 Effect of a possible rotation

The fact that $b$ takes values of opposite signs for the narrow and broad components shows that we are globally dealing with expansion rather than with rotation: rotation would be about an axis making a small angle with the line of sight and projecting on the $x$ axis; but the whole $y<0$ hemisphere would be shifted the same way, blue or red, and the whole $y>0$ hemisphere would be shifted the opposite way, red or blue, for both the narrow and broad components. The sign of $b$ implies that the axis of the star makes a small positive angle $\varphi$ with the line of sight, away from the $y$ axis, namely pointing to the $\cos\!\varphi<0$ hemisphere, as can be seen directly in the right panels of Figure 4.

The small values found for the phase-shift $\psi_{0}$ show that the rotation of 20∘ applied to the system of coordinates was correct. We note however that $\psi_{0}\sim-4$∘ for the narrow component and $\sim+3$∘ for the broad component, a barely significant difference as we estimate the uncertainty on these numbers as $\sim 10$∘. Yet, the difference of $\sim 7$∘$\pm 14$∘ may reveal a small rotation of the narrow component with respect to the broad component. A rotation of mean velocity $V_{rot}$ about the star axis contributes a term $V_{rot}\sin\!\psi\sin\!\varphi$ to $V_{z}$ instead of $V\cos\!\psi\sin\!\varphi$ for expansion, adding up to $\sqrt{V^{2}+V^{2}_{rot}}\cos[\psi-\tan^{-1}(V_{rot}/V)]\sin\!\varphi$. A difference of phase shift is therefore a direct measure of $V_{rot}/V$; namely the observed phase shift, 7∘$\pm$14∘, translates into $V_{rot}/V=0.12\pm 0.26$ or, conservatively, in an upper limit of $\sim$0.38 for the ratio $|V_{rot}|/V$. To the extent that the broad component can reasonably be assumed to display no rotation, this is a direct estimate of the amount of rotation that the equatorial outflow can accommodate.

3.5 Dependence on latitude of the outflow

The value of $<|V_{z}|>$ is larger for the CO(1-0) than for the CO(2-1) broad component, 5.3 km s*-1* instead of 4.6 km s*-1*  as expected from the shape of the Doppler velocity spectra (Figure 2 right). It is a factor $\sim 15$ larger than the amplitude of the narrow component oscillation; for an inclination angle $\varphi\sim 10$∘, this implies a ratio between the average expansion velocity of the bipolar outflow and the equatorial expansion velocity of the order of $15\sin$(10∘)$\sim 2.6$. As the average expansion velocity of the equatorial outflow is $\sim 1.9$ km s*-1* for $\varphi\sim 10$∘, it means an average expansion velocity of the bipolar outflow of 2.6$\times$1.9 $\sim 4.9$ km s*-1*. Moreover, the Doppler velocity spectra extend up to $|V_{z}|\sim 11$ km s*-1*, namely a factor $\sim 2.2$ larger than the average expansion velocity of the bipolar outflow for $\varphi=10$∘, implying an important increase of the expansion velocity when moving away from the equator to the poles. These estimates are in agreement with what found by Nhung et al. (2015a).

3.6 Radial dependence of the effective emissivity

The $R$-dependence of the integrated intensity of the narrow and broad components is shown separately for CO(1-0) and CO(2-1) emissions in the left panels of Figure 6. Beyond $R=2.5$ arcsec, fits of the form $I_{0}\times 10^{-R/R_{0}}$ give results listed in Table 3. The CO(2-1) distribution is steeper than the CO(1-0) distribution: on average, the ratio between the values of $R_{0}$ obtained for CO(2-1) and CO(1-0) is 0.55 (Figure 6 right). The ratio between the values obtained for the broad and narrow components is 0.71. As $R$ is approximately equal to $r\cos\!\alpha$, where $\alpha$ is the stellar latitude, the latter ratio gives the scale of the mean angle made by the bipolar outflow with the equatorial plane: $\cos^{-1}(0.71)=45$∘ (an isotropic and uniform outflow gives 0.5, namely 60∘). At least qualitatively, this implies that the broad component is associated with an outflow reaching maximum emissivity at intermediate stellar latitudes.

The former ratio, between CO(2-1) and CO(1-0) emissions, was discussed by Nhung et al. (2018a) and interpreted as a temperature effect. The temperature-dependent factors $Q$ that scale the emissivity are of the form $Q=Q_{0}(2J+1)(2.8/T)\exp(-E_{J}/T)$; for CO(1-0) or respectively CO(2-1), $J=1$ or 2, $E_{J}=5.5$ K or 16.6 K and $Q_{0}=7.4\ 10^{-8}$ s*-1* or $7.1\ 10^{-7}$ s*-1*. This can be checked directly on the narrow component for which $r\sim R$: there is no need for de-projection. As illustrated in the right panel of Figure 6, a radial dependence of the temperature of the form $T=45\,r^{-0.72}$ K ($r$ in arcsec) reproduces the observed CO(2-1) to CO(1-0) ratio, $\sim 16\exp(-11.1/T[K])$, for $R>2.5$ arcsec.

Below 2.5 arcsec, the observed intensity of the narrow component is steeper than what a simple scaling of $R_{0}$ between CO(2-1) and CO(1-0) on one hand and between broad and narrow components on the other would imply. This translates into a wiggle on the $R$-dependence of the intensity ratio, clearly seen in the right panel of Figure 6. However, in this region, significant absorption effects may be expected and the narrow and broad components should not be studied in isolation but jointly.

4 The narrow component

4.1 Doppler velocity and intensity fluctuations

As the narrow component is associated with the equatorial outflow and is seen nearly face-on, much can be learned about its properties without having to de-project the data cube into space. The results displayed in Figure 5 and listed in Table 2 have shown that the mean Doppler velocity is well described by a $\cos\!\psi$ term with amplitude $<\!\!V\!\!>_{0}=0.30$ (0.36) km s*-1* for CO(1-0) (CO(2-1)); in order to reveal possible additional features, it is therefore pertinent to introduce, for each pixel, a quantity $\delta{V(x,y)}=<\!\!V_{z}\!\!>-<\!\!V\!\!>_{0}\cos\!\psi$. The distribution and sky maps of $\delta{V}$ are displayed in Figure 7 for CO(1-0) and CO(2-1) emissions separately. The mean and rms values of $\delta{V}$ are respectively 0.01 and 0.26 km s*-1* for CO(1-0) and 0.05 and 0.32 km s*-1* for CO(2-1). The agreement between CO(1-0) and CO(2-1) observations gives confidence in the quality of the result.

The sky maps display resolved lumps that are very similar between CO(1-0) and CO(2-1) emission, giving confidence in their reality: the kinematics of the narrow component is unambiguously shown to consist of lumps of slightly different Doppler velocities. Such lumpiness must, somehow, be associated with a related lumpiness of the morphology. In order to find out how, we define a function $\mit{\Phi}(x,y)=F(x,y)/\!\!<\!\!F\!\!>_{R}$ where $F(x,y)$ is the intensity measured in pixel $(x,y)$ and $<\!\!F\!\!>_{R}$ is its value averaged over $\psi$. In practice, we use parameterizations of $<\!\!F\!\!>_{R}$ in the form of a sum of two exponentials as given in the Appendix (Table 5), which give good fits to the observations (Figure 6). The distribution of $\mit{\Phi}$ is displayed in Figure 8 for CO(1-0) and CO(2-1) separately together with the sky maps. For both CO(1-0) and CO(2-1) it is well described by a Gaussian having mean 1 and $\sigma=0.36$. As in the case of $\delta{V}$ the maps of $\mit{\Phi}$ display features that are very similar in the CO(1-0) and CO(2-1) data, the effect of the smaller beam size being clearly apparent on the CO(2-1) map. This gives confidence in the method, showing in particular that the reference $<\!\!F\!\!>_{R}$ distributions are properly evaluated and introduce no significant bias. In both CO(1-0) and CO(2-1) emissions, an arc of enhanced intensity is seen to extend from $R\sim 5$ to $\sim 7$ arcsec in the north-eastern half of the sky and a deep depression is visible in the south-western quadrant.

To illustrate and quantify the similarity between CO(1-0) and CO(2-1) observations, we display in Figure 9 maps of

[TABLE]

Both $\delta{V_{-}}$ and $\mit{\Phi}_{-}$ are seen to take small values. The results of Gaussian fits to the distributions of $\mit{\Phi}_{\pm}$ and $\delta{V_{\pm}}$ are listed in Table 4 and illustrated in Figure 10. In both cases the difference quantities, $\mit{\Phi}_{-}$ and $\delta{V_{-}}$, are about half as wide as the averaged quantities, $\mit{\Phi}_{+}$ and $\delta{V_{+}}$. A small but significant anti-correlation is found between $|\delta{V_{+}}|$ and $\mit{\Phi}_{+}$ with $\mit{\Phi}_{+}=1.10-0.57|\delta{V_{+}}|$ ($\delta{V_{+}}$ in km s*-1*), as illustrated in the right panel of Figure 10.

The dependence on $R$ of $\delta{V}$ is displayed in the left panel of Figure 11, giving evidence for a strong modulation with maxima at $R\sim 0$, 3 and 7 arcsec for both CO(1-0) and CO(2-1) emissions. The dependence on $R$ of $\delta{V_{+}}$ is shown in the middle panel with a sine wave fit having amplitude of $\pm 0.16$ km s*-1*, a period of 3.8 arcsec and a negligible phase-shift. Note that one might fear that such an oscillation would have been absorbed in the division by $<\!\!F\!\!>_{R}$ when dealing with the intensity; this is not the case, however, as $<\!\!F\!\!>_{R}$, being defined as shown in Figure 6 and Table 5, is free of such an oscillation.

4.2 Spirals and arcs

The features displayed on the sky map of $\mit{\Phi}_{+}$ are reminiscent of the spiral studied by Homan et al. (2018) and the rings reported by Winters et al. (2007). The question then arises of the relation between these and of their relation to the radial oscillations displayed by $\delta{V_{+}}$. In order to quantify the spiral assignment, we parameterize a possible spiral as $R_{sp}=R_{0}+\frac{1}{2}h\psi/\pi$ and use as a model $\mit{\Phi}_{+}(R,\psi)=1+\mit{\Phi}_{0}\cos[2\pi(R-R_{sp})/h]$. The best fit is obtained for $\mit{\Phi}_{0}=0.13$, $R_{0}=0.87$ arcsec and $h=4.0$ arcsec. Similarly, we define circles of radius $R_{circle}=1.9$ (modulo 3.8) arcsec meant to match the zeroes of $\delta{V_{+}}$ and model $\delta{V_{+}}$ as $\delta{V_{+}}=0.16\sin[2\pi(R-R_{circle})/3.8]$. We restrict the analysis to projected distances from the star, $R$, in excess of 1.5 arcsec, where the spiral and circular assignments are more reliable. Figure 12 displays the sky maps obtained from the resulting models (a spiral for $\mit{\Phi}_{+}$ and circles for $\delta{V_{+}}$) together with the dependence of $\mit{\Phi}_{+}$ and $\delta{V_{+}}$ on both $(R-R_{sp})/h$ and $(R-R_{circle})/3.8$ arcsec.

The result is clear: the spiral model fits $\mit{\Phi}_{+}$ but does not fit $\delta{V_{+}}$, the circular model fits $\delta{V_{+}}$ but does not fit $\mit{\Phi}_{+}$. However, the best-fit spiral is but a crude representation of the observations and differs from that observed by Homan et al. (2018): it has a pitch of 4.0 arcsec instead of 2.8 arcsec (note that Homan et al. erroneously quote an inter-arm distance of only $\sim 1$ arcsec in their paper); the reason is that it gives a better fit to the observations over the whole circle $R<8$ arcsec while the spiral described by Homan et al. covers only $R<5$ arcsec. This underlines the arbitrariness inherent to such spiral assignments when the observed features are insufficiently clear.

As the claim for evidence for a spiral enhancement is the central argument developed by Homan et al., we illustrate this point further in Figure 13 that displays sky maps of the intensity integrated over $|V_{z}|<0.3$ km s*-1*, the central velocity bin used by these authors in their Figure 5. Both CO(1-0) and CO(2-1) maps are shown, the spiral appearance being clear in both cases. The small width of the velocity interval used in Figure 13 is equal to the amplitude of the oscillation of the Doppler velocity as a function of position angle $\psi$; as this may be a cause of bias, we repeat the exercise by subtracting from the Doppler velocity a correction term of $0.3\cos\!\psi$ km s*-1* before selecting $|V_{z}|<V_{cut}$, where $V_{cut}$ takes values of 0.1, 0.3, 0.5 and 0.7 km s*-1*. The result is displayed in Figure 14: the effect of the correction is negligible. Finally, Figure 15 displays maps of the intensity integrated over the Doppler velocity interval $|V_{z}|<0.3$ km s*-1* for the broad and narrow components separately.

Bringing together the information contained in Figures 13 to 15, a spiral-like enhancement is visible on the narrow component for both CO(1-0) and CO(2-1) emissions but the broad component displays no significant feature. Correcting for inclination has no significant effect and the narrow component enhancement is sharper for narrower velocity intervals. Radial enhancements are seen to link the spiral arms at irregular angular intervals of typically 50∘ to 60∘; they may equally well be interpreted as successive enhancements of the intensity along the spiral arms. Different choices of colour scale produce different images revealing enhancements of apparently different shapes. As commented by Homan et al. (2018), the spiral-like enhancement is far from an ideal Archimedes’ spiral, at strong variance with observations made, for example, on carbon star LL Pegasi (Morris et al., 2006).

Finally, we use the CO(1-0) and CO(2-1) maps of Figure 13 to draw a curve, $S$, that runs along maximal intensity, without attempting to find its equation, and compare it with the maps of $\mit{\Phi}_{+}$ and $\delta{V_{+}}$ displayed in Figure 9. Pixels having their centre at less than 0.7 arcsec from $S$ are delineated in Figure 16; they have average values of $\mit{\Phi}_{+}$ and $\delta{V_{+}}$ of respectively 1.14 and $-$0.02 km s*-1*. As $\mit{\Phi}_{+}$ is the intensity divided by its average over position angle at the relevant value of $R$, it tends to be enhanced at slightly larger values of $R$ than $S$ is; but, apart from this difference, $S$ matches reasonably well the enhanced arc of $\mit{\Phi}_{+}$ emission. On the contrary, $S$ spans a broad range of values of $\delta{V_{+}}$, as could be expected from the circular pattern displayed by this quantity.

4.3 Line widths and flaring of the equatorial outflow

Having evaluated the average dependence of $<\!\!V_{z}\!\!>$ on $R$ and $\psi$ we correct the Doppler velocity spectra accordingly in each pixel, using a correction term of the form $\Delta{V_{z}}=a\cos\!\psi+b\sin(2\pi R/3.8)$. We find that $a$ and $b$ can be set to 0.30$\pm$0.05 and $-0.20\pm 0.05$ km s*-1* for both CO(1-0) and CO(2-1) emissions to obtain minimal line widths, in agreement with the results obtained above. Fitting Gaussians to the line profiles of CO(1-0) and CO(2-1) emissions gives $\sigma$ values of respectively 0.62 and 0.67 km s*-1* before correction, decreasing to respectively 0.51 and 0.55 km s*-1* after correction: the correction is found to significantly decrease the line width, by $\sim 18$%. If, instead of applying the correction we re-centre each Doppler velocity spectrum on its mean value in each pixel separately, we obtain $\sigma$ values of respectively 0.44 and 0.47 km s*-1*, showing that after correction there is not much room left for an additional systematic contribution to the line width. Indeed, as a function of position angle $\psi$, the line widths fluctuate by less than 0.03 km s*-1* with respect to their mean. Table 6 in the Appendix summarizes the results.

The corrected line widths receive contributions from several factors, instrumental resolution, thermal broadening, turbulence and flaring being expected to be the most important. Instrumental resolution ($\sigma$) contributes $\sim 0.1$ km s*-1* and thermal broadening a factor $\sqrt{2kT/M_{CO}}$ where $T$ is the temperature, $k$ Boltzman constant and $M_{CO}$ the mass of the CO molecule, meaning 0.16 km s*-1* at $T=45$ K. Together, they contribute therefore $\sim 0.2$ km s*-1*, leaving some 0.4 km s*-1* to be accounted for by turbulence and flaring, the uncertainty on this number not exceeding 0.1 km s*-1*. We retain therefore the value of $0.4\pm 0.1$ km s*-1* for the ($\psi$-independent) contribution of flaring and turbulence. Defining flaring (FWHM) as $a_{flar}=2.35\times\sigma_{flar}=2.35\times\mbox{Rms}(\sin\!\alpha)$ and recalling that $V\sin\!\varphi\sim 0.33$ km s*-1*, the flaring contribution to the line width is, to leading order in $\varphi$, $(0.33/\sin\!\varphi)\,\sigma_{flar}$. Blaming on it the whole $0.4\pm 0.1$ km s*-1* means therefore $\sigma_{flar}=(1.2\pm 0.3)\sin\!\varphi$. For $\varphi$=10∘, $\sigma_{flar}=(0.20\pm 0.05)$ giving an upper limit of $\sim 0.6$ on $a_{flar}$ corresponding to a latitude of $\pm 17$∘. Note that, to leading order in $\varphi$, expansion alone contributes to the broadening due to flaring, the rotation velocity being always normal to the star axis. Note also that if effects such as turbulence contribute significantly to the line width, the flaring contribution and the limit obtained on the flaring angle will accordingly be simply smaller.

The line widths measured on and in between the spiral arms, using only the $\cos\!\psi$ term as correction and requiring $R$ to exceed 1.5 arcsec, differ by $\sim 0.06$ km s*-1* as do the line widths measured on and in between the concentric circles (see Table 6): large $|\delta{V_{+}}|$ values (corresponding to $|R-R_{circle}|>0.95$), or large $\mit{\Phi}_{+}$ values (corresponding to $|R-R_{sp}|<h/4$) tend to produce narrower line profiles. This barely significant difference, if blamed on turbulence, would mean a contribution of $\sim 0.24$ km s*-1*. Line widths ($\sigma$) measured at less than 0.7 arcsec distance from the spiral-like curve $S$ are 0.40 km s*-1* and 0.42 km s*-1* for CO(1-0) and CO(2-1) emissions respectively, namely $\sim 0.05$ km s*-1* lower than average, again a barely significant difference.

4.4 Summary

In summary, the picture of the narrow component that emerges from this analysis is of an equatorial outflow having a small inclination $\varphi$ with respect to the sky plane. On average, it expands radially with a uniform velocity of $\sim(0.33\pm 0.03)/\sin\!\varphi$ km s*-1*, namely 1.9$\pm$0.2 km s*-1* for $\varphi$=10∘. An upper limit of 0.38 has been obtained for the ratio of a possible rotation velocity to the radial expansion velocity. Lumps of small velocities normal to the sky plane have been detected in excess to the contribution of the expansion velocity; they are at the scale of 0.24 km s*-1* (1$\sigma$) and display a clear modulation with $R$ of amplitude $\pm$0.16 km s*-1* and period 3.8 arcsec. The intensity displays fluctuations at the level of $\pm$36%, however not clearly correlated with the velocity fluctuations, the latter following an approximate concentric circular pattern and the former being dominated by an arc of enhanced intensity in the north-eastern half of the sky and a deep depression in the south-western quadrant, in qualitative agreement with the spiral-like feature described by Homan et al. (2018). A small but significant anti-correlation is detected between the amplitude of the intensity fluctuations and the amplitude of the small velocity components detected in excess to the contribution of the expansion velocity. Line widths of $\sim$1.2 km s*-1* (FWHM) have been measured, after subtraction of the effect of the inclination of the star equator with respect to the plane of the sky, consistent with a very small contribution of other effects such as flaring, instrumental resolution and thermal broadening. The dependence on $|\cos\!\psi|$ of the flaring contribution to the line width has provided an estimate of an upper limit on the flaring angle, $\pm$17∘ at half maximum for $\varphi$=10∘. The study has failed to reveal new features that could strengthen the evidence for the significance of the spiral observed in the low $V_{z}$ channel maps. On the contrary it has revealed as significant a circular feature in the maps of the mean Doppler velocity, apparently uncorrelated with the spiral. As both are near the limit of sensivity of the observations, caution must be exercised when proposing an interpretation.

5 The broad component

Contrary to the narrow component, the broad component, which is seen pole-on, needs to be de-projected in space to reveal the details of its morphology and kinematics. To do so, we normally use a radial expansion velocity of the form $V=V_{0}(1-\lambda\cos\!2\alpha)$, where $\alpha$ is the stellar latitude, to de-project the effective emissivity following the method developed in Nhung et al. (2018a, b). We use a polar expansion velocity $V_{pole}=V_{0}(1+\lambda)=10.3$ (10.7) km s*-1* for the blue(red)-shifted part of the spectrum (see Section 5.3 below) and a prolateness parameter $\lambda=0.7$, meaning an equatorial velocity $V_{eq}=V_{0}(1-\lambda)=1.9$ km s*-1*. The inclination of the star axis with respect to the line of sight is fixed as $\varphi$=10∘.

The low value adopted for the equatorial radial velocity deserves some comments. It is obtained from the arguments developed in Section 3.3 that measure the radial wind velocity at equator as (0.33$\pm$0.03)/$\sin\varphi$ km s*-1*, and from the choice of 10∘ for the value of $\varphi$. The escape velocity from a solar mass star is 2.7 km s*-1* at $r=250$ au (where broad and narrow components can be described separately, see below) and reaches 1.9 km s*-1* at $r=490$ au, meaning 5.6 arcsec. To compensate gravity at $r=250$ au requires a rotation velocity 2.6 times larger than the maximum acceptable according to the arguments developed in Section 3.4. It is therefore reasonable to conclude that the equatorial wind velocity obtained in Section 3 is below escape velocity over an important fraction of the radial range explored here. In the region where the wind is being accelerated, the competition between the mechanism of acceleration and the gravity of the star governs the gas dynamics. It is again reasonable to expect that the terminal velocity has been essentially reached at $r\sim 250$ au. A realistic model of the wind velocity would need to take all these arguments in due account, in particular the effect of gravity causing a decrease of the radial velocity at the equator. However, a reliable treatment requires an understanding of the dynamics governing small distances to the star, where other tracers than CO molecules are necessary, which we leave for a later publication including analyses of SiO and SO2 emissions (Tuan-Anh et al. , in preparation). In the present work, the form of the radial wind velocity used for de-projection is meant to simply help with the evaluation of the effective emissivity and ultimately the flux of matter, which it does well (see Section 6.3).

5.1 Broad-narrow separation

We remark that our definition of the broad component includes part of the equatorial outflow, obtained by interpolation of the blue- and red-shifted wings of the Doppler velocity spectrum below the narrow central peak. As a result, there exists a region in space where both components merge and where the distinction between them becomes meaningless. This is illustrated in Figure 17 that displays the dependence on $R$ of the de-projected emissivity of the broad component near the equator, averaged over $z^{\prime}=\pm 0.5$ arcsec (measured along the star axis), and compares it with the intensity $F_{narrow}$ measured for the narrow component. It shows dominance of the narrow component for values of $R$ exceeding typically 1 arcsec. We recall that for $\varphi$=10∘ $\sigma_{flar}\sim 0.2$, implying that such a $z^{\prime}$ interval (FWHM) is covered at a distance $R\sim 2.1$ arcsec. In practice, we shall retain 2 arcsec as the distance from the star where the two components merge and where the separation between them becomes meaningless.

5.2 General properties

The de-projected emissivity, $\rho$, multiplied by $r^{2}$, is displayed in Figure 18 (left) for CO(1-0) and CO(2-1) emissions separately as a function of the distance to the star, $r$, together with polynomial fits, $F_{1}(r)$ (Table 7 in Appendix). The functions $\rho r^{2}/F_{1}(r)$ are displayed in Figure 18 (middle) as a function of stellar latitude $\alpha$, together with polynomial fits $F_{2}(\sin\!\alpha)$ (Table 7). They increase from the equator up to latitudes of the order of $\pm$45∘ where they reach a maximum before decreasing smoothly toward the poles. The reliability of this result is critically discussed in Section 5.3 below.

The normalised effective emissivity, $\rho^{*}(x,y)=\rho(x,y)r^{2}/[F_{1}(r)F_{2}(\sin\!\alpha)]$ is displayed in Figure 18 (right) as a function of stellar longitude, $\omega$. It is nearly constant for both CO(1-0) and CO(2-1) emissions, providing evidence for axi-symmetry when averaged over other variables. The region $\omega>180^{\circ}$ is in very slight excess with respect to the region $\omega<180^{\circ}$. In each case, $\rho^{*}$ is averaged over the other two variables with the conditions $0.5<r<8$ arcsec and $|\sin\!\alpha|<0.9$.

5.3 Dependence on latitude and polar depression

It is important to assess the reliability of the evidence for a polar depression of the effective emissivity presented in the preceding section. As was remarked by Nhung et al. (2018a) the shape of the Doppler velocity spectrum near its end points is closely related to the latitudinal dependence of the effective emissivity near the stellar poles. In particular, a locally isotropic and $r$-independent radial wind associated with a spherical emissivity proportional to $1/r^{2}$ produces a Doppler velocity spectrum proportional to $[1-(V_{z}/V)^{2}]^{-1/2}$ when integrated over a circle on the sky plane. This means sharp horns at the ends of the spectrum, which are absent from the present observations: the CO(2-1) spectrum displays no horn and the CO(1-0) spectrum displays horns much less sharp than produced by the above law. The small inclination of the star axis with respect to the line of sight can be ignored in this context. A natural interpretation is a depression of the effective emissivity near the stellar poles, which may be associated with a decrease of density and/or an effect of temperature and/or an effect of absorption. However, one might argue that this result is an artefact of the form of wind velocity used for de-projection: we address this issue below.

In addition to a marked difference between the CO(1-0) and CO(2-1) Doppler velocity distributions (Figure 2 right), a notable feature of the broad component is the asymmetry between the red-shifted and blue-shifted brightness, which had been previously noted by Nhung et al. (2015a). It is illustrated in Figure 19 for CO(1-0) and CO(2-1) emissions separately, together with the ratio of the respective spectra. The latter displays a remarkably smooth shape if one excludes end effects where noise is important. A Gaussian fit in the interval $-9<V_{z}<9$ km s*-1* gives a mean at 0.5 km s*-1* and a $\sigma$ of 7.3 km s*-1*. We evaluate the position of the blue end-point at $\sim-10.2$ km s*-1* and that of the red end-point at $\sim 10.6$ km s*-1*. The precise location of the end points of the Doppler velocity spectrum fixes the polar wind velocity just above their values, but not higher. Using too large a polar velocity, larger than the end point velocity by a fraction $\varepsilon$, causes the de-projected effective emissivity to cancel at latitudes exceeding $\sin^{-1}(1-\varepsilon)\sim\pi/2-\sqrt{2\varepsilon}$: $\epsilon=5$% means a cut-off latitude of 72∘.

Another notable property of the Doppler velocity spectra concerns their evolution when moving away from the star. While being approximately independent of position angle, they narrow down when $R$ increases, as illustrated in Figure 20. This is not surprising, larger values of $R$ probing regions of the outflow farther away from its axis. More surprising, however, is the appearance or enhancement of horns at the ends of the spectra in the $R$ intervals between 1 and 3 arcsec. A spherical emissivity distribution would produce instead horns that decrease in amplitude and sharpness when $R$ increases: this, again, suggests that the emissivity is depressed in the polar regions. To further illustrate the evolution of the spectrum horns with $R$ we define a parameter $\xi$, meant to describe how sharp and enhanced they are, as the ratio between the content of the spectrum in the 2 km s*-1* wide interval below the end point and in the interval covering between 5 and 3 km s*-1* below the end point. The dependence of $\xi$ on $R$ is displayed in the right panel of Figure 20 for CO(1-0) emission. The evidence for peaking just below $R\sim 2$ arcsec is very strong. The effect is also present but weaker for the CO(2-1) spectrum. The ratio between the CO(2-1) and CO(1-0) spectra remains well behaved in separate $R$ intervals and evolves from a nearly flat distribution near the star to a distribution reaching a maximum at low $|V_{z}|$ for large values of $R$.

As a further illustration of the evidence for polar depression of the emissivity, we show in Figure 21 the dependence on stellar latitude $\alpha$ of the broad component de-projected emissivity for wind velocity distributions of different width around the poles. The more collimated is the wind velocity, the closer to the poles is the maximum of the emissivity, but the polar depression is always present.

5.4 Inhomogeneities

The distributions of $\rho^{*}$ are shown in the left panel of Figure 22, together with Gaussian fits having (mean, rms) values of (0.98, 0.21) for CO(1-0) and (0.95, 0.30) for CO(2-1). Averaging between CO(1-0) and CO(2-1), the lumpiness is even smaller than for the narrow component: $\sigma=0.26$ instead of 0.36. In order to quantify this statement, we also show in Figure 22 the distributions of the quantities $\rho^{*}_{+}=\frac{1}{2}(\rho^{*}_{\rm{CO(2-1)}}+\rho^{*}_{\rm{CO(1-0)}})$ and $\rho^{*}_{-}=\frac{1}{2}(\rho^{*}_{\rm{CO(2-1)}}-\rho^{*}_{\rm{CO(1-0)}})$. Gaussian fits give (mean, $\sigma$)=(0.97, 0.24) for $\rho^{*}_{+}$ and (0, 0.14) for $\rho^{*}_{-}$: the dispersion of $\rho^{*}_{+}$ owes as much to the difference between CO(1-0) and CO(2-1) data as to its intrinsic dispersion, $\sim 0.19=\sqrt{0.24^{2}-0.14^{2}}$. At variance with the narrow component, the broad component displays more global uniformity and less similarity between CO(1-0) and CO(2-1) data, preventing a reliable identification of significant inhomogeneity of the $\rho^{*}_{+}$ distribution.

5.5 Summary

In summary, the broad component has been shown to be physically distinct from the narrow component at distances from the star in excess of 2 arcsec or so. At distances smaller than 1 arcsec, both components are undistinguishable and must be studied jointly. At variance with the narrow component, the broad component displays significant differences between CO(1-0) and CO(2-1) emissions, and smaller intensity inhomogeneity. Moreover, it shows a small but significant asymmetry between its blue-shifted and red-shifted parts. The effective emissivity is axi-symmetric. The wind velocity, assumed to be small near the equator in order to join smoothly with the equatorial outflow, must increase toward the poles in order to produce Doppler velocities reaching $\sim 11$ km s*-1*. The absence of strong horns at the extremities of the measured Doppler velocity spectra generally implies significant depressions of the effective emissivity near the stellar poles, causing it to reach maximum at intermediate latitudes of the order of $\sim 45$∘ to 60∘.

6 Joint description

6.1 Aim and method

The aim of the present section is to draw a unified picture of the morphology and kinematics of the circumstellar envelope, as realistic and reliable as possible. The preceding sections have provided detailed descriptions of its main components, a narrow component associated with an equatorial outflow in slow expansion and a broad component associated with a bipolar, biconical outflow in more rapid expansion. However, the effect of temperature and absorption has been ignored; accounting for it requires a radiative transfer calculation and a joint description of the morpho-kinematics of the central region, $r<\sim 2$ arcsec, where the narrow and broad components merge.

The two-component description has been very effective in clarifying the properties of the morpho-kinematics of the circumstellar envelope. At distances from the star in excess of $\sim 2$ arcsec, it has been shown to correspond to two different regions of space, each having its own identity. However, this does not imply that both components are separated by a gap; on the contrary, we shall see that they join smoothly. Moreover, it is important to be conscious of the arbitrariness inherent to the de-projection of radio observations, which has been discussed in detail by Diep et al. (2016) and Nhung et al. (2018b). In particular, once a specific form has been assumed for the field of wind velocities, effective emissivities that reproduce exactly the observed data cubes can be easily de-projected: what matters then is to identify the main features that can be reliably ascertained and to discuss critically their reality.

In the remainder of the present section, we first use the two-component description to discuss the effects of temperature, evaluate its distribution in space and calculate the flux of matter and the mass-loss rate after de-projection of the broad component; we then use a joint description to evaluate the effect of absorption, and, taking it in due account, we calculate the distribution of temperature, density and flux of matter in stellar coordinates; finally, we discuss the results obtained.

6.2 Temperature

Ignoring absorption, the CO(2-1) and CO(1-0) brightnesses are expected to be in a ratio $Q_{2-1}/Q_{1-0}=16\exp(-11.1/T[K])$ where 16 is the ratio between the respective $Q_{0}(2J+1)$ factors and 11.1 K the difference between the respective $E_{J}$ temperatures, as discussed in Section 3.6. Namely, ignoring absorption, the brightness ratio $R_{T}=f_{\rm{CO(2-1)}}/f_{\rm{CO(1-0)}}$ provides a measure of the temperature: $T[{\rm K}]=-11.1/\ln(R_{T}/16)$ in each data cube element separately. The dependence of $T$, calculated in pixels of 0.9$\times$0.9 arcsec2 and Doppler velocity bins of 0.6 km s*-1*, on $R$, $\psi$ and $V_{z}$ is displayed in Figure 23. It is approximately independent of $\psi$ and takes values mostly below 20 K. It decreases with $R$ from $\sim 30$ K at $R=1$ arcsec to $\sim 10$ K at $R=8$ arcsec. It is higher on the broad component than on the narrow component, reaching maxima at intermediate values of $|V_{z}|$. In what follows, we ignore the 9 central pixels having $R<0.5$ arcsec and require as usual $R$ not to exceed 8 arcsec.

As a check of the method, we correct the CO(1-0) and CO(2-1) data cubes for temperature: we divide each of them by the associated $Q$ factor calculated using the temperature measured as described in the preceding paragraph. The result is displayed in the upper left panel of Figure 24. A small difference between the CO(1-0) and CO(2-1) corrected spectra can be seen at the red-shifted end; it is due to uncertainties associated with the temperature evaluation near the end points.

The temperature-corrected data cubes are then de-projected in space using an axi-symmetric, $r$-independent radial velocity of the form $V=V_{eq}+(V_{pole}-V_{eq})\sin^{2}\!\alpha$ with $V_{eq}$ and $V_{pole}$ as defined in Section 5.1: $V_{eq}=1.9$ km s*-1* and $V_{pole}=10.3(10.7)$ km s*-1* for the blue(red)-shifted part of the spectrum. De-projection measures now the density rather than the effective emissivity. It is used to average the temperature as a function of $r$ and $\alpha$ instead of $R$ and $V_{z}$. The result is illustrated in Figure 24, showing that a good fit is obtained using a form $T=40(1-0.25\cos\!4\alpha)\exp(-r/5.52)$ with $T$ in Kelvin and $r$ in arcsec (an exponential gives a better fit than a power law). The relative uncertainty is $\sim 13$% as obtained from the $\sigma$ of a Gaussian fit to the residuals.

6.3 Flux of matter and mass-loss rate

The data cube, averaged between CO(1-0) and CO(2-1) emissions after correction for temperature as described in the preceding section, is now split into its narrow and broad components following the method described in Section 3.2. We use the results obtained in Sections 4 and 5 to construct as simple as possible a model of their morphology and kinematics in the region where they are distinct, $R>2$ arcsec. For both components, the wind velocity is taken of the form used in the preceding section (independence from stellar longitude $\omega$ is assumed throughout) and the inclination angle $\varphi$ is fixed at 10∘.

From the temperature corrected intensity of the narrow component, assuming its density to be a Gaussian in $\sin\!\alpha$ of $\sigma=\sigma_{flar}=0.20$ (see Section 4.3), we obtain its dependence on $R$ or, to a good approximation, on $r$. Its product by $r^{2}$ is displayed in the left panel of Figure 25. It has an average value of 29 CO molecules cm*-3*arcsec2 and shows that the equatorial outflow has been developing at approximately constant mass-loss rate during the past 2 to 3 thousand years. However, it displays significant deviations with respect to such an approximation, at the typical level of $\pm 10$%, with maxima at $\sim 3.5$ and 5.2 arcsec and minima at $\sim 2$ and 8 arcsec.

The dependence on $r$ and $\alpha$ of the broad component de-projected density multiplied by $r^{2}$ is displayed in the middle panels of Figure 25. It is approximately constant for $r>2$ arcsec, at a level of 22 CO molecules cm*-3*arcsec2, showing that the bipolar outflow, as in the case of the equatorial outflow, has been developing at approximately constant mass-loss rate during the past 2 to 3 thousand years. Its dependence on $\alpha$ displays maxima at intermediate latitudes.

The dependence on stellar latitude $\alpha$ of the flux of matter $\rm d\mit{\dot{M}}/\rm d\mit{\Omega}$ per solid angle $\rm d\mit{\Omega}=2\pi\cos\!\alpha\,\rm d\alpha$ can be obtained from the product of the density by $r^{2}$ and by $2\pi V\cos\!\alpha$ where $V$ is the wind radial velocity. The result is displayed in the right panel of Figure 25 for the broad and narrow components separately and for their sum. Assuming a solar-like abundance ratio, $\rm{CO/H}\sim 2.5\,10^{-4}$, the total mass-loss rate is measured this way to be $0.9\,10^{-7}$ M*⊙yr-1*, with the contribution of the narrow component being only 10% of the total. This result is independent of what happens at short distances from the star and needs only to be corrected for absorption.

In Section 5, we commented on the low value of the equatorial wind velocity and on its use for de-projection and for the evaluation of the effective emissivity and flux of matter; we were then claiming that concerns related to the equatorial wind velocity being below escape velocity over a significant fraction of the explored radial range were irrelevant in the present context. To illustrate this claim we display in panels b) and c) of Figure 25 the results of a de-projection using as $V_{eq}$ the escape velocity evaluated for a solar mass star, $V_{esc}$ [km s*-1*]$=42/\sqrt{r[\rm au]}$. The difference with the $V_{eq}$ result is indeed very small and does not significantly affect any of the conclusions made in the present Section.

6.4 Absorption

In all previous considerations, absorption was ignored. It was implicitly assumed that its effect would be small, would not vary too much across the data cube and would not differ too much between CO(1-0) and CO(2-1) emissions. We are now in a position to evaluate it for each line and in each pixel independently. In each pixel, the relation between $V_{z}$ and $z$ is defined by choosing the standard wind velocity distribution ($V_{eq}=1.9$ km s*-1* and $V_{pole}=10.3$ (10.7) km s*-1* for the blue(red)-shifted part of the spectrum); the temperature $T$ and the density $d$ are then well defined as functions of $z$, or $\sin\!\alpha=z^{\prime}/r$ or $V_{z}$. In order to avoid too low signal-to-noise values, we use large data cube elements: 0.9 arcsec wide in $x$ and $y$ and 0.6 km s*-1* wide in Doppler velocity. We moreover assume that the wind reaches terminal velocity at $r=0.5$ arcsec and increases in proportion to $r$ below this value; in practice, this is largely irrelevant to the arguments of the present article, which do not depend significantly on the morpho-kinematics in the close environment of the star.

We start from the results of Section 6.2, including de-projected effective emissivity for each line and an estimate of the temperature in the limit of zero absorption (Figures 23 and 24). We then use the de-projected effective emissivity, $\rho$, to calculate the brightness, this time accounting for absorption (ignoring absorption gives the same values as we started from). This is done in each pixel and for each line separately by scanning along the line of sight in small $z$ steps of width $\Delta{z}=0.004$ arcsec, starting from the red-shifted end ($z=30$ arcsec). At each step, we calculate $V_{z}$ and add to the corresponding bin of the brightness distribution the contribution of the new $z$ step, $\rho\Delta z$, while subtracting the attenuation of the contributions of the earlier $z$ steps by a factor $f_{abs}\rho\Delta z$ where $f_{abs}=F_{J}[\exp(\lambda_{J})-1]$; here, $F_{J}$ is a constant inversely proportional to the third power of the line frequency and $\lambda_{J}$ depends on temperature as $\lambda_{J}=2E_{J}/[T(J+1)]$. We obtain this way a brightness $f^{\prime}=Af$, smaller than the measured brightness $f$ by a factor $A$ evaluated in each data cube element separately. Its distribution in the $R$ vs $V_{z}$ plane and its projection on the $R$ and $V_{z}$ axes are displayed in Figure 26 for each line separately. Absorption is measured at typical levels of 10 to 30%. It is larger for the narrow than for the broad component because in the former case $V_{z}$ varies little when $z$ scans along the line of sight and line emission keeps self-absorbing. The difference between CO(1-0) and CO(2-1) absorptions is small, it rarely exceeds 10%.

6.5 Results

Having evaluated the effect of absorption and found it relatively small and well-behaved, we obtain new evaluations of the temperature and density using absorption-corrected flux densities, $f^{*}=f/A$. To avoid introducing undesired fluctuations, we use as $A$ the product of its projected distributions on the $R$ and $V_{z}$ axes, as displayed in Figure 26, properly normalized. This is equal to the value of $A$ evaluated in each data cube separately to within 6% (rms). The temperatures $T^{*}$ evaluated from the ratio $R^{*}_{T}$ between the CO(2-1) and CO(1-0) values of $f^{*}$ differ from the former values $T$ to the extent that the absorptions of CO(1-0) and CO(2-1) emissions are different. The resulting distributions of temperature and density, $T^{*}$ and $d^{*}$, averaged over stellar longitude $\omega$, are displayed in the star meridian plane, $z^{\prime}$ vs $R^{\prime}$, in Figures 27 and 28; also shown are their projections on the $r$ and $\alpha$ axes. By construction, they describe equally well CO(1-0) and CO(2-1) observations. Absorption is small enough to make further iteration unnecessary: using the $T^{*}$ and $d^{*}$ distributions displayed in Figures 27 and 28 to perform a radiative transfer calculation, we reproduce the measured data cubes to excellent precision.

The temperature (Figure 27), averaged over stellar latitude, decreases with $r$ nearly exponentially, by a factor of $\sim 2$ over $\sim 2$ arcsec, more precisely as $T^{*}\sim 109\exp(-r/3.1)$ with $T^{*}$ in Kelvins and $r$ in arcsec; the power law best fit gives an exponent equal to unity, namely $T^{*}\sim 106/r$; averaged over $r$, the temperature displays maxima at intermediate stellar latitudes, being slightly warmer in the red-shifted than in the blue-shifted hemisphere.

The density (Figure 28) has been evaluated for three values of the inclination angle $\varphi$ of the star axis with respect to the line of sight, changing $V_{eq}$ accordingly: $V_{eq}=(0.33\pm 0.03)/\sin\!\varphi$ km s*-1*. It displays an excess of $\sim 30$% at distances from the star smaller than $\sim 2$ arcsec and then decreases slightly faster than $1/r^{2}$. Such a decrease is expected as the result of dissociation of the CO molecules by interstellar UV radiation. However, the present data cannot measure its value precisely. As was illustrated in Figure 25 b, it depends on the form used to describe the equatorial wind velocity, which is arbitrary as long as it averages to some 2 km s*-1* over the explored radial range. We note that estimates of the decrease have been given by Mamon et al. (1988) and Groenewegen (2017); while the former predicts a decrease of $\sim$9% at $r=5$ arcsec and $\sim$33% at $r=10$ arcsec, the latter predicts instead a negligible effect: 1% at 5 arcsec and 8% at 10 arcsec. The dependence on $\alpha$ displays maxima at intermediate latitudes but the value at the equator depends on $\varphi$: the smaller $\varphi$ the larger the enhancement of the equatorial density, trivially scaling with $V_{eq}$ (as $\rho/f=\mbox{d}V_{z}/\mbox{d}z$ scales with $V_{eq}$ at the equator).

The mass-loss rate, which was evaluated as $0.9\,10^{-7}$ M⊙ yr*-1* in Section 6.3 ignoring absorption, becomes ($1.6\pm 0.4$) $10^{-7}$ M⊙ yr*-1* when using the absorption corrected values of the temperature and density, $T^{*}$ and $d^{*}$. The uncertainty is estimated by including effects of absorption and inclination angle. The dependence of the flux of matter on stellar latitude $\alpha$ is displayed in Figure 29.

Several uncertainties are attached to these results. In addition to those inherent to the arbitrariness of the form adopted for the wind velocity, the evaluation of the temperature becomes less reliable at higher values. To illustrate this point, we display in Figure 30 the dependence on $R_{T}$ of temperature dependent quantities of relevance to the radiative transfer calculation: the temperature $T$, the absorption coefficients, $f_{abs(1-0)}$ and $f_{abs(2-1)}$ and the factors $Q_{(1-0)}$ and $Q_{(2-1)}$ that control the temperature dependence of the emissivity. For convenience, we normalize all quantities to their values at $R_{T}=1$; the steep dependence on $R_{T}$ near the ends of the allowed range (0 to 16) implies large uncertainties on the evaluation of the temperature: for $T=30$ K, $R_{T}=11$ but $T$ increases to 40 K when $R_{T}$ increases by only 10%: large temperatures are therefore affected by very large uncertainties. This has a strong effect on the evaluation of the emission probability but relatively little effect on the evaluation of the absorption because the $f_{abs}$ coefficients become very small at high temperatures, the exponential factor $\exp(E_{J}/T)$ acting in opposite directions in the expressions for $Q$ and $f_{abs}$.

The choice of a radial wind velocity independent from $r$, while sensible, is arbitrary. Deviation from a radial form is best discussed in terms of a possible rotation component, which was done in Section 3.4. A possible dependence of the radial velocity on $r$ would obviously affect the results. However, we have seen (Figures 25 and 28) that the radial fluctuations of the mass-loss rate and of the density calculated assuming an $r$-independent wind velocity are small. As wind velocity and mass-loss rate are linearly related, fluctuations of one must be expected to be accompanied by fluctuations of the other. Namely, one might use for de-projection a wind velocity distribution that would display small fluctuations in $r$ mimicking those observed on the mass-loss rate; one would then obtain a new density that would be just as reliable as that displayed in Figure 28. Indeed, both wind velocity and density are likely to display related fluctuations in $r$; the important point is that both are expected to be small. The density enhancement observed at $r<\sim 2$ arcsec is independent of the choice of inclination angle $\varphi$ and would imply a positive wind velocity gradient in this region if the flux of matter were to be constant.

Taking all above arguments into account, we estimate that, on average, temperatures are measured with a precision of $\sim\pm 25$% and densities with a precision of $\sim\pm 15$%.

7 Summary and discussion

7.1 Summary

The present analysis of ALMA observations has confirmed the overall picture that had been previously drawn of the circumstellar envelope of EP Aqr, displaying approximate axi-symmetry about an axis making an angle of $\sim 10$∘ with the line of sight and projecting on the sky plane some 20∘ west of north. The Doppler velocity spectra have been separated into two components that have been shown to be associated with distinct entities: a narrow central component associated with an equatorial region expanding at low velocity $V_{eq}\sim(0.33\pm 0.03)/\sin\!\varphi$ km s*-1* and a broad component associated with a bipolar outflow expanding at a velocity increasing with stellar latitude $\alpha$ up to 10 to 11 km s*-1* near the poles. The presence of faint wings at approximate noise level possibly extending up to some $\pm 16$ km s*-1*  as suggested by the observation of SiO emission by Homan et al. (2018), cannot be excluded. An upper limit of 38% of the expansion velocity has been placed on a possible rotation velocity of the equatorial outflow. Both components merge at an angular distance of $\sim 2$ arcsec from the star ($\sim 200$-250 au). A joint description of both the narrow and broad components, for both CO(1-0) and CO(2-1) emissions, has been presented, taking temperature and absorption effects in due account. The values taken by the temperature and density have been evaluated in the star meridian plane and the uncertainties attached to these have been discussed.

The flaring of the equatorial outflow was shown to be small from the narrow width of the line, at the scale of the inclination angle $\varphi$, which is estimated not to exceed $\sim 20$∘. Corrected for inclination of the star axis with respect to the line of sight, thermal broadening and instrumental resolution, line widths in the equatorial outflow reach $\sim 1.2$ km s*-1* FWHM, consistent with a flaring of $\pm 17$∘ FWHM for $\varphi=10$∘. Both CO(1-0) and CO(2-1) emissions of the equatorial outflow display remarkably similar intensity fluctuations, at the level of $\sim\pm 36$%, dominated by an eastern enhancement in the form of an arc and a south-western depression, together having the appearance of a spiral. Evidence was presented for a radial modulation of the Doppler velocity of very low amplitude, $\sim 0.16$ km s*-1*, and a period of $\sim 3.8$ arcsec. Apart from a small anti-correlation of their amplitudes, no clear correlation has been found between intensity and velocity fluctuations of the equatorial outflow. The identification of a spiral in the intensity map (Homan et al., 2018) has been thoroughly discussed.

The bipolar outflow is smoother than the equatorial outflow, with density fluctuations at the level of $\sim\pm 26$%, but reveals both an asymmetry between blue and red sides and significant differences between CO(1-0) and CO(2-1) emissions. Evidence has been presented for the emission to reach a maximum at intermediate star latitudes, with a significant depression near the poles. The reliability of this conclusion has been discussed. The distribution of temperature in space has been evaluated from the ratio between CO(2-1) and CO(1-0) emissions. A radiative transfer calculation has been performed to evaluate the absorption, which was found at the level of 20% to 30% on average, not exceeding $\sim 40$%. Accounting for it, temperature $T^{*}$ and density $d^{*}$ averaged over stellar longitude have been evaluated in the meridian plane of the star as functions of stellar latitude $\alpha$ and distance $r$ to the star. Apart from fluctuations of the mass-loss history evaluated at the level of $\pm 10$%, the mass-loss rate is found approximately constant at the level of ($1.6\pm 0.4$) 10*-7* M⊙ yr*-1* over the past 2500 years. Such fluctuations are described in terms of density in the present model because the expansion velocity has been assumed to be free of them but they could as well be blamed on velocity fluctuations with a smooth density and are most likely the effect of both density and velocity fluctuations of the flux of matter in the recent history of the star. The temperature decreases exponentially with $r$, approximately as $109\mbox{[K]}\exp(-r\mbox{[arcsec]}/3.1)$ or $106\mbox{[K]}/r\mbox{[arcsec]}$ and is maximal at intermediate latitudes. Apart from a $\sim 30$% excess close to the star ($r<\sim 2$ arcsec), the density decreases slowly with distance as approximately expected from UV dissociation. Its equatorial value depends trivially on the inclination angle of the star axis, small inclinations implying large equatorial wind velocities. Uncertainties on the measurements of temperature and density are estimated at the levels of 25% and 15% respectively.

7.2 State of the art

Our current understanding of the mass-loss mechanism of AGB stars is still relatively vague (for a recent review see Höfner & Olofsson, 2018). The basic picture is of isotropic, uniform and constant slow winds, with mass-loss rates in the range of 10${}^{-7\,\rm{to}\,-5}$ M⊙ yr*-1*, triggered by pulsation-induced shock waves, creating favourable conditions for the formation of dust grains that absorb the stellar light and transfer the resulting momentum to the surrounding gas by colliding with it. Acceleration occurs at distances from the star not exceeding a few tens of stellar radii, resulting in wind velocities at the scale of a few 10 km s*-1*. However, in the case of oxygen-rich stars like EP Aqr, serious doubts have been voiced on the ability of such a mechanism to produce sufficient acceleration (Woitke, 2006a), the dust grains being formed lacking the absorption power necessary to produce it. It has instead been proposed that micron-size iron-free and transparent silicate grains close to the star are accelerated by photon scattering (Höfner, 2008) and confirmation of this hypothesis has been subsequently obtained (Norris et al., 2012). Recently, McDonald & Zijlstra (2016) and McDonald et al. (2018) have proposed convincing arguments in favour of a scenario according to which the mass-loss episode would be triggered predominantly by pulsations rather than by radiation pressure on dust, in excellent agreement with the model predictions of Winters et al. (2000) (their B-regime that applies to short period, low luminosity stars); such a scenario would particularly apply to EP Aqr, which, with its low pulsation period of 55 days, is, according to McDonald et al., just starting to lose mass.

In any case, while the standard picture provides a satisfactory zero order approximation to the mass-loss mechanism, the physics behind deviations from it, whether bipolar outflow, equatorial enhancement or inhomogeneity in the form of arcs and lumps, is poorly understood. If we exclude spectacular events such as the formation of a super-wind at the end of the AGB phase (Schröder et al., 1999) or of thin spherical detached shells around some carbon-rich stars (Olofsson et al., 2000; Schöier et al., 2005; Maercker et al., 2010, 2016) which are apparently unrelated to the case of EP Aqr, two main candidates for breaking the symmetry of the zero order picture have been considered and studied in some detail: the presence of large convection zones in the star atmosphere and the attraction of a stellar companion or massive planet.

Major convection cells in the star cause the mass loss to occur in patches in its upper atmosphere, producing an initially clumpy wind (Freytag et al., 2017). Recent observations of the stellar atmospheres of AGB stars such as W Hya (Vlemmings et al., 2017) or R Sculptoris (Wittkovski et al., 2017) have triggered a surge of interest in this topic. To which extent such lumpiness is preserved over large time scales is unclear. Similarly, while the effect of the short period pulsations is expected to be washed away when looking at the morphology of the circumstellar envelope at large distances from the star, the role played by episodic violent events such as He flashes is also unclear (Olofsson et al., 1990). An extensively studied example of a circumstellar envelope displaying the imprint of its mass-loss history is Mira Ceti (Ramstedt et al., 2014; Nhung et al., 2016), with fragmented distinct lumps moving out radially from the star.

The co-existence in the circumstellar envelope of radial winds of different velocities can generate important inhomogeneity from their interaction. Such interacting winds have been studied extensively, in particular at the end of the AGB, before entering the Planetary Nebula phase (Schröder et al., 1999), but also in less extreme configurations (Villaver et al., 2002). In the superwind case, the strong temperature dependence of the micro-physics and dust-formation chemistry dominates the effect of gravity, suggesting that qualitatively similar but quantitatively more modest effects might be present in younger AGB stars. Of possible relevance to EP Aqr is the interaction between the AGB slow wind and a polar jet emitted by a companion accreting matter from it (Soker & Rappaport, 2000; Livio & Soker, 2001; García-Arredondo & Frank, 2004).

The main effect of the presence of a stellar companion or massive planet has been shown to be twofold (Kim & Taam, 2012): if the companion is massive enough for the reflex motion of the mass-losing star to be significant, the wind velocity is the sum of the orbital velocity of the star and the radial expansion velocity of the gas, producing helical spiral patterns having a pitch that decreases with stellar latitude; moreover, the companion attracts some gas toward the plane of its orbit, thereby producing an equatorial enhancement in its wake. How significant these effects may be depends on the parameters describing the companion and its movement, but the less massive the companion, the smaller they are. A signature of such effects is the presence of spiral, or spiral-like arc enhancements of the gas density (Mastrodemos & Morris, 1998, 1999; Kim & Taam, 2012; Homan et al., 2015). A few remarkable examples have been identified, all from carbon stars, leaving no doubt on their interpretation. Among these are LL Pegasi (Kim et al., 2017), R Sculptoris (Maercker et al., 2012) or CIT6 (Kim et al., 2013). Apart from the case of EP Aqr, no convincing example has yet been found among M or S oxygen rich stars, although it is tempting, in the absence of other compelling explanation of the presence of arcs of enhanced emissivity, to identify some of these with spiral arms (Ramstedt et al., 2014; Mayer et al., 2014). In the present case (Homan et al., 2018), however, the evidence for the presence of a spiral is much weaker than it is for carbon stars.

In the case of carbon stars simple hydrodynamic models (Woitke, 2006b) give evidence for instabilities resulting from the episodic formation of dust above the stellar surface sufficient to break the initial spherical symmetry and produce incomplete dust shells or dust arcs. At least qualitatively similar phenomena may be expected to occur in the envelope of M stars, the intrinsic instabilities inherent to astrophysical gases under the influence of a distant source of radiation being of a rather basic and general nature (Woitke et al., 2000).

The presence of a companion is not the only mechanism that has been invoked to explain the formation of an equatorial enhancement. The equatorial outflow around EP Aqr, keeping the imprint of the recent mass-loss history, might also be caused by the rotation of the star atmosphere, as described by Dorfi & Höfner (1996). The observation by Homan et al. (2018) of the emission of SO2 molecules gives indeed evidence for significant rotation at distances from the star well below 1 arcsec or so. Finally, we must mention the possible influence of magnetic fields (Vlemmings, 2013).

7.3 Discussion

When discussing the bipolar outflow detected from EP Aqr, we must keep in mind that the departures from the zero order picture that we are talking about are relatively modest. In particular, the equatorial wind velocity, $\sim 2$ km s*-1* is well below typical terminal wind velocities and below escape velocity over a significant fraction of the explored range. While of no serious concern for the results of the de-projection and the conclusions summarized in Section 7.1 above, this is of high relevance to the understanding of the mass loss dynamics. Not many oxygen-rich AGB stars have been studied to a same degree of detail as EP Aqr but for at least one of these, RS Cnc, a similar low equatorial expansion velocity has been observed (Nhung et al., 2015b). At variance with EP Aqr, RS Cnc has its axis making a large angle, $\sim 45$∘, with respect to the line of sight: the evaluation of the equatorial velocity, $\sim 2$ km s*-1*, is not significantly affected by uncertainties on its precise value. A meaningful description of such a low velocity wind implies a deeper understanding of the mechanism of acceleration, and particularly of the radial range over which it is efficient, than currently available. In particular, it requires including information from observations using other tracers, such as SiO and SO2, that explore shorter distances from the star. Further comments are therefore postponed to a forthcoming publication (Tuan-Anh et al. , in preparation).

An isotropic mass-loss shaped by a bipolar wind would produce a density that is maximal at the equator, which is far from being the case. The present results imply therefore that, even before acceleration, the flux of matter is anisotropic and peaks at intermediate latitudes. Indeed, as elaborated in Nhung et al. (2018a), whatever form is adopted for the wind velocity, a spherical flux of matter is incompatible with the present observations. This implies in turn complex mechanisms for both mass-loss and acceleration, none of the simple scenarios alluded to in the preceding paragraph producing such morpho-kinematics on its own. In particular, the pulsations that trigger the wind, corresponding to the radial first overtone (Wood, 2015, and references therein) are expected to be mostly isotropic although some non-radial $l=1$ contributions are known to be often present.

We remark that qualitatively the latitudinal dependence of the wind velocity, temperature and density resembles that observed in the Red Rectangle (Anh et al., 2015; Bujarrabal et al., 2016): a biconical enhancement of the density and temperature, peaking at intermediate latitudes, an equatorial expansion velocity of 1.6 km s*-1* together with a rotation velocity at the scale of 1 km s*-1* and a bipolar outflow velocity reaching $\sim$7.6 km s*-1*. An important difference is the temperature, typically 5 times higher in the Red Rectangle than in EP Aqr. Yet, the Red Rectangle is a proto-planetary nebula while EP Aqr is a rather young AGB star. The Red Rectangle is believed to have evolved from oxygen-rich to carbon-rich (Waters et al., 1998) and is known to be a binary, with a probable White Dwarf companion. Binarity is generally held responsible for shaping the envelope of the Red Rectangle, the companion accreting from the main star in an equatorial disc and ejecting along the line of poles a pair of high velocity (a few 100 km s*-1*), high temperature (at the scale of $\sim$1000 K) jets. However, no agreement has yet been reached on the precise mechanism (Martínez González et al., 2015); in particular, the companion polar jets may be narrow in precession about the binary axis with a period of 17.6 years (Velázquez et al., 2011) or broad along the binary axis (Thomas et al., 2013).

The similarities between the two stars suggest that binarity may be the main shaping mechanism of the EP Aqr envelope. If one takes seriously the evidence for a spiral enhancement claimed by Homan et al. (2018) and thoroughly discussed in Section 4.2, the separation of 3 arcsec between arms corresponds, for an expansion velocity of 2 km s*-1* to 850 years, at strong variance with the Red Rectangle where the orbital period is at the scale of a single year (Velázquez et al., 2011). For a total mass (main star and companion) of 2 M⊙, this means a distance between them of $(850^{2}\times 2)^{1/3}=112$ au or $\sim 1$ arcsec. The Red Rectangle hosts instead a close binary with a separation at the scale of only 1 au. In summary, the similarities remarked between the morpho-kinematic properties of the CO emission of the circumstellar envelopes of EP Aqr and the Red Rectangle contrast with the major quantitative differences between their binary properties and between their mass-loss rates, two orders of magnitude higher for the Red Rectangle than for EP Aqr: if a pair of jets were shaping the polar depressions of EP Aqr, their detection would probably be difficult (Witt et al., 2008).

To the extent that the fate of the expanding gas envelope is largely determined by what happens in the small central region, $r<\sim 1$ arcsec, a detailed knowledge of its morpho-kinematics, temperature and chemistry is obviously essential to understand the morpho-kinematics at larger distances. None of the models of the mass loss mechanism at stake in oxygen-rich AGB stars can be confirmed or denied by the results of the present work but such a thorough analysis of the properties of the circumstellar envelope at distances exceeding $\sim 250$ au severely constrains their main features. To definitely choose with confidence between the currently proposed models requires at least including what can be learned of the dynamics at shorter distances from the emission of tracers such as SiO and SO2. Even so, the under-determination of the problem of de-projection of radio observations restricts the reliability of the conclusions that can be drawn and imposes particularly critical and rigorous scrutiny.

ACKNOWLEDGEMENTS

We are deeply grateful to the referee, Professor Albert Zijlstra, for very helpful comments on the manuscript. We thank Dr Iain McDonald for comments on pulsation-induced acceleration and Dr Pierre Lesaffre for his interest in our work and useful comments. This paper makes use of the following ALMA data: 2016.1.00026.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC 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 work was supported by the Programme National Physique et Chimie du Milieu Interstellaire (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. The Hanoi team acknowledges financial support from VNSC/VAST, the NAFOSTED funding agency, the World Laboratory, the Odon Vallet Foundation and the Rencontres du Viet Nam. This research is funded by Graduate University of Science and Technology under grant number GUST.STS.DT2017-VL01.

Appendix A

Figure 31 summarizes the variables used in relation to the sky and stellar systems of coordinates used in the article. The velocity vector is written V, the Doppler velocity $V_{z}$ being its component on the $z$ axis.

Tables 5 to 7 provide additional information that is being referred to in the text.