De-projection of radio observations of axi-symmetric expanding circumstellar envelopes
P. T. Nhung D. T. Hoai P. Tuan-Anh P. Darriulat T. Le Bertre J. M. Winters P. N. Diep N. T. Phuong


firstpage–lastpage
[
    [
    December 30, 2023
=====================





  The problem of de-projection of radio line observations of axi-symmetric expanding circumstellar
  envelopes is studied with the aim of easing their analysis in terms of physics models.
  The arguments developed rest on the remark that, in principle, when the wind velocity
  distribution is known, the effective emissivity can be calculated at any point in space.
  The paper provides a detailed study of how much this is true in practice. The wind velocity
  distribution assumed to be axi-symmetric and in expansion, is described by four parameters:
  the angles defining the orientation of the symmetry axis, an overall velocity scale and
  a parameter measuring the elongation (prolateness) of the distribution. Tools are developed
  that allow for measuring, or at least constraining, each of the four parameters. The use
  of effective emissivity as relevant quantity, rather than temperature and density being
  considered separately, implies important assumptions and simplifications meaning that the
  approach being considered here is only a preliminary to, and by no means a replacement
  for, a physics analysis accounting for radiative transfer and hydrodynamics arguments.
  While most considerations are developed using simulated observations as examples, two
  case studies (EP Aqr, observed with ALMA, and RS Cnc, with NOEMA) are presented that
  illustrate their usefulness in practical cases.

  Key words: stars: circumstellar matter, AGB and post-AGB; methods: data analysis











§ INTRODUCTION




Recent years have seen high quality observations of molecular line emissions from evolved
stars become available. Such observations, in particular from NOEMA (NOrthern Extended Millimeter Array)
and from the Atacama Large Millimeter/submillimeter Array (ALMA), offer a spatial and spectral resolution calling
for analysis methods making the best possible use of it. The present work addresses the
case of expanding circumstellar envelopes of evolved stars, particularly stars populating
the Asymptotic Giant Branch <cit.>.
Such envelopes have shapes that often evolve from spherical to
axi-symmetric morphology, providing the seed for possibly more irregular configurations
later observed in post-AGB stars and Planetary Nebulae. The physics governing the breaking
of spherical symmetry is currently the subject of active research. In the case of binaries,
attraction from the companion has been shown to play an important role. Yet, many unanswered
questions, such as the role played by magnetic fields, remain to be elucidated. 


The analysis of radio observations of molecular line emissions requires, as a preliminary,
a de-projection in space. This is a largely under-determined problem: only two out of three
position coordinates are measured, those in the sky plane; the position along the line of
sight is unknown. And only one out of three velocity components is measured, that along
the line of sight, from the observed Doppler shift.  In a recent paper <cit.> 
we made general considerations on the problem of de-projection, with particular emphasis on
the differentiation between expansion and rotation and on the use of Position-Velocity (PV)
diagrams. We were then addressing issues related to both the physics of proto-stars and of
evolved stars. In the present work, we concentrate instead on expanding circumstellar envelopes.
The aim is to shed new light on the problem of de-projection in as simple terms as possible
and to understand in depth the difficulties that its solution needs to face. To do so, we
deliberately ignore complications such as arising from optical thickness or from the possible
presence of rotation, and more generally from any form of departure from exact axi-symmetry
and exact radial expansion. We mostly exploit the simplification offered by the constraint of
axi-symmetry (one relation) and from the hypothesis of radial expansion (two relations),
helping with the solution of the problem of de-projection.


The considerations that follow have no ambition at replacing the physics analysis required
by an in-depth understanding of the physics mechanisms at play. They are simply meant as a
useful preliminary step, providing helpful tools and possibly inspiring considerations on the issue of de-projection.  


§ THE FRAMEWORK



We use coordinates (x,y,z) attached to the sky plane and (x',y',z')  attached to the star (Figure <ref>).
The z axis is parallel to the line of sight, pointing away from Earth, while x is pointing
east and y north. The z' axis is the symmetry axis of the star morphology and kinematics,
making an angle φ with the line of sight. Its projection on the sky plane makes an angle
θ (position angle) with the y axis, where θ is the angle between the x axis and the x'
axis, the latter taken to be in the (x,y) plane. The transformation relations between the two systems
of coordinates read

    x'   = xcosθ+ysinθ
    
        y'   =(-xsinθ + ycosθ) cosφ + zsinφ
    
        z'   = -(- xsinθ + ycosθ ) sinφ + zcosφ
        
        x   =x'cosθ-(y'cosφ-z'sinφ)sinθ
    
        y   =(y'cosφ - z' sinφ) cosθ + x' sinθ
    
        z   =y'sinφ + z'cosφ










However, in much of what follows, we redefine the y axis as the projection of the star axis
on the sky plane, which is equivalent to setting θ=0. In this case x=x' and the
transformation relations between (y,z) and (y',z') read:

    y'=ycosφ + zsinφ        		y=y'cosφ - z' sinφ
    
    		 z'=-ysinφ + zcosφ       	z=y' sinφ + z' cosφ

In order to illustrate our arguments as simply as possible we use a model (in the remainder
of the article we refer to it as “the simple model”) in
which the wind velocity is radial and independent of the distance from the star with a
dependence on star latitude α of the form

    V=V_polesin^2α + V_eqcos^2α=V_0( 1 - λcos2α)

with the polar velocity V_pole=V_0(1+λ) and the equatorial velocity V_eq=V_0(1 - λ)
(Figure <ref> left). For λ between 0 and 1 one obtains prolate velocity distributions (bipolar
outflow) and for λ between -1 and 0 oblate velocity distributions (equatorial outflow),
λ=0 corresponding to an isotropic (spherical) velocity distribution. Here V_0 defines the
velocity scale but its precise value is irrelevant; in practice, we set it at 5 km s^-1 for the
purpose of illustration. The Doppler velocity at space point (x,y,z) is

    V_z=V sinζ

where ζ   is the angle between the (x,y,z) direction and the plane of the sky.
We also define R=√(x^2+y^2)=z/tanζ and r=√(R^2+z^2)=z/sinζ=R/cosζ.

From Relations <ref>, <ref> and <ref>, using  sinα=z'/r and replacing
z' by its expression in Relation <ref> we can express the Doppler velocity as a function
of ζ, φ and ψ, defined as the position angle measured with respect to the projection
of the star axis on the sky plane, y=Rcosψ:

    V_z    = V_0 (z/r)[1 - λ (1 - 2sin^2α)]
       =V_0(z/r)[1 - λ + 2 λ(z cosφ  - y sinφ)^2/r^2]


    V_z  =    V_0 sinζ [ 1 - λ + 2 λ ( sinζcosφ - cosζcosψsinφ)^2]
    
      =    V_0 [( 1 - λ ) sinζ + 2 λsin^3ζcos^2φ
       - 4 λsin^2ζcosζsinφcosφcosψ
        + 2λsinζcos^2ζsin^2φcos^2ψ] 5a


Taking the derivative along the line of sight (R and ψ being fixed) and using the identities
d/dz=[d/d(sinζ)][d(sinζ)/dz]=r^-1 (cos^2 ζ) d/d(sinζ) and d(cosζ)/d(sinζ)= - tanζ we obtain:

    rdV_z/dz=V_0[(1 - λ)cos^2ζ + 2λ P]  5b

with

    P=   3cos^2φsin^2ζcos^2ζ
       -cosψsin2φcosζsinζ(2cos^2ζ-sin^2ζ) 
       + cos^2ψsin^2φcos^2ζ(cos^2ζ-2sin^2ζ) 5c


Radio astronomy measurements are in the form of a data-cube with elements f(x,y,V_z) measuring the brightness
along the line (x,y) normal to the sky plane at Doppler velocity V_z. It is convenient to define
the effective emissivity at space point (x,y,z) as 

    ρ(x,y,z)=f(x,y,V_z)dV_z/dz

and the measured intensity as 

    F(x,y)= ∫f(x,y,V_z)dV_z=∫ρ(x,y,z)dz




As mentioned in the introduction, the assumption of pure radial velocity and of axi-symmetry should help
with the de-projection of the effective emissivity using Relation <ref>. How much this is true in
practice is the subject of the present article.








A number of symmetries are apparent in Relations <ref>. Changing φ in -φ  and ψ in
180^∘ - ψ leaves V_z invariant, so does also changing ψ in -ψ; changing
φ in 180^∘ - φ and ζ in -ζ changes V_z in - V_z. Accordingly, it is
sufficient to limit the ranges of φ and ψ to respectively [0^∘, 90^∘] and [0^∘, 180^∘]. 



A remarkable consequence of the above relations is the ability to estimate the r-dependence of the
effective emissivity independently from the form chosen for the wind velocity (V_0,λ,φ)
as long as V does not depend on r (Figure <ref> right). Indeed, consider a line passing by the star
and having direction (ψ,ζ), namely making an angle ζ with the plane of the sky and
projecting on it at position angle ψ; this line corresponds to a single value of V_z, independent
of both R and r as long as the radial expansion velocity is constant on it, namely independent of r.
Therefore, as long as both the real wind velocity, V(ψ,ζ), and the wind velocity used for
de-projection, V'(ψ,ζ) are both independent of r (but depend on ψ and ζ, generally
in different ways)  the de-projected line is also a straight line, having the same projection on the
plane of the sky; its direction (ψ',ζ') is related to the direction (ψ,ζ)  by the
relations ψ'=ψ and V'(ψ,ζ')sinζ'=V(ψ,ζ)sinζ.  De-projection simply
transforms ζ into ζ' and the ratio z'/z is independent of R. As a result, the true
effective emissivity ρ and the de-projected effective emissivity ρ' have the same dependence
on R on each of these two lines, (ψ,ζ) for the former and (ψ,ζ') for the latter:
they are proportional to a same data-cube element f(x,y,V_z), independently from R . To the extent
that the true wind velocity and the wind velocity used for de-projection are not too different, the
r-dependence of the de-projected effective emissivity is similar to the r-dependence of the true
effective emissivity, both being dominated by the R dependence of the Doppler velocity spectrum.


§ GENERAL CONSIDERATIONS

De-projection implies using Relation <ref> to calculate the effective emissivity ρ(x,y,z) from the measured brightness f(x,y,V_z) by associating to each value V_z of the measured Doppler velocity spectrum a point (x,y,z) in space. This is only possible if the relation giving V_z as a function of z can be inverted into a relation giving z as a function of V_z. In general, an extremum of the dependence of V_z on z will generate in its vicinity two values of z for a same value of V_z. In principle, this should prevent de-projection as one does not know how to share the brightness measured at Doppler velocity V_z between the two corresponding space points; this issue will be discussed in some detail in Section <ref>. As the relation between z and sinζ is one-to-one, z=Rtanζ=Rsinζ / √(1-sin^2ζ), the extrema of V_z vs z are the same as of V_z vs sinζ.






In order to obtain some insight into this question, we consider a sample of uniformly distributed (λ,φ,ψ) triplets. Figure <ref> (left) displays the values of sinζ at which an extremum of V_z vs z is found; they concentrate near sinζ=±1. To see it more directly, we display in Figure <ref> a set of representative functions sinζ vs V_z. For sinζ=±1,  V_z=± V_0(1+λcos2φ) and when sinζ departs from ±1, V_z increases or decreases depending on the sign of λsin2φcosψ. Globally, V_z must increase with z; if it starts in the wrong direction, it needs to turn back and one obtains an extremum. When such extrema are in the vicinity of sinζ=±1, they are not too harmful for de-projection: they simply mean that each end region of the Doppler velocity spectrum has to be assigned a broad range of sinζ values in the vicinity of the line of sight. However, the farther away they are from sinζ=± 1, the larger the fraction of the Doppler velocity spectrum that becomes unsuitable for de-projection. In a majority of cases (68%) there is no “harmful” extremum, defined as having |sinζ|>0.9.


The right panel of Figure <ref> displays the region of the (λ,φ) plane
where no “harmful” extremum occurs, namely where dV_z/dz does not cancel in the
|sinζ|<0.9 region. In this region, as clearly illustrated in Figure <ref>,
large values of |V_z| are associated with large values of |sinζ| and
de-projection can be unambiguously performed as long as V_0 is large enough.
In the complementary ambiguous region, on the contrary, the larger values of |V_z|
are associated with intermediate values of |sinζ|, the larger values of
|sinζ| being now associated with lower values of |V_z|: the Doppler
velocity spectrum is folded on itself, the larger values of |V_z| being
associated with two values of |sinζ|, making de-projection ambiguous
and unreliable. The strong qualitative difference between the two regions will
be seen to play a major role in the arguments developed in the present study. 




§ POSITION ANGLE OF THE PROJECTION OF THE STAR AXIS ON THE SKY PLANE

The invariance of Relations <ref> when ψ changes sign implies an exact symmetry of the data-cube with respect to the plane of position angle θ, perpendicular to the sky plane and containing the star axis and the line of sight. It is indeed a consequence of axi-symmetry and is expected to apply for any axi-symmetric model, not just the simple model used here for illustration. In principle, the symmetry plane, and therefore the value of θ, can be simply found by minimizing the quantity

    χ^2_rot=∑ [(f(x,y,V_z) - f(x^*,y^*,V_z)]^2/(Δ f)^2
 
where (x^*,y^*) is the symmetric of (x,y) with respect to direction θ:	

    x^*=-(x cos2θ + y sin2θ )     		y^*=-(x sin2θ - y cos2θ)

and where Δ f is the uncertainty on the f measurement. In the present section we use a simulation having a wind of the form given in Relation <ref> and we simply take as Δ f the quadratic sum of the rms deviations of f from its mean in the vicinity of each of (x,y,V_z) and (x^*,y^*,V_z).


As an illustration of the procedure, Figure <ref> (left) shows the dependence of χ ^2_rot on θ-θ_0 for λ=0.5 and various values of φ, with θ_0 being the value of θ used to produce the simulated effective emissivity. The angle θ is undefined in two obvious cases: for an isotropic velocity distribution, λ=0, and for a star axis parallel to the line of sight, φ=0; the latter case is a trivial effect of geometry, the elementary solid angle being dΩ=sinφ dθ dφ rather than simply dθ dφ. Indeed, the dependence of χ^2_rot on θ - θ_0 is observed to display a steep minimum at 0 as long as λ and φ are not too close from zero.


 

In order to assess the accuracy of the θ measurement, we calculate the relative increase
S_rot of χ^2_rot associated with a shift of ± 1^∘ from θ=θ_0:
S_rot=1/2[χ^2_rot(θ=θ_0 - 1^∘)+χ^2_rot(θ=θ_0 + 1^∘)]/χ^2_rot(θ=θ_0).
The larger S_rot and the farther away from unity, the better defined is the value of
θ. Figure <ref> (right) displays the dependence of S_rot on λ
and φ. Its actual value when applied to real observations depends on the relevant
uncertainties, usually caused by the lumpiness of the effective emissivity rather than by
noise; but its behaviour in the (λ,φ) plane remains essentially the same as
found here. However, a systematic rather than random deviation from axi-symmetry may affect
the measurement of θ, as was already discussed in the case of rotation by <cit.>.
However, a systematic violation of axi-symmetry by the effective emissivity or by the wind
radial expansion velocity of the form 1 + ϵcos(ω-ω_0) does not shift the
value of θ but simply broadens the minimum of χ^2_rot, because, contrary to rotation,
it distorts in a same way the red-shifted and blue-shifted hemispheres.



§ MAGNITUDE OF THE WIND VELOCITIES: WIDTH AND OFFSET OF THE DOPPLER VELOCITY SPECTRUM

In general, the radial wind velocity V is confined between two finite values; in the case
of the simple model, these are V_eq and V_pole, respectively V_0(1-λ) and
V_0(1+λ). However, at a given point (x,y) in the sky plane, the measured Doppler
velocity varies between two values V_zmin and V_zmax that are not simply related to the
above. Yet, when choosing a velocity distribution with which to de-project the effective
emissivity, it is essential to have some idea of its scale, meaning the value of V_0 in
the case of the simple model. In principle, having chosen a pair (λ,φ) for
de-projection, the values of V_zmin/V_0 and V_zmax/V_0 that they generate in each
pixel are known. One should then choose for V_0 the ratio between the measured values
of V_zmin and V_zmax and the model values of V_zmin/V_0 and V_zmax/V_0. The
values of V_zmin/V_0 and V_zmax/V_0 depend on the pixel.
In general, different pixels produce different values of V_0 and their mean, or better their maximum, should be
retained for de-projection. Note that changing z in -z and cosψ in -cosψ
changes V_z in -V_z but leaves dV_z/dz invariant; therefore it changes V_zmin
in -V_zmax and V_zmax in -V_zmin: it is sufficient to confine ψ to the
[0^∘,  90^∘] interval. Figure <ref> displays the maps of V_zmin/V_0 and V_zmax/V_0
in the (λ,φ) plane for three intervals of ψ, each 30^∘ wide,
covering between 0^∘ and 90^∘. Qualitatively, the dependence on ψ is weak.
The main features reflect the effect of the width of the Doppler velocity spectrum,
which is large on the descending and small on the ascending diagonal. On the former,
from  (λ,φ)= (+1,0^∘) to (-1,90^∘), namely from a bipolar outflow
parallel to the line of sight to an equatorial outflow having its axis in the sky plane,
V_zmin/V_0 varies between -1 and -2 and V_zmax/V_0 between 1 and 2. On the latter,
from (λ,φ)=(-1,0^∘) to (+1,90^∘), namely from a bipolar to an equatorial
outflow, both in the sky plane, V_zmin/V_0 varies between -1 and 0 and V_zmax/V_0 between 0 and 1.




When de-projecting the effective emissivity using a (λ,φ) pair of parameters,
agreement between the values obtained in each pixel for V_0 as a function of ψ and R is a
useful indicator of the suitability of the particular (λ,φ) pair to describe
the observations. The ratio Q=V_rms/V_mean between their rms value V_rms and their mean
value V_mean, calculated over the whole image, when too large, can be used to reject
unsuitable wind configurations: the study of the
dependence on position angle of the width and offset of the Doppler velocity spectrum is
not only a tool to obtain an evaluation of the scale of the space velocity (here V_0)
but also, in principle, to reject unsuitable values of the (λ,φ) pair. This
is illustrated in Figure <ref>, which displays the dependence of Q on (λ,φ) for
a few typical simulated wind configurations. In all cases a steep minimum of Q is obtained
in a narrow region of the (λ,φ) plane containing the values used in the simulation,
suggesting that an important fraction of the (λ,φ) plane could be eliminated
by simply requiring Q not to exceed some threshold. In practice, however, as will be seen
in Section <ref>, the minimum of Q is much less steep for real than for simulated data.

We remark that if V_0 is slightly overestimated, its largest values will de-project in
a region of the observed Doppler velocity spectrum where there are no data and will accordingly
set the de-projected effective emissivity to zero. On the contrary, if V_0 is slightly
underestimated, the larger values of the observed Doppler velocity spectrum will be ignored.
It is therefore important, in practice, to make sure that the obtained value of V_0 is optimal
and, if necessary, to fine-tune it.







§ MEASURING THE PROLATENESS PARAMETER Λ AND THE INCLINATION Φ OF THE STAR AXIS WITH RESPECT TO THE LINE OF SIGHT

The results obtained in the preceding two sections are largely independent from the particular
form of the dependence of the effective emissivity on stellar latitude. They provide reliable
estimates of the scale V_0 of the space velocity distribution and of the position angle θ
of the projection of the star axis on the sky plane. Moreover, they eliminate regions of the
(λ,φ) plane that are unsuitable for de-projection. We are then left with two parameters,
λ and φ, to be measured in the case of the simple model. In the general case,
several parameters will be necessary to describe the wind velocity in place of the pair (V_0,λ)
and λ must be seen as a measure of the effective prolateness of the wind velocity distribution.


In the present section, we exploit the constraint resulting from the requirement of axi-symmetry
of the effective emissivity to help with the measurement of the (λ,φ) pair. This
constraint is less useful in some cases than in others. For example, in the case of a wind velocity
having its axis parallel to the line of sight any pair (λ,φ=0) used in de-projection
will produce an axi-symmetric effective emissivity, independently from the value assumed for λ.

 

 In order to understand under which conditions the constraint of axi-symmetry is strong, we use a simple measure of the amount of axi-symmetry of the de-projected effective emissivity about axis (θ,φ): 
 
    χ^2_axi=∑[ρ(r,α,ω)  - <ρ(r,α)>]^2/(Δρ)^2
	
 where ω is the stellar longitude (x'  = rcosαcosω, y' = rcosαsinω
 and tanω = y'/x'); ρ (r,α,ω) is the effective emissivity de-projected using
 a wind configuration of axis (θ,φ), effective prolateness λ and velocity scale
 V_0 estimated from the procedure described in the preceding section;  is its
 mean value at (r,α), averaged over longitude ω; Δρ is the uncertainty
 attached to the evaluation of ρ. The sum extends over the whole space over which measurements
 are available. The evaluation of Δρ=Δ f |dV_z/dz|, where Δ f is the uncertainty on the de-projected
 data-cube element, uses here the same estimate of  Δ f as used for χ^2_rot in Section <ref>, namely the
 rms deviation of f from its mean in the vicinity of the de-projected space point. In real cases,
 however, the definition of Δ f is delicate and needs to account for both experimental uncertainties
 and uncertainties attached to the de-projection; we discuss this point in more detail in Section <ref>.
 In practice the calculation of χ^2_axi proceeds as follows: one chooses a circle having the
 star axis as axis and defined by its position (R',z') in the star frame and scans over it by
 varying the stellar longitude ω, each time calculating the space position (x,y,z) and the associated
 value of V_z for the values of (V_0,λ,φ) of the axi-symmetric wind used in de-projection.
 
 The main weakness of this method is its mishandling of cases where two different values of z/r
 are associated with a same bin of Doppler velocity. Such bins contain contributions from each of
 the two regions but the de-projection algorithm wrongly assigns their total content to each of
 the two regions, generating double-counting. As was already remarked in Section <ref>, one does not
 know how to share it between the two regions and such bins are unsuitable for de-projection.
 We have seen in Section <ref> that for z/r=±1, V_z=± V_0(1 + λcos2φ) and that
 when z/r departs from ±1, V_z increases or decreases depending on the sign of
 λsin2φcosψ. Rather than simply ignoring the ambiguous intervals of Doppler
 velocity between V_zmin and -V_0(1+λcos2φ)>V_zmin and between
 V_0(1+λcos2φ)<V_zmax and V_zmax, we simply assign to each of the two associated
 z values half the de-projected emissivity, thereby avoiding double-counting.

 Figure <ref> illustrates the dependence of χ^2_axi on parameters λ and φ
 for the same simulated wind configurations as displayed in Figure <ref>. The effective emissivity
 used to produce the simulated data-cube is of the form 1/r^2 for 0.5<r<5 arcsec with no dependence on
 stellar latitude in order not to bias the evaluation of χ^2_axi or, more precisely, not to complicate
 its interpretation. For r>5 arcsec, the effective emissivity is taken to cancel and for r<0.5 arcsec  it is
 taken to be constant.
 
 In all cases χ^2_axi is minimal at the simulated value of (λ,φ) but in some cases
 a broader region of the (λ,φ) plane is observed to be equally acceptable. The minimum
 is better behaved along the descending diagonal (wind velocities out of the sky plane) than along
 the ascending diagonal (wind velocities near the sky plane): in the former case χ^2_axi is
 more efficient than Q to constrain the (λ,φ) pair while in the latter case Q is more
 efficient than χ^2_axi.

 

 Using φ=0 for de-projection deserves a special comment; in this case, a given circle (R',z')
 is made of points having all the same value of z=z' and the same value of R=R', therefore a same
 value of V_z and a same value of dV_z/dz: χ^2_axi is simply a measure of the axi-symmetry
 of the data-cube about the z axis for Doppler velocities |V_z|<(1+λ)V_0. In particular,
 when the real wind and the wind used for de-projection are both axi-symmetric with respect to the
 z axis, with respective prolateness λ_true and λ_deproj, χ^2_axi=0 is of
 no help to constrain the value of λ_true. The same is true of the study of V_zmin,
 V_zmax and Q, in the approximation where the wind velocity and mass loss rate do not depend
 on r: whatever the value of λ the extrema of the velocity spectrum are reached in all
 pixels at a same value of V_zmax and V_zmin=-V_zmax. The result of the de-projection depends
 on the location of λ_true and λ_deproj with respect to the value -0.2 that
 separates ambiguous from unambiguous cases. When both λ_true and λ_deproj are
 larger than -0.2, the maximum of |V_z| is reached at the poles and the value of V_0 used
 for de-projection is the true value multiplied by (1+λ_true)/(1+λ_deproj).  However,
 when one of the λ parameters takes values that vary between -0.2 and -1, the maximum
 of |V_z| is reached at values of |sinζ| that vary between 1 and 1/√(3), the
 corresponding values of |V_z| varying between 0.8V_0 and 4V_0/[3√(3)]∼ 0.77V_0
 with a minimum at λ=-1/2 where |sinζ|=1/√(2) and |V_z|=1/√(2)V_0∼ 0.71V_0.
 This near independence on λ (∼± 3V_0/4) of V_zmin and V_zmax when φ=0
 and λ<-0.2 could already be seen in Figure <ref>.
 In general, when using φ=0 and λ<-0.2 for de-projection of a (λ_0,φ_0) data-cube
 located in the unambiguous region will require the use of a V_0 value significantly larger than the true value
 and will deproject the large values of |V_z|, associated with large values of |sinζ|, to intermediate
 values of |sinζ| resulting in relatively low values of χ^2_axi.


§ CASE STUDIES: RS CNC AND EP AQR
 
In the preceding sections, arguments were developed using simulated rather than real observations. In practice, the morphology of circumstellar envelopes of evolved stars are far from being as smooth and well-behaved as those simulated in the present study. The question of the practicability of using the arguments and the tools developed in the preceding sections remains therefore open at this stage. While each particular case must be considered separately, we find it useful to devote the present section to two case studies in order to get some idea of the nature and magnitude of the difficulties that one may have to face. 




 §.§ RS Cnc

The first star selected for this purpose is RS Cnc, an S-type AGB star (CSS 589 in Stephenson's
1984 catalogue) in the thermally-pulsing phase <cit.>, which can be considered as
representative of its family, the morphology of its envelope being quite clumpy but not excessively.
Analyses of Plateau de Bure observations of its CO(1-0) and CO(2-1) emissions have been published
earlier <cit.>. Here, we use CO(2-1) observations
from the upgraded Plateau de Bure interferometer (NOEMA) having a beam size (FWHM) of 0.44×0.28
arcsec^2 and a spectral resolution of ∼ 0.2 km s^-1. The noise level per data cube element
is 1.6 mJy beam^-1 for pixels of 0.07× 0.07 arcsec^2. In order to improve the signal to noise
ratio, we limit the study to the region R<4 arcsec and we group pixels by 3×3=9, meaning
0.21× 0.21 arcsec^2; we use as experimental uncertainty on the brightness the quadratic
sum of the noise and of 20% of the measured value. Figure <ref> displays sky maps of the measured
intensity multiplied by R and of the mean Doppler velocity, together with the integrated Doppler velocity spectrum.


As already apparent from the map of the mean Doppler velocity, the projection of the star axis on
the sky plane is nearly north-south oriented. Indeed we find that χ^2_rot is minimal for
θ_0 ∼ 7^∘ and we display in Figure <ref> (left) the dependence on θ-θ_0
of its normalized value (divided by the number of degrees of freedom).  The value of χ^2_rot
at minimum, ∼ 4, is the result of the relative lumpiness of the measured brightness. The uncertainty
on θ is accordingly poorly defined; as an indication of how well θ is measured we
quote as effective uncertainty the value associated with a 10%  increase of χ^2_rot with
respect to minimum, namely θ =7^∘±10^∘. Indeed, masking measured brightness values
below noise level (one sigma) minimizes χ^2_rot at θ∼ 10^∘ instead of ∼ 7^∘. 


Having obtained an evaluation of  θ we rotate the data-cube by 7^∘ about the line of sight in order
to have effectively θ=0 and proceed with the evaluation of the scale V_0 of the wind velocity and the
exclusion from the (λ,φ) plane of regions unsuited for de-projection as was done in Section <ref>
using simulated data. The dependence of Q on λ and φ is displayed in the central panel of
Figure <ref>. At first sight, the resemblance with the upper-middle panel of Figure <ref>
(λ=0.4, φ=45^∘) is striking; however, the minimum is now considerably less steep than
it was for simulated data and extends over a broad region that excludes wind velocities closely confined
to the sky plane, whether bipolar or equatorial outflows. The resulting V_0 values are displayed in the
right panel of Figure <ref> and range between ∼ 4 and ∼ 15 km s^-1, the lower values being associated with
prolate bipolar outflows near the line of sight or oblate equatorial outflows having their axis near the
sky plane, namely wind configurations in the unambiguous de-projection region. The larger values are instead
confined to the ambiguous de-projection region of the (λ,φ) plane. 


A major cause of error in the evaluation of V_0 is the difficulty to measure accurately V_zmin and/or V_zmax
in pixels where the brightness is close to noise level at the extremities of the Doppler velocity spectrum. In
calculating in each pixel the values of V_zmin and V_zmax we used brightness in excess of two noise σ's
and checked that the result was essentially unaffected when using instead one or three noise σ's. The left
panels of Figure <ref> display distributions of the values obtained for V_0 in each pixel (two values per pixel)
for some representative values of the (λ,φ) pair. Ideally, if the selected wind configuration
describes the data well, the V_0 distribution must have the shape of a narrow peak centered at the proper V_0
value. Such a peak is indeed visible for some values of the (λ,φ) pairs, but it is broad and
accompanied by a low V_0 tail: exploiting the information contained in the dependence over the sky plane
of V_zmin and V_zmax is clearly more difficult when dealing with real data than it is with simulated
data. In pixels where the brightness is close to noise level at the extremities of the Doppler velocity spectrum
|V_zmin| and/or |V_zmax| are under-evaluated and the same is therefore true for V_0. Pixels at larger
values of R, being associated with lower intensities, are more likely to be of that kind. Indeed, as can
be seen in Figure <ref>, the distribution of V_0 obtained inside the circle R<3 arcsec is narrower
and less contaminated by a low V_0 tail than inside the circle R<4 arcsec. The distributions displayed
in the left panels of Figure <ref> suggest that a figure of merit revealing the presence of a peak in
the V_0 distribution might be more efficient than Q in exploiting the information carried by V_zmin and V_zmax.
We use a simple algorithm to evaluate, for each value of the (λ,φ) pair the peak to tail ratio
of the V_0 distribution, P/T. The distribution of P/T over the (λ,φ) plane is displayed
in Figure <ref> (right). It is consistent with the information carried by Q but is much more selective:
a reasonable cut Q>0.36 is only a factor 1.06 above minimum while an equivalent cut P/T<2 is a factor ∼ 4
above minimum and a factor 2.5 below maximum, allowing for safely rejecting regions of the (λ,φ) plane having P/T<2.   




In order to further restrict the acceptable region of the (λ,φ) plane, we still need to exploit
the constraint of axi-symmetry as was done in Section <ref> using simulated observations. This, however, cannot
be done reliably in the regions of ambiguous de-projection (right panel of Figure <ref>). As these have a large
overlap with the regions disfavoured by the analysis of the (λ,φ)-dependence of V_zmin and
V_zmax, we may exclude them from the analysis by rejecting regions containing a “harmful” extremum of
the V_z vs z relation as defined in the right panel of Figure <ref>, namely regions where dV_z/d_z
cancels in the |sinζ| < 0.9 interval. When an extremum occurs at |sinζ| > 0.9, we share
the de-projected emissivity equally between the two associated values of sinζ. The distribution of
χ^2_axi in the (λ,φ) plane is displayed in Figure <ref> together with the boundaries
associated with regions containing a “harmful” extremum and with the contours associated with Q=0.36 and P/T=2.
The main contribution of χ^2_axi is to disfavour wind configurations close to spherical with an axis
at intermediate inclination with respect to the sky plane, leaving a relatively narrow region of acceptable
bipolar outflows in the upper-left quadrant of the (λ,φ) plane with λ between ∼ 0.3
and ∼ 0.8 and φ between ∼ 15^∘ and ∼ 45^∘.  



The joint analysis of CO(1-0) and CO(2-1) emissions presented in <cit.> gave
(V_0, λ, θ, φ)=(5.0±0.2 , 0.50±0.02, 9^∘± 6^∘, 38^∘± 2^∘) ,
however allowing for some r-dependence of the expansion velocity. This corresponds indeed to values of
Q<0.36, P/T>2 and χ^2_axi<10 favoured by the present study. We display in Figure <ref>
(right) the dependence of the de-projected emissivity, multiplied by r^2 on respectively r, ω
and z'/r=sinα for three different values of the (λ,φ) pair in the favoured region,
(λ,φ)= (0.5, 40^∘), (0.6, 30^∘) and (0.7, 20^∘) respectively. All three distributions
are averaged in the sphere r<4 arcsec. A remarkable result is the independence of the r-dependence of
the de-projected effective emissivity on the value of (V_0, λ, φ) used for de-projection,
a result that had been anticipated and explained in Section <ref> using Relations <ref>. However, the longitudinal
dependence differs significantly from uniform, with an excess in the 180^∘<ω<360^∘ hemisphere
compared with 0<ω<180^∘. The latitudinal dependence displays a clear asymmetry with respect to
the star equator.  Table <ref> summarizes the results.

Commenting further on these results goes beyond the scope of the present study. Physics arguments need now
to be used for interpreting the observed behaviour of the de-projected emissivity and its relation with the
probably much too simple form assumed for the wind configuration. Yet, the present results are of considerable
help in constructing a physics model of the morphology and kinematics of the circumstellar envelope and have
provided a deep insight in the constraints that such a model has to obey.







 §.§ EP Aqr

The second star used for illustration is EP Aqr, an oxygen-rich M type AGB star that is probably at the beginning of its evolution on the thermally pulsing phase <cit.> in spite of observations of trailing gas <cit.> suggesting a mass loss episode at the scale of 10^4 to 10^5 years. Recently, observations of ^12CO(1-0) and ^12CO(2-1) emissions using the IRAM 30-m telescope and the Plateau de Bure Interferometer have been reported <cit.>. Here, we use CO(1-0) and CO(2-1) observations made in Cycle 4 of ALMA operation <cit.>. The beam size (FWHM) is respectively 0.78×0.70   and 0.33×0.30 arcsec^2 and the noise respectively 8 and 7 mJy beam^-1 for pixels of respectively 0.2× 0.2 and 0.1×0.1 arcsec^2 and Doppler velocity bins of 0.2 . Similar but different ALMA observations, including also SiO and SO_2 emission, have recently been presented by <cit.>.




We restrict the analysis to the sky plane region having R<8 arcsec. Figure <ref> displays sky maps of the deviation from unity of the ratio Δ(x,y)=F(x,y)/<F(x,y)> between the measured intensity F(x,y) and its average <F(x,y)> over position angle. It gives evidence for approximate isotropy in the sky plane, suggesting that the wind is either spherical or axi-symmetric with the symmetry axis near the line of sight. Indeed, χ^2_rot is found to display a very broad minimum, centered at θ_0 ∼  -20^∘ for CO(1-0) emission and at θ_0 ∼ -31^∘ for CO(2-1) emission; its dependence on θ - θ_0 is shown in Figure <ref>. Using the same criterion as for RS Cnc to measure the uncertainty on its measurement, namely a 10% increase of χ^2_rot, gives θ= -20^∘± 37^∘ and -31^∘± 43^∘ respectively, consistent with a common value of -25^∘± 28^∘. The agreement between the CO(1-0) and CO(2-1) results suggests that the deviation from sphericity is real and the broad distribution of χ^2_rot as a function of θ-θ_0 indicates that the star axis must make a small angle with the line of sight. Also shown in Figure <ref> (middle panels) is the dependence of Q on λ and φ. It remains nearly constant over a very large region of the (λ,φ) plane. The reason is that V_zmin and V_zmax are nearly invariant over the sky plane. In order to quantify better this statement, we integrate the Doppler velocity spectrum between R=1 arcsec and R=8 arcsec and display it in 18 bins of position angle ψ, each 20^∘ wide. The high sensitivity that results allows for measuring in each bin V_zmin and V_zmax with a precision of ∼ 0.1 . We obtain this way measurements of the width W_vz and the offset Δ_vz of the Doppler velocity spectrum in each bin of ψ. Their mean and rms values for CO(1-0) and CO(2-1) emission respectively, relative to the spectrum width, are
<W_vz>=21.30  and 21.61 , Rms(W_vz)/<W_vz>=2.1 10^-2 and 2.9 10^-2 and Rms(Δ_vz)/<W_vz>=1.2 10^-2 and 2.0 10^-2 respectively. The latter two quantities are compared with the predictions of the simple model in the right panels of Figure <ref>. Compared with Q they are more reliable discriminants in rejecting regions of the (λ,φ) plane but leave a very broad acceptable region along the descending diagonal. The situation here is very different from what was found for RS Cnc: it would not have helped to use a discriminant such as P/T because the V_0 distribution is extremely narrow in the low Q region of the (λ,φ) plane.

Figure <ref> displays the dependence of χ^2_axi on λ and φ. We use as uncertainty on the measured brightness the quadratic sum of 20% of its value and of the noise. For both CO(1-0) and CO(2-1) emissions, large values of φ are excluded; taking into account the constraint imposed by the study of V_zmin and V_zmax leaves a large region in the upper left quadrant of the map acceptable for de-projection.  The result obtained earlier by <cit.>, (V_0, λ, θ, φ)=(6.0 , 0.67, -36^∘, 13^∘) is contained in this region. We select three representative (λ,φ) pairs for purpose of illustration: (0, 0^∘), (0.35, 6^∘) and (0.7, 12^∘). The corresponding parameters are listed in Table 2 and the dependence on r, ω and sinα of the associated de-projected effective emissivity multiplied by r^2 inside the sphere r<8 arcsec is displayed in Figure <ref>. The difference between the CO(1-0) and CO(2-1) radial distributions is understood as being caused by the different temperature dependence <cit.>. The smaller λ, the larger V_0 and the stronger the concentration of the effective emissivity near the equatorial plane. Indeed one evolves from a bipolar outflow for λ = 0.7 with no significant equatorial enhancement of effective emissivity to a spherical wind producing a strong enhancement of effective emissivity at equator for λ = 0. Only physics arguments can help choosing between the possible wind configurations. In particular, the presence of a companion star or massive planet in the equatorial plane would probably favour an enhancement of the effective emissivity while not making the wind velocity much deviate from spherical. A bipolar wind associated with an isotropic effective emissivity would be more difficult to explain. However, such arguments are well beyond the scope of the present article.






§ SUMMARY AND CONCLUSION

The above analysis of the problem of de-projection of radially expanding axi-symmetric circumstellar envelopes has provided a deep insight into its main features and has devised tools that help with its solution. However, it does not offer a substitute to a detailed analysis that takes into account the physics of the mechanisms at play, such as hydrodynamical constraints and radiative transfer considerations. It is only meant to shed light on some of its intricacies and to set the frame for preliminary considerations that can help the construction of a realistic physical model. The effective emissivity used throughout the paper to describe observations is a convenient quantity that offers simplicity but hides the difficulty of disentangling the effects of temperature from those of density. Its assumed axi-symmetry is often a good approximation but is only valid for envelopes that are sufficiently optically thin and the assumption of radial expansion, implying the absence of rotation, is often violated in the later stages of the star evolution on their way to planetary nebulae. It is therefore essential to keep in mind the limited scope of the results of the present work.






With these caveats in mind, we summarize below the main results that have been obtained:


i) The position angle θ of the projection of the star axis on the sky plane has been shown to minimize a quantity, χ^2_rot, independently from the particular form taken by the effective emissivity; its minimization, as long as it is made over the whole data-cube and not simply on its projection on the sky plane, makes optimal use of the available information. Only in cases where the problem has no solution, either because the star axis is close to the line of sight or because the wind velocity and the effective emissivity are nearly isotropic, is θ ill-defined.   


ii) A good evaluation of the scale of the wind velocity -V_0 in the simple model-has been obtained from a study, in each pixel, of the dependence on position angle ψ of the end points V_zmin and V_zmax of the Doppler velocity spectra. To this effect, the rms deviation of the velocity scale relative to its mean over the sky plane, Q, was found to be sufficient, in the case of simulated data, to quantify the agreement with a model. Its minimization helps in eliminating regions of the prolateness (λ) versus inclination (φ) plane which are unsuitable for de-projection. The method is particularly efficient when the wind velocity is confined near the sky plane, either as strongly prolate with axis close the sky plane or as strongly oblate with axis close to the line of sight. However, in the case of real data, a more careful analysis is needed in order to exploit the richness of the information contained in the dependence of V_zmin and V_zmax on λ and φ. A careful evaluation of V_zmin and V_zmax, taking the noise level in proper account, is mandatory and algorithms allowing for separating the V_0 peak from a background of improper values may be helpful.  

iii) Having obtained sensible estimates of θ and V_0, one can further constrain the (λ,φ) pair by imposing axi-symmetry on the de-projected effective emissivity. To this effect, a quantity χ^2_axi has been constructed, which is minimal for maximal axi-symmetry.
Contrary to the minimization of Q, the minimization of χ^2_axi is inefficient when the wind velocity is confined near the sky plane. In such cases, the evaluation of the optimal value of the (λ,φ) pair relies more on the constraints imposed by the distribution on the sky plane of the width and offset of the Doppler velocity spectra than on the constraints imposed by the requirement of axi-symmetry.

iv) Having obtained estimates of the orientation of the star axis, θ and φ, of the scale of the wind velocity, V_0, and of the effective prolateness of its distribution, λ, we are then in a position to de-project the effective emissivity as long as a single value of z/r is associated with each bin of measured Doppler velocity. However, the effective emissivity cannot be directly de-projected in regions of z/r that are associated with ambiguous velocity bins. Such regions are particularly important when the wind velocity is confined in the vicinity of the sky plane, implying that the Doppler velocity distribution is folded on itself and confined to lower values. In such cases, the observations measure mostly the projection of the effective emissivity on the sky plane, namely the integrated brightness (or intensity). 

v) To a good approximation, the r-dependence of the de-projected effective emissivity is obtained independently from the choice of the wind configuration used for de-projection as long as the wind velocity does not depend on r. 


vi) As an illustration, we have presented two case studies of CO emission of AGB stars. They were chosen to be representative of typical observations rather than of the best space and spectral resolutions and sensitivity available today. In both cases, the method and tools developed in the present article have been shown to select efficiently wind configurations suitable for de-projection. Published analyses of the relevant observations have proposed models that are indeed favoured by these results. However, they also suggest exploring other wind configurations, of different inclination and prolateness, which may deserve being considered in the framework of a physics analysis. More importantly, they provide a deep understanding of the constraints imposed on a physics model and of how unique is the region of the (λ,φ) plane ultimately selected.

vii) The main contribution of the present work may be the insight it has provided on the issue of the under-determination of the problem of de-projection. It has underlined the importance of being conscious that a broad family of wind configurations, illustrated in blue in the right panel of Figure <ref>, can be used in principle to de-project the observed brightness data-cube into an effective emissivity at each point in space. In practice, however, constraints imposed by the requirement of axi-symmetry of the de-projected effective emissivity and of the need to populate in each pixel the totality of the observed Doppler velocity spectrum (but no more) have been found to complement each other; taken together, they are efficient in restricting the domain of the (λ,φ) plane acceptable for de-projection. A good approximation to the r-distribution of the de-projected emissivity has been obtained under the hypothesis of constant velocity but physics considerations must then be used to decide how best to combine a possible velocity gradient with such r-dependence. More generally, the constraints imposed by hydrodynamics on the relation between density, temperature and velocity will be determinant in deciding on the “best” physics model.



§ ACKNOWLEDGEMENTS

We thank Professor Pierre Lesaffre for a careful reading of the manuscript and very useful and pertinent comments that helped improving significantly its content. 
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. It also makes use of observations carried out with the IRAM NOEMA Interferometer and the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This research is funded by Graduate University of Science and Technology under grant number GUST.STS.DT2017-VL01.
Financial and/or material support from the Vietnam National Space Center, the National Foundation for Science and Technology Development (NAFOSTED), the World Laboratory and Odon Vallet fellowships is gratefully acknowledged.


          





















99
[Cami et al.2000]Cami2000
  Cami, J., Yamamura, I., de Jong, T., et al., 2000, A&A, 360, 562
[Cox et al.2012]Cox2012
  Cox, N. L. J., Kerschbaum, F., van Marle, A.-J., et al., 2012, A&A, 537, A35
[Diep et al.2016]Diep2016
  Diep, P. N., Phuong, N. T., Hoai, D. T. et al., 2016, MNRAS, 461, 4276
[Hoai et al.2014]Hoai2014
  Hoai D.T, Matthews, L.D., Winters, J.M. et al., 2014, A&A, 565, A54
[Höfner & Olofsson2018]Hofner2018
  Höfner,S. & Olofsson, H., Astron Astrophys Rev, 2018, 26, 1
[Homan et al.2018]Homan2018
  Homan, W., Richards, A., Decin, L. et al., 2018, arXiv:1804.05684v2, accepted for publication in A&A
[Le Bertre & Gérard2004]Le Bertre2004
  Le Bertre, T., & Gérard, E., 2004, A&A, 419, 549
[Le Bertre et al.2016]Le Bertre2016
  Le Bertre, T., Hoai, D.T., Nhung, P.T. and Winters, J.M., 2016, Conférence de la Société Française d'Astronomie et d'Astrophysique, L. Cambrésy et al. Eds.
[Lebzelter & Hron1999]Lebzelter1999
  Lebzelter, T. & Hron, J., 1999, A&A, 351, 553
[Libert et al.2010]Libert2010
  Libert , Y., Winters, J.M., Le Bertre, T. et al. 2010, A&A, 515, A112
[Nhung et al.2015a]Nhung2015a
  Nhung, P.T., Hoai, D.T, Winters J.M., et al., 2015a, RAA, 15, 713
[Nhung et al.2015b]Nhung2015b
  Nhung, P.T., Hoai, D.T., Winters, J.M. et al., 2015b, A&A, 583, A64 
[Nhung et al.2018]Nhung2018
  Nhung, P.T., Hoai, D.T., Tuan-Anh, P. et al., 2018, submitted for publication in RAA
[Stephenson1984]Stephenson1984
  Stephenson, C. B., 1984, Publ. Warner & Swasey Obs., 3, 1
[Winters et al.2003]Winters2003
  Winters, J.M., Le Bertre, T., Jeong, K.S. et al., 2003, A&A, 409, 715
[Winters et al.2007]Winters2007
  Winters, J.M., Le Bertre, T., Pety, J.& Neri, R., 2007, A&A, 475, 559

















