© ESO, 2017

1. Introduction

Low and intermediate mass stars evolve up the asymptotic giant branch (AGB) as they reach the end of their lives. This phase is characterised by the development of a dense and cool molecular envelope scattered with (sub-)microscopic dust particles, which the star expels into space. Recent observations of the circumstellar environments (CSEs) of AGB stars with high-resolution instruments have revealed a rich spectrum of structural complexities. These include bipolar structures (e.g. Balick et al. 2013), arcs (e.g. Decin et al. 2012; Cox et al. 2012), shells (e.g. Mauron & Huggins 2000), clumps (e.g. Bowers & Johnston 1990; Ohnaka et al. 2016; Khouri et al. 2016), spirals (e.g. Mauron & Huggins 2006; Mayer et al. 2011; Maercker et al. 2012; Kim et al. 2013), non-rotating tori (e.g. Skinner et al. 1998), and bubbles (Ramstedt et al. 2014). The most significant source of asphericity in evolved stars are thought to be magnetic fields (e.g. Pérez-Sánchez et al. 2013; van Marle et al. 2014) and binarity (Soker 1997; Huggins 2007). The importance of the latter phenomenon cannot be understated, as the multiplicity frequency of the main sequence (MS) progenitors of AGB stars is found to exceed 50 percent (Raghavan et al. 2010; Duchêne et al. 2013).

L₂ Puppis is a semi-regular pulsating variable (P = 141 days, (Kholopov et al. 1985; Bedding et al. 2005)) whose M5III spectral type corresponds to an effective temperature of T_eff ≈ 3500 K. It has radial velocity relative to the Local Standard of Rest (lsr) of v_lsr = 33.3 km s^-1. At a distance of only 64 pc (π = 15.61 ± 0.99 mas; van Leeuwen 2007) L₂ Pup is the second nearest asymptotic giant branch (AGB) star next to R Doradus (π = 18.31 ± 0.99 mas), and it is ≈30% closer than Mira. Recent continuum observations by Kervella et al. (2014, 2015) show a compact dust disk around the central star, with an inner rim of 6 AU and an outer edge of 13 AU, at an inclination of approximately 82°.

In November 2015, we obtained spatially resolved ALMA band 7 observations of the CSE of L₂ Pup using the most extended, 16-km configuration (see Sect. 2). The data shows a differentially rotating gas disk, oriented in a nearly edge-on position, adding yet another morphology to the above list. The very long baselines allowed us to reach an angular resolution of 15 × 18 mas. The high spectral resolution in the molecular data enabled us to accurately map the motion of the gas in the equatorial regions of the disk in the first paper in this series (Kervella et al. 2016). The Keplerian motion of the gas in the inner regions of the disk have been used to accurately determine the central mass of the L₂ Pup system. In addition, a significant asymmetry in the continuum emission of the east and west sides of the disk points towards the presence of a binary companion.

In the present work we model the observed emission of the ¹²CO and ¹³CO rotational transition J = 3−2 around L₂ Pup, which allows us to make predictions on the (thermo)dynamical and morphological properties of the circumstellar gas. This provides valuable constraints on the possible physical or thermo-chemical formation mechanisms producing such a morphological complexity, as well as insights on the local physics that dominate the inner wind. In addition, revealing and understanding the local inner wind dynamics sheds light on the processes that dictate the further evolution into post-AGB stars and planetary nebulae, whose morphologies have been extensively documented and have been found to be highly aspherical. Therefore, an in-depth understanding of the finer physical details governing the inner regions of complexity-harbouring AGB CSEs will aid in the quest to understand the missing morphological link between AGB stars and their progeny.

2. ALMA data description

For an extensive description of the ALMA observing set-up as well as an in-depth discussion of the continuum residuals see Kervella et al. (2016). We summarise the most relevant information here.

L₂ Pup was observed on 2015 Nov. 5 by ALMA for project code 2015.1.00141.S (PI Kervella), starting at UT 08:39:25 (epoch 2015.8448, MJD = 57 331.361) and ICRS position RA 07:13:32.475914 Dec – 44:38:17.91396. There were 45 antennas present for all or most of the observations, providing baselines from 0.09 to 16 km. However, the coverage inside 350 m is very sparse, giving a maximum angular scale of about 250 mas for reliable imaging. Standard ALMA observing and data reduction procedures were used, including phase-referencing. The flux scale is accurate to about 7% and the absolute astrometric error is about 7 mas; comparisons within this data set are much more accurate, limited only by the signal-to-noise ratio. The line-free continuum channels of L₂ Pup were used to self-calibrate the continuum and the solutions applied to all target data, and then the continuum was subtracted. Image cubes were made for the lines present in the data, and were adjusted to constant systemic velocity v_sys in the target frame.

In this paper, we present the ¹²CO and ¹³CO cubes. These have a spectral resolution of ~0.22 km s^-1 and a rms noise (off-source) of 2.5 and 2.6 mJy beam^-1, respectively. The synthesised beam size is ~, where the exact value depends on frequency. A 2.5 mas pixel size and an image size of was used.

In the following sections we describe the morphological properties of the gas emission, as deduced from velocity channel maps, and position-velocity diagrams (PVDs). We identify the patterns within the inner arcsec as consistent with the emission produced by a rotating disk of gas, and describe the fine features present in the data in detail. All descriptions relating to velocity are made relative to v_lsr = 33.3 km s^-1, i.e. the estimated systemic velocity.

The intensity maps of the ¹²CO and ¹³CO emission as a function of velocity are visually represented in Figs. A.1 and C.1, respectively. The figures are comprised of 49 panels, selected from the data cubes at even intervals in velocity space. The velocity range of the panels span a domain selected such as to best show the spatial extent of the emission.

2.1. ¹²CO J = 3−2 emission

The channel maps (Fig. A.1) reveal horizontally elongated emission features, which range up to radial angular distances of 300 mas away from the central star, and vertical heights up to 300–350 mas. These features show a very clear east-west spatial offset separation between blue and red-shifted emission. In addition, the highest (absolute) velocities probe the regions closest to the central star, while the lower (absolute) velocities probe the outer regions. These features clearly point to a compact disk with a global velocity field that has the characteristics of differential rotation, seen edge-on. This claim is further supported by making integrated-intensity maps, or stereograms, of the data (Fig. 1), where the spatial and spectral separation between both sides of the disk are clearly apparent.

The high-velocity emission of the inner disk probes only a very small vertical region. Shifting to velocities closer to the stellar velocity the extent of this vertical signal increases. This suggests that the disk is flared. However, the emission does not seem to support a stereotypical flaring pattern (i.e. an increasing rate at which the scale-height grows with radius). Rather, the ratio of the scale height over distance seems to decrease with respect to distance from the centre.

Although expected to be present, the inner rim of the gas disk is not clearly visible in the data. This is mainly due to a combination of the ALMA beam size and the fact that the tilt of the disk obscures a significant portion of the inner rim. Nevertheless, the channel maps show a rather wide initial scale height, indicating that the inner rim of the gas disk may be relatively extended (of the order of a few AU) or even inflated or puffed up (Dullemond et al. 2001; Dullemond & Dominik 2004).

The emission around central velocity appears more diffuse and extended. This emission is probably dominated by the extended CSE of L₂ Pup. However, owing to the poor uv coverage of the smallest baselines of the configuration, length scales above 250 mas have been poorly sampled. In addition, the maximum recoverable scale of the configuration is ~450 mas. Thus, some of the extended emission is probably resolved out and the remaining flux produces low-level artifacts, limiting the dynamic range in these channels. Kerschbaum & Olofsson (1999) detected ¹²CO J = 3−2 with the JCMT (beam size 13′′), with a peak of ≈59 Jy (using S(Jy)/T_mb(K) = 15.6), and estimated an expansion velocity of 3 km s^-1; their sensitivity was of order 1 Jy and the extreme-velocity wings and disk emission was not detected. Figure 3 shows spectra extracted from the ALMA data in apertures of the indicated radii. Thus, the JCMT dectected almost ten times as much emission as ALMA from the extended CSE.

Figure 3 also shows how, using a 250-mas aperture similar to the size of the disk (i.e. all the flux is recovered), there is a drop in intensity that is close to the systemic velocity, and using an aperture comparable to the size of the star itself shows absorption. This is seen in the corresponding channels of Fig. A.1 where the emission becomes very faint. This is probably due to self-absorption along the line of sight to the star by the outer regions of the disk, which must therefore be optically thick and too cool to emit significantly at 345 GHz.

On smaller scales, which are visible on both the blue- (−4.5 to −3.0 km s^-1) and red-shifted (3.0 to 5.0 km s^-1) sides, the emission is slightly brighter above and below the midplane, compared to a narrow horizontal strip passing through the star. These features may be caused by local heating of the inner disk walls directly facing the stellar radiation field. This would also explain the puffed-up inner rim. We discuss this hypothesis in Sect. 6.2.1.

The PVDs created by positioning the slit parallel to the disk plane (Fig. 2) permit the quantification of the spatial offset increase as a function of velocity. These trends have been extensively analysed by Kervella et al. (2016), resulting in a very accurate determination of the central mass of the system (0.659 ± 0.052M_⊙). The PVDs nicely exhibit the lack of compact emission around zero velocity. However, some diffuse emission can be seen around zero velocity up to large spatial offsets. This is the poorly sampled, large-scale emission discussed above. Comparing both diagrams it can be seen that the high-velocity tails do not span the same width in velocity space. The tail (of significant signal-to-noise) of the narrow slit PVD is longer by about 5 km s^-1. This suggests that the material of the inner disk has higher tangential velocities close to its equator, compared to the gas at higher vertical offsets. To quantify this velocity gradient, we constructed 11 PVDs with a vertical separation of 10 mas. We extracted the maximal velocity at each offset by measuring the length of the tails of each PVD, employing the criterion that the signal must be twice the rms noise value, outside the line. Indeed, this operation permitted us to recover the trend shown in Fig. 4. One possible cause for this vertical profile is that one may expect the upper disk regions to strongly interact with the AGB outflow, causing hydrodynamical perturbations that would impact the velocity. We discuss this further in Sect. 6.1.2.

2.2. ¹³CO J = 3−2 emission

Morphologically, the qualitative descriptions outlining the ¹³CO emission patterns (Fig C.1) are in line with the ¹²CO emission (Fig. A.1). In principle no major intrinsic differences are to be expected. The ¹³C/¹²C ratio (see e.g. Hinkle et al. (2016) for estimates in AGB stars) is likely to be reflected in the ¹²CO/¹³CO ratio, although not necessarily in a one-to-one relation because of isotope-selective photodissociation and chemical fractionation; therefore we expect the ¹³CO emission to have a lower optical depth than that of the ¹²CO emission. Hence, the absolute amount of observed emission is substantially lower in the ¹³CO emission, pushing the emission from the low-density regions of the disk below the sensitivity levels of ALMA. This is best seen in the apparent height of the ¹³CO disk, which is noticeably smaller than the ¹²CO disk (see Fig. 1). In addition, the ¹³CO disk is mostly optically thin, as it is still visible around the central velocity, contrary to the ¹²CO disk (see Fig. C.1). Combined, these effects strongly dilute and blur the finer details present in the ¹²CO emission maps, resulting in a much smoother emission distribution.

There is, however, a strong asymmetry between the low-velocity blue-shifted side compared to the low-velocity red-shifted side. This asymmetry is best seen by comparing the emission in the −3.5 to −2.0 km s^-1 channels with the 1.0 to 3.5 km s^-1 channels. In these channels, the blue-shifted emission conforms with the expected patterns of a rotating disk. However, the red-shifted emission seems to be substantially perturbed, causing a deviation from the expected emission morphology. Emission around the equatorial regions does not seem to be present, while an emission tendril is seen to reach out to the north-east. This feature subsides beyond 3.5 km s^-1. Although one would expect the presence of a binary companion (Kervella et al. 2016) to generate local perturbations, the projected velocities at which the perturbation is seen are too low to be related to the companion system, which, due to its proximity to the central star, can be expected to only substantially affect the highest velocities. We therefore interpret this feature as a hydrodynamical inhomogeneity in the outer boundaries of the disk.

The ¹³CO PVDs (Fig. 5) show a high degree of similarity with the ¹²CO diagrams. In fact, the only clear difference is the length of the high-velocity tail, which seems to be identical for both the narrow-slit and the wide-slit PVDs. This can again be attributed to the molecular abundance fraction. The confined vertical extent of the ¹³CO emission means it mainly probes the Keplerian equatorial regions. The high-velocity tails of the PVDs of the ¹³CO emission are thus mostly insensitive to the sub-Keplerian velocities at higher vertical spatial offsets.

3. Computational methods

In this section we briefly describe the computational methods and strategy used to model the L₂ Pup CO line emission. We perform 3D radiative transfer via the LIME code. The result of the radiative transfer calculations are subsequently subjected to the simulation observation algorithms of CASA to generate 3D synthetic data that can be directly compared to the observed L₂ Pup data.

3.1. LIME

We used the full 3D, non-local thermodynamical equilibrium (NLTE), infrared (IR) and submillimeter radiative transfer (RT) code LIME for the radiation transport calculations; an in-depth description of the methods and algorithms utilised by the code can be found in Brinch & Hogerheijde (2010). The model is sampled by 10⁵ grid points, of which half are distributed logarithmically (gradually refining the grid towards its centre), and half are distributed randomly. In addition, the position of the grid points is weighted by relative density, further increasing the mesh refinement in the high-density regions of the model. Another 10⁴ grid points are positioned at the edge of the numerical domain, representing the points where the RT calculations are finalised. This total of 11 × 10⁵ grid points are Delaunay-triangulated to form approximately 7 × 10⁵ tetrahedral cells, which are subsequently Voronoi-tesselated and inside which the physical conditions are assumed to be constant. The physical set-up is described in terms of density, temperature, molecular abundance, macro- and micro-scale velocity fields. In order to detemine the mean intensity field, the RT equations are solved along the Delaunay lines connecting the neighbouring grid points up to the edge of the numerical domain. Once the RT calculations are finalised and the level populations have converged, the simulation is ray traced from a specific vantage point. It is important to note that the code solves the RT equations exclusively in order to converge the level populations for a given initial physical set-up. This user-defined physical set-up is thus not modified as the level populations converge. This, in turn, means that self-consistency of the physical model is only guaranteed if it was verified prior to the application of LIME. A final remark is that LIME was developed primarily for the modelling of molecular emission. The manner in which dust is taken into account in the calculations is rather rudimentary, and only suitable for first order dust modelling. However, because the CO molecule has such a low electric dipole moment, the contributions of the diffuse radiation field (from e.g. the dust) to the level excitation are not very important, diminishing the available free parameters and thus the modelling difficulty. The spectroscopic CO data of the LAMDA database (Schöier et al. 2005) were used; the collisional rates were taken from Yang et al. (2010). Because of its low electric dipole moment, the surrounding radiation field is not expected to affect the CO emission much. Nevertheless, both dust and the stellar radiation field were taken into account in the radiative-transfer calculations. The dust composition has been adopted from Kervella et al. (2015), with the dust density distribution following the gas density distribution and a gas-to-dust mass ratio of 100. The stellar radiation field is approximated by a black body with a temperature of 3500 K and a luminosity of 2500 L_⊙.

3.2. CASA

After retrieving the intrinsic emission distribution from LIME, we post-processed it with CASA (McMullin et al. 2007) to simulate an observation with ALMA. The actual general observation conditions and instrumental set-up were adopted as input parameters for the simulations, and are shown in Table 1. This also implies we used the same antenna locations to give consistent visibility plane coverage between the observations and simulations.

3.3. Modelling strategy

We wish to make it clear that our best-fit models are not achieved via the exploration and statistical analysis of a grid of models covering the more than 20 free parameters available to us. Such an analysis cannot simply be performed with present-day computational facilities. The best-fit model we obtained was created by making reasonable assumptions and educated guesses on the general physical properties (based on previous studies and the nature of the observed emission patterns), and subsequently attempting to converge to the reproduction of the emission characteristics via iterative adjustment of these parameters. Thus, when we use the term “best fit” in this paper, we do not refer to a statistically significant result that follows from extensive and detailed post-processing, but rather to the set of parameters for which the simulated emission shows a good visual resemblance to the emission patterns obtained by ALMA.

Owing to these computational limitations, we devised a strategy for the determination of the physical disk properties via RT modelling. We decided to start by focusing the analysis on the ¹²CO emission. The reasons for this which all stem from the expected lower ¹³CO molecular abundance compared to the ¹²CO abundance (Hinkle et al. 2016), are listed below:

• ¹³CO over H₂ abundance ratios are not well known, while ¹²CO over H₂ ratios are better constrained.

• The lower abundance of ¹³CO results in a flatter and dimmer disk, in which the high degree of complexity in the thermodynamical, kinematical and morphological properties of the disk are less pronounced.

• Photodissociation effects, if present, have a larger impact on the overall morphology of the ¹³CO disk. We make no attempt to quantify the magnitude of photodissociation effects on the gas, so using the ¹²CO emission as a starting point for the modelling drastically reduces this systematic error.

• We can expect the ¹³CO emission to be optically thin, causing temperature and absolute density to be fully degenerate. Using the optically thick ¹²CO emission untangles both the temperature and the denisty by strongly reducing the influence of the absolute density.

After we determined the set of physical conditions that best reproduce the ¹²CO emission patterns, we refined the model using the ¹³CO emission, making sure to retain consistency with the ¹²CO emission.

The extended emission around central velocity (visible in Fig. A.1, discussed in Sect. 2.1) is probably not part of the disk system, and we therefore do not include it in the modelling effort.

4. Physical structure of the gas disk

4.1. Dominant and turbulent velocity field

A first analysis of the molecular emission in the ALMA data of L₂ Pup by Kervella et al. (2016) has shown that the differential rotation of the equatorial regions of the inner gas disk (within a radius of 6 AU) follows Keplerian rotation. This analysis was primarily based on the ²⁹SiO emission, which has the highest signal-to-noise value of the data set. The tangential velocities of the inner disk could be used to estimate the central mass of the system, which was found to be 0.66 ± 0.05M_⊙. Beyond this Keplerian region (in the outer disk, i.e. for radii greater than 6 AU) the velocity field sharply transitions to a sub-Keplerian regime, where the radial velocity decrease is proportional to r^-0.85. The radius at which the tangential velocity transitions from Keplerian to sub-Keplerian seems to correspond to the inner rim of the dust disk, as determined by Kervella et al. (2015), indicating that the dust disk may play an important kinematical role in the steady-state hydrodynamical disk system. Following these findings, we separated our velocity model into two regimes, the inner disk (r < 6 AU), and the outer disk (r> 6 AU).

4.1.1. Inner disk

The tangential gas speeds follow Keplerian rotation: ${\mathit{v}}_{{\mathit{r}}_{\mathit{xy}}\mathit{,}\mathrm{inner}}\mathrm{=}\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{r}}_{\mathit{xy}}}}\mathit{,}$(1)where ${\mathit{r}}_{\mathit{xy}}\mathrm{=}\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}$ is the cylindrical radial distance coordinate, G is the universal gravitational constant, and M_∗ is the mass of the central object.

We approximated the vertical velocity profile (Fig. 4) in the inner disk. We define the polar angle as $\mathit{\theta }\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}\mathrm{=}\arccos \left(\frac{\mathit{z}}{\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}}\right)\mathit{,}$(2)assuming the coordinates (x,y,z) make up an ortho-normal set of axes, and the equatorial plane of the disk coincides with the (x,y) plane. To reproduce the vertical velocity structure of the inner disk (Fig. 4), the equatorial tangential velocity field of the inner disk is multiplied with a factor ζ, expressed as $\mathit{\zeta }\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{1}}\mathrm{)}\mathrm{+}\left(\frac{{\mathit{f}}_{\mathrm{1}}}{{\mathrm{e}}^{{\mathit{s}}_{\mathrm{1}}\mathrm{(}\mathit{\theta }\mathrm{-}\mathrm{(}\mathit{\pi }\mathrm{-}\mathit{\delta }\mathrm{\left)}\mathrm{\right)}}\mathrm{+}\mathrm{1}}\right)\mathrm{-}\left(\frac{{\mathit{f}}_{\mathrm{1}}}{{\mathrm{e}}^{{\mathit{s}}_{\mathrm{1}}\mathrm{(}\mathit{\theta }\mathrm{-}\mathit{\delta }\mathrm{)}}\mathrm{+}\mathrm{1}}\right)\mathit{,}$(3)where f₁ is the factor by which the velocity drops vertically, s₁ is the slope steepness of the drop, and δ is the polar angle at which the drop occurs. The best-fit parameters are summarised in Table 2. This function ζ is visually represented in Fig. 6. We hypothesise about the physical origin for the ζ function in the Sect. 6.1.2.

4.1.2. Outer disk

A second (sub-Keplerian) velocity field for the equatorial tangential velocities in the outer gas disk (r> 6 AU) is also based on the preliminary findings of Kervella et al. (2016) and is expressed as ${\mathit{v}}_{{\mathit{r}}_{\mathit{xy}}\mathit{,}\mathrm{outer}}\mathrm{=}\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{\mathrm{6}\mathrm{AU}}}{\left(\frac{\mathrm{6}\mathrm{AU}}{{\mathit{r}}_{\mathit{xy}}}\right)}^{\mathrm{0.85}}\mathrm{·}$(4)We found that introducing an additional vertical velocity profile in the outer disk better reproduces the outer edges at high spatial offsets of the wide-slit PVD in Fig. 2. Equation (4) can thus be improved by multiplication with a factor α, expressed as $\mathit{\alpha }\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{2}}\mathrm{)}\mathrm{+}{\mathit{f}}_{\mathrm{2}}\exp \left(\frac{\mathrm{-}{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{s}}_{\mathrm{2}}^{\mathrm{2}}}\right)\mathrm{·}$(5)In this expression f₂ represents the percentage factor by which the equatorial tangential velocity is reduced and s₂ is a measure for the rate at which the velocity drops to (1−f₂)v_rxy,outer. The parameter values found to best reproduce the emission patterns are listed in Table 3.

The α factor simply reduces the tangential velocities of the regions with high vertical offsets by 10%. The continuous transition from the equatorial tangential velocities to the reduced tangential velocities follows a Gaussian trend, with a characteristic width of 1.5 AU. This means that the tangential velocities are fully reduced to 90% of their equatorial values at a vertical distance of around 3 AU from the equator.

4.1.3. Small-scale random motions

We used the emission gap around the systemic velocity to determine the velocity field of the small-scale random motion (SSRM). Assuming the gap is indeed created by the absorption of the ¹²CO emission by the cold outer regions of the disk, the width of the gap can be used as a proxy for the SSRM. The emission gap spans a width of 1.7 km s^-1. Assuming such motions to follow a Gaussian distribution, such a gap thus corresponds to a SSRM standard deviation of v_ssrm = 0.55 km s^-1. Because we have no means of probing the SSRM for the regions closer to the star, we assume it to be constant throughout space. However, one can expect the SSRM to increase more deeply inside the disk. We discuss this further in Sect. 6.1.2.

4.2. H₂ density constraints

The morphology of the disk is properly captured by assuming the following expression: $\mathit{\rho }\mathrm{\left(}{\mathit{r}}_{\mathit{xy}}\mathit{,\phi ,z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\left(\frac{{\mathit{r}}_{\mathit{xy}}}{{\mathit{r}}_{\mathrm{c}}}\right)}^{\mathit{p}}\exp \left[\frac{\mathrm{-}{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}\mathit{H}\mathrm{(}{\mathit{r}}_{\mathit{xy}}{\mathrm{)}}^{\mathrm{2}}}\right]\mathit{,}$(6)with r_c the inner rim of the gas disk (in the equatorial plane), ρ₀ the H₂ gas density at r_c, p the rate with which the density subsides radially, and H(r_xy) the Gaussian scale height of the disk, expressed as $\mathit{H}\mathrm{\left(}{\mathit{r}}_{\mathit{xy}}\mathrm{\right)}\mathrm{=}{\mathit{H}}_{\mathit{c}}{\left(\frac{{\mathit{r}}_{\mathit{xy}}}{{\mathit{r}}_{\mathrm{c}}}\right)}^{\mathit{h}}\mathit{,}$(7)where H_c is the initial scale height at the location of the inner rim in the midplane of the disk, and h a measure for the rate at which the scale height increases. The radial distance variable r_xy is expressed as the cylindrical radial coordinate $\sqrt{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}}$. This density expression follows from the solution of a thin, vertically isothermal, non-self-gravitating disk in hydrostatic equilibrium. Although the data indicate that the density profile is more complex, this simplified expression for the profile is sufficiently accurate for the modelling. We reflect upon this in Sect. 6.3.2.

Kervella et al. (2016) have shown that the gas in the inner disk (<6 AU) is undergoing Keplerian rotation. This allows us to make an estimate of the radius of the inner rim using the largest observed velocities in the disk signal. The data shows the faintest beginnings of a disk signal approximately 17 km s^-1 away from the central velocity. Assuming a central mass of 0.659 M_⊙ (Kervella et al. 2016), this tangential velocity translates to the inner rim of the gas being located at approximately r_c = 2.0 AU.

The evolution of the height of the gas disk as a function of radial distance in the disk plane could be readily deduced from the shape of the disk emission in the low-velocity channel maps (Fig. A.1). They show that the disk has a degree of flaring. However, contrary to the archetypal appearance of a flaring disk, the rate of flaring of the disk of L₂ Pup decreases as a function of distance from the central star. This trait is reflected in the value of the parameter h, which has been determined to be less than unity. Finally, combining the location of the inner rim with the flaring rate, the height of the emission in each channel is best reproduced for an initial scale height of 1.5 AU.

The radial behaviour of the density is captured in the parameters ρ₀ and p. The value of ρ₀ is empirically determined by the iterative efforts targeting the convergence of model versus ¹²CO data emission values in the innermost regions of the disk. The parameter p determines the rate at which the radial density drops, but this parameter cannot be determined from modelling the ¹²CO emission. This is due to the overall optically thick nature of the disk, which results in a substantial insensitivity of the radiative transfer model to the parameter p. The optically thin ¹³CO emission permits the evaluation of this parameter p. An extensive exploration of the p parameter space has revealed that the data is best modelled by simply assuming a radial power law drop characterised by p = −3.1.

The best-fit parameters, which follow from the above-mentioned considerations, are listed in Table 4.

4.3. Disk temperature

The ¹²CO disk must also be optically thick in its more inner regions because it is optically thick at the disk’s outer edges. Thus, the ¹²CO emission patterns are mainly produced on the surface boundary confining the gas emitting at a certain frequency along the line of sight. The emission produced by the bulk of the gas lies hidden behind this surface. Therefore, the ¹²CO emission probes the temperature of the gas in this surface layer very accurately and is mostly insensitive to the distribution of mass. Hence, because of the extremely high spectral resolution of the data, these surfaces are only spatially separated by a small distance, and the temperature can therefore be accurately and sufficiently sampled all throughout the disk, permitting the full temperature profile to be precisely determined from the ¹²CO emission.

4.3.1. Radial profile

As seen in Fig. A.1, the ¹²CO emission close to the central star is brightest, indicating relatively high gas temperatures. The emission drops substantially within the first hundred mas, after which it drops much slower, suddenly plummeting as the outer edge of the emission is reached. This suggests a three-staged temperature profile.

Analytically, the temperature profile producing the best-fit model is expressed as $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{T}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{(}{\mathit{T}}_{\mathrm{z}}\mathrm{-}{\mathit{T}}_{\mathrm{p}}\mathrm{)}\exp \left(\mathrm{-}\frac{{\mathit{r}}_{\mathrm{xy}}^{\mathrm{2}}}{\mathrm{2}{\mathit{w}}_{\mathrm{1}}^{\mathrm{2}}}\right)\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{-}\left(\frac{{\mathit{T}}_{\mathrm{p}}}{\mathit{\pi }}\right)\arctan \left(\frac{\mathit{r}\mathrm{-}\mathit{D}}{{\mathit{w}}_{\mathrm{2}}}\mathrm{-}\frac{\mathit{\pi }}{\mathrm{2}}\right)\mathit{,}\end{array}& & \end{array}$(8)where T_z is the temperature at the inner rim of the disk, T_p is the temperature of the plateau, w₁ is a measurement for the rate of the first temperature decrease, D is the distance at which the second major drop in temperature is located, and w₂ is a measure for the rate of the second temperature drop. Table 5 lists the parameter values that were been found to produce the best fit of the data. The temperature profile is shown in logarithmic scales in Fig. 7. The temperature levels-off to a near-constant value of 500 K at a distance of approximately 6 AU, which corresponds with the inner rim of the dust disk observed by Kervella et al. (2015). This indicates that the dust disk plays an important thermodynamical role in the disk system. The temperature of the inner walls of the gas disk (facing the direct radiation of the stellar photosphere) was determined to be around 2500 K. We discuss this further in Sects. 6.2.1 and 6.2.2.

The obtained radial temperature expression can also be approximated by a sequence of power laws, expressed as T(r) = T_i(r/r_i)^ϵ. Here, T_i stands for the temperature at distance r_i, and ϵ is the slope of the power law. Table 5 lists the power law sequence that best fits the continuous model (shown as black dotted lines in Fig. 7).

4.3.2. Vertical profile

The high angular resolution of the data allows the identification of some structural complexity in the inner disk. As discussed in Sect. 3.2.2, the channel maps seem to exhibit two patches of increased emission in the inner disk both on the blue- and the red-shifted sides. We interpret this signal as the local heating of the walls of the inner disk facing the central star. We modelled this feature as follows: ${\mathit{T}}_{\mathrm{z}}\mathrm{=}{\mathit{T}}_{\mathrm{w}}\mathrm{-}\mathrm{(}{\mathit{T}}_{\mathrm{w}}\mathrm{-}{\mathit{T}}_{\mathrm{eq}}\mathrm{)}\exp \left(\frac{\mathrm{-}{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{T}}^{\mathrm{2}}}\right)\mathit{,}$(9)where T_w is the temperature of the star-facing disk wall, T_eq is the temperature of the equatorial disk plane, and σ_T is the Gaussian one-sigma size of the cool equatorial disk plane. In order to reproduce the correct vertical emission contrast in the inner disk we need the temperature of the central disk plane to be approximately at a temperature of 1000 K. Table 5 lists the parameters that were been found to best fit the data.

4.4. ¹²CO and ¹³CO molecular abundance

As stated before, molecular abundance and absolute density are degenerate. We opted not to express the ¹²CO abundance as a function of the spatial coordinates, but rather as a constant value. We found no studies that explicitly focus on the molecular abundance of ¹²CO in L₂ Pup have been found in the literature. Additionally, the handful of studies describing the general parameters of disks around evolved stars indicate the ¹²CO abundance to be around 10^-4 (e.g. Bujarrabal et al. 2013; Bujarrabal & Alcolea 2013). Studies of ¹²CO abundance values inside protoplanetary disks around young stars also seem to reproduce this value (e.g. France et al. 2014). In addition, because of its high binding energy, carbon monoxide is generally almost chemically unreactive (e.g. Gobrecht et al. 2016).

We have thus chosen to set the ¹²CO/H₂ ratio to the constant value of 10^-4 for this model; this value which lies within the reasonable range of expected ¹²CO abundance values for AGB stars (e.g. Decin et al. 2010; De Beck et al. 2012).

The ¹³CO abundance is easily acquired from the values of the physical properties of the inner rim determined from the ¹²CO analysis as boundary condition. The ¹³CO molecular abundance was iteratively adapted until the emission of the inner rim of the ¹³CO model matched the data. This resulted in an abundance estimate of 10^-5, resulting in a ¹²CO over ¹³CO abundance of 10.

5. Radiative transfer results

The internal temperature and velocity structure of the disk around L₂ Pup is visually represented by the schematic diagram shown in Fig. 8. This physical set-up is propagated through the RT code LIME, which produces velocity-dependent maps of intrinsic line emission from the species in question. The CASA observation simulation scripts transform the intrinsic maps, by introducing observational and instrumental effects, so they can be compared one-to-one with the ALMA data. We simulated the observables for an almost edge-on appearance of the disk, at an inclination of 82 degrees.

These predicted channel maps are represented visually in Figs. B.1 and D.1. Here one can see that the overall morphology and general features present in the ALMA data have been well reproduced. One must keep in mind that all emission in the ALMA data which is smooth on scales greater than 250 mas has been insufficiently sampled in uv space, generating the (seemingly) randomly scattered emission patches that are most strongly present around systemic velocity in the ¹²CO maps. The bulk of this emission probably makes no part of the disk system, and it is therefore not present in the model data.

It is quite difficult to compare the model and data channel maps one-to-one. We shall thus focus the rest of the comparative discussion on the PVDs, which are shown in Figs. 9–11. Figure 9 nicely illustrates the quality of the ¹²CO model along zero spatial offset, which corresponds to the disk midplane. All major and minor features present in the data have been reproduced. We succeeded in replicating the edges of the emission signal in velocity space, the absolute velocity magnitude of the offset peaks, the width(in velocity space) of the emission gap around the central velocity and the internal emission distribution of the PVD. The curvature of the signal, which characterises how the outer edge of the disk emission varies as a function of velocity has also been modelled properly. Also the absolute emission values have been accurately reproduced. The major difference is found when comparing the signal near the velocities closest to v_sys. The outer regions of the disk are probably quite turbulent on different length scales, causing the offset of the lowest velocity regions to decline in a more irregular fashion compared to the analytical model.

The limitations of the model become fully apparent in the full disk-width slit PVDs, shown in Fig. 10, where important discrepancies between the data and the model appear. The absolute level of emission around zero spatial offset is slightly overestimated by the model. This is probably caused by the introduction of a vertical temperature dependence to model the temperature difference between the walls of the disk and its equatorial regions. The model most probably overestimates the size of the hottest regions, contributing to the overestimated average emission around zero offset. There is also an emission distribution asymmetry between the east and west wings of the data, which is of unknown origin, but can perhaps be attributed to the presence of a binary companion, as suggested by Kervella et al. (2016). Close to the central velocities one can clearly identify emission in the data that is not present in the model. We believe this emission originates completely from the poorly reconstructed large-scale emission from the CSE outside the disk, observed by the (poorly covered) shortest baselines of the configuration.

The final PVD comparison we consider is the full disk-width slit PVD of the ¹³CO data, shown in Fig. 11. We do not focus on the narrow-slit PVD of the ¹³CO emission as it is virtually identical to the discussion on the ¹²CO narrow-slit diagram. The emission distribution of the wide-slit PVD has been well reproduced by the model, and so has the absolute level of the emission. However, a clear difference is seen in the extent and location of the maximum spatial offset peaks. The peaks of the model are very pointy and long, and are not as widely separated in velocity space as the offset peaks in the data. We believe this difference can be attributed to photodissociation effects, which could affect the outermost regions of the ¹³CO disk. Accounting for this effect could thus seriously reduce the maximal spatial extent of the offset peaks in the model PVD. In addition, the location of the peaks in velocity space would then also shift to higher projected velocities.

6. Discussion

6.1. Parameter sensitivity

Because we cannot perform an extensive and in-depth exploration of the parameter space available to us (see Sect. 3.3), we cannot quantify the quality of the model in terms of statistically significant errors. Hence, we dedicate this section to a discussion of the uncertainties in the estimated parameter values.

The kinematical properties of the disk can be directly measured from the data. We are thus confident that both the equatorial and vertical velocity profile are accurately constrained and the parameter values that have been deduced directly from the velocity field can be trusted. Furthermore, the geometrical properties of the system can also be directly extracted from the data. We do not expect major uncertainties here either.

Carbon monoxide has a very low electric dipole moment causing the density, temperature, and molecular fractional abundance (DTA) parameters to be completely degenerate. Fortunately, the modelled data consisted of both optically thick and thin emission. In case of the optically thin ¹³CO lines, one cannot pin down the exact combination of temperature and density properties of the disk and the abundance of this molecule. In case of the optically thick ¹² lines the degeneracy of temperature and density profile is mostly lifted, and the temperature profile can be determined. However, the assumed value for the ¹²CO molecular abundance remains an uncertainty. The adopted value, from reasonable and justifiable arguments, directly and linearly affects the disk density. For example, an increase in the molecular abundance by a factor 2 should result in a decrease in density by the same factor. Hence, the largest uncetainty relates to the determined value for ρ₀. But, as argued in Sect. 4.4, we do not expect the true abundance fraction to deviate substantially from the assumed one. The radial density power law should not be affected by the abundance estimate, as it is mostly dependent on the assumed temperature profile.

The optical thickness of the ¹²CO emission, and the high spectral resolution of the data, has permitted us to very accurately constrain the equatorial temperature profile. And assuming the adopted ¹²CO abundance fraction is correct, then the retrieved equatorial density parameters and ¹³CO abundance value are also very tightly constrained. However, where the ¹²CO emission becomes optically thin, the subsequently constrained parameter values likely have much larger uncertainties. In the inner disk, the low density regions of high vertical offset are cut off by the shearing of the stellar outflow. Here the ¹²CO emission does not become optically thin, and the free parameters are thus expected to be tightly constrained. At the top and bottom edges of the outer disk, the ¹²CO emission definitely turns optically thin. This locally reinstates the full DTA degeneracy, therefore breaking down the assumption that the temperature can be directly measured from the emission. However, the regions in question are small compared to the full disk size. Thus, increased local uncertainties on the derived parameter values probably does not affect the final result by too great an amount.

6.2. Velocity field

6.2.1. Sub-Keplerian regime of the outer disk

The sub-Keplerian regime of the disk has a power law dependence of the tangential velocity that is proportional to r^-0.85 (Kervella et al. 2016). This suggests the presence of a mechanism that slows down the gas and whose effectiveness increases with radius. One explanation for such a deviation from Keplerian rotation is a radial force imbalance. However, in this scenario the sub-Keplerian nature of the outer disk should result in the appearance of an inward motion of the gravitationally attracted material. Such an inward velocity component has not been detected in the emission (see Homan et al. 2016, for an overview of the expected emission patterns caused by such an additional inward velocity component). This means that either the projected inward velocity component is negligibly small compared to the spectral resolution, or a mechanism exists that generates an additional outward force which balances the gravitational pull, yet increases in strength with distance. A possible origin for such a force would be the radial stratification of dust grain sizes within the disk. Such a dependence, with larger grains closer to the star, can produce a radial radiation pressure profile that is able to alter the effective gravity felt by the dust grains throughout the disk (see upcoming paper). Assuming a tight mechanical coupling between dust and gas, this scenario should then act as a cumulative reduction of the effective gravity felt by the surrounding gas as a function of radial distance.

6.2.2. Vertical velocity profile of the inner disk: the ζ and α functions.

We hypothesise that the ζ profile (Fig. 6), detailing the vertical distribution of tangential velocities in the inner disk, originates from hydrodynamical boundary instabilities at the disk/radial outflow interface. Media subjected to velocity shear will develop Kelvin-Helmholtz instabilities (KHI; Funada & Joseph 2001), which destroy any original organised velocity field, causing turbulent features of high complexity to arise. For a physical set-up like the inner disk of L₂ Pup, the timescale for KHI to arise is of the order of years (Funada & Joseph 2001). The unstable regions of the system will converge to a maximum entropy state via mixing when left to operate for a long time (about 10³ × KHI timescale). For the L₂ Pup system this would be of the order of thousands of years, which is much shorter than the AGB lifetime. The resulting velocity field inside this mixed region is an averaged distribution over all velocities that were part of the interaction. Thus, a radial outflow (with no tangential velocity component) shearing against a rotating gas disk creates KHI, which average out to an effective decreased rotation speed. Similarly, one can expect a small radial velocity component to be present in this “boundary region”. However, because the radial component would be small, and because of the orientation of the radial field with respect to the observer, the effect would potentially only be visible in the channel maps close to central velocity. The emission in these channel maps is hidden behind the optically thick, cold outer disk edge. Thus, we are unable to confirm the presence of this small radial velocity component in the regions with a high vertical offset of the inner disk.

The origins of the α profile, which modifies the vertical profile of the rotational velocity in the outer disk, can be attributed to the same effect. However, because the relative velocities between the disk and the outflow are smaller in the outer disk compared to the inner disk, the magnitude of the effect is reduced.

6.2.3. Velocity field of the small-scale random motion

We derived the magnitude of v_ssrm from the (lack of) spectral signatures of the gas located in the outermost regions of the disk. We do not possess a complete description of how one can expect the SSRM of a medium to depend on the physical properties of that medium. Therefore, to limit the number of assumptions made, we approximated the SSRM as constant throughout the disk. However, this assumption is rather unlikely. The rms thermal speed of molecules at a certain temperature is given by ${\mathit{v}}_{\mathrm{th}\mathit{,}\mathrm{rms}}\mathrm{=}\sqrt{\frac{\mathrm{3}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{m}}}\mathit{,}$(10)where k_B is the Boltzmann constant, T is the temperature of the gas, and m is the molecular weight. The rms thermal speed reaches values up to 1.5 km s^-1 for CO molecules at 2500 K, which is significantly greater than the v_ssrm = 0.5 km s^-1 value used in this paper. In addition, as can be inferred from Fig. 12, the thermal velocity in the cool outer regions only makes up about 50% of the deduced SSRM. This is a strong indication that the small-scale velocity in the disk consists of additional components besides the thermal rms speed and that our assumed v_ssrm probably underestimates its true physical value inside the disk.

6.3. Temperature structure

6.3.1. Vertical temperature profile of the inner disk

The modelling strongly suggests the requirement of a vertical temperature profile describing the extent by which the inner disk walls facing the direct stellar radiation field are heated. The maximum temperature of the gas in these regions was determined to be approximately 2500 K. This can be readily explained by thermal coupling of the gas with the stellar radiation field. The disk walls directly facing the stellar radiation heat up to high temperatures, while shielding the equatorial regions of the disk, keeping them at cooler temperatures.

6.3.2. Radial temperature profile

The shape of the radial temperature profile (Fig. 7) seems rather peculiar. However, the regime transition between the inner disk and the outer disk seems to strongly indicate that the presence of dust plays an important thermodynamical role. We hypothesise the following explanation. In the inner disk, the temperature profile is simply set by the dissipation of heat, provided by the photospheric resonance photons. From 6 AU onward, dust is prominently present in the disk. Thus, from this radius outward a substantial increase in sensitivity to the photospheric continuum emission exists. Assuming a tight mechanical coupling between dust and gas (which the velocity profile suggests), the dust provides a continuous source of heat from inside the disk. This lessens the radial rate of temperature decrease and acts as a thermal buffer. This happens for as long as an interaction between the stellar radiation field and the dust in the disk occurs. However, finally, at a specific radius one can expect the dust to become fully optically thick (radially), prohibiting the further propagation of energy via photospheric continuum photons, causing a final temperature drop at large radii.

6.4. Density structure

6.4.1. Disk mass

From the density structure of the disk we may calculate its mass. Integrating over an infinite volume the mass of the disk is found to be approximately 2.2 × 10^-4M_⊙, or about 0.23 Jupiter masses. We analysed the sensitivity of this value to the determined density parameters. Assuming each parameter would have its true physical value between plus and minus 25% of the determined value, then the disk mass would lie between 0.65 × 10^-4M_⊙ and 9.7 × 10^-4M_⊙. The determined disk mass is most sensitive to the exponent of the radial power law. If we assume the same spread on the parameter values as before, but now keep the radial power fixed, the spread on the disk mass diminishes significantly. It would lie between 1.15 × 10^-4M_⊙ and 3.4 × 10^-4M_⊙. These values are in agreement with the estimates of Kervella et al. (2014, 2015).

6.4.2. Density profile

As stated in Sect. 4.3, the density description conforms with the solution of a thin, vertically isothermal, non-self-gravitating disk in hydrostatic equilibrium. Although on first sight the disk does not seem to be congruent with these assumptions, a deeper consideration of the deduced disk physics can prove otherwise. The channel maps (Figs. A.1 and C.1) exhibit a rather broad disk with a height of around 300 mas. We find that the observations are well fitted by a disk with a Gaussian scale height (Sect. 4.3.2), which means that the bulk of the disk is in fact located in a thinner slab around the equator. Around 75% of the total disk material is located in a 5 AU (78 mas) wide slab. The mass of the disk is approximately 1.5% of the central mass, making its self-gravity negligible with respect to the primary gravitational source. The largest part of the disk has been found to be vertically isothermal. This property has only been found to break down in the star-facing wall of the inner disk. The resultant features, however, are only present in the thin interaction region at the disk-outflow interface. The effect of this local heating on the overall thermal disk structure is negligible. In fact, this vertical temperature profile of the inner disk could be the cause for the h = 0.2 radial scale height growth to be lower than expected (<1). As for vertical hydrostatic equilibrium (VHE), this is very difficult to deduce from the observations. One can only assume that if the disk is indeed a stable structure present around the AGB star, then VHE must be satisfied. It therefore seems that the assumptions contained within the analytical density description may very well be satisfied.

The value p = −3.1, which characterises the slope of the radial density power law, seems very steep compared to the typical slopes presented in the literature. Typical density slope values for disks around post-AGB stars and planetary nebulae are of the order of p = −2.0 (e.g. Bujarrabal & Alcolea 2013; Bujarrabal et al. 2015, 2016). Protoplanetary disks seem to possess even flatter profiles with −2.0 ≤ p ≤ −1.0 (e.g. Hartmann et al. 1998; Zhang et al. 2014; Podio et al. 2014; Kirchschlager et al. 2016). Rings of asteroids in dynamic resonance with a planet, such as in the asteroid belt of our solar system, have an even lower value of p ~ −0.1 (Ryan et al. 2009).

6.5. Evolutionary considerations

6.5.1. Disk origin

The scenario that the disk is a remnant from an earlier evolutionary phase is improbable as protoplanetary disks disperse fairly early on in the main sequence phase (e.g. Kimura et al. 2016). Therefore, the disk must have been formed during the later evolutionary stages of the star. A scenario for a recent formation of a disk may be that of a latitude dependent mass loss in the AGB phase. However, it seems highly unlikely that such a mechanism produces a disk system that is as vertically compact as observed here. Moreover, it needs a mechanism to prevent the material escaping from the system.

One of the most plausible scenarios for the formation of an equatorial density enhancement in the CSE of an AGB star is by means of a binary companion (Theuns & Jorissen 1993; Soker 1994; Mastrodemos & Morris 1999; Kim 2011; Mohamed & Podsiadlowski 2012). Such a companion funnels part of the wind material in its gravitational field, expelling it as a compact stream along the orbital plane. Kervella et al. (2014, 2015, 2016) have suggested the presence of a binary companion around L₂ Pup, based on precise measurements of continuum asymmetries. The fact that the inner rim of the gas disk (determined in this work) coincides with the location of the observed candidate companion of L₂ Pup in literature supports this hypothesis.

The presence of a binary companion is also supported by angular momentum considerations. From the density distribution and velocity structure we calculate the total angular momentum contained within the gas disk to be ${\mathit{L}}_{\mathrm{tot}}\mathrm{=}\mathrm{3.1}\mathrm{\times }{\mathrm{10}}^{\mathrm{42}}{\mathrm{m}}^{\mathrm{2}}\mathrm{kg}\mathit{/}\mathrm{s}\mathit{.}$(11)A first analysis by Kervella et al. (2016) suggests that the star L₂ Pup had an initial mass comparable to the mass of the Sun. Recent work has shown the Sun to contain approximately 10⁴¹ m²kg/s of angular momentum (Iorio 2012). In addition, recent astroseismologic analyses of light curves of solar-type stars indicate only minor differences in internal rotation structure (Benomar et al. 2015; Schunker et al. 2016). Thus, making the assumption that the total angular momentum of L₂ Pup on the main sequence was comparable to that of the Sun, we find that the disk around L₂ Pup contains approximately 30 times the total amount of available angular momentum inside the star. The presence of a binary companion could easily explain this angular momentum discrepancy. By assuming the companion follows Keplerian rotation, and by assuming an angular momentum transfer efficiency of 30% (Nordhaus & Blackman 2006; Nordhaus et al. 2007) from the companion to its surroundings, the total angular momentum in the disk could be explained by the presence of a 1 Jupiter mass companion. Extending the “sensitivity analysis” on the mass of the disk to this discussion (assuming the same uncertainty range of the determined density parameters and assuming zero uncertainty on the velocity field), the mass of the companion could be expected to lie within the 0.25 to 5.5 Jupiter mass range. This estimation is in agreement with the mass range of the companion proposed by Kervella et al. (2016). These estimates are lower limits, as one can expect a portion of the main sequence angular momentum to have been lost in the stellar wind during previous evolutionary stages.

In Paper I, Kervella et al. (2016) discuss the likelihood of the companion planet to have formed inside the disk. After all, the disk contains enough material to build one or more planets. However, as an additional consideration, Pelupessy & Portegies Zwart (2013) argue that the formation of planets inside evolved disks is rather unlikely. Timescale considerations seem to support this. Miller Bertolami (2016) claim that thermally pulsing AGB lifetimes are typically less than one Myr, while typical planet formation timescales are orders of magnitude larger (Alibert et al. 2005, 2013; Hori & Ikoma 2011). In addition, at present no mechanism, other than binarity, can explain the presence and formation of differentially rotating disks around evolved stars, besides binarity. Therefore, the planetary binary companion was probably present around the star before the onset of the AGB.

Assuming the companion was accompanying the star throughout its life, the disk could either be made up of mainly AGB wind material, or could contain a substantial portion of evaporated planet material. The evolution of the star to the giant phase would have substantially heated the planet, possibly even causing evaporation and outgassing. Using our determined ¹²CO over ¹³CO ratio as a proxy for the ¹²C/¹³C abundance fraction, we find the first scenario more probable. A ¹²C over ¹³C ratio of 10 (see Sect. 4.4) is typical of oxygen-rich AGB gas. In comparison with typical values in our Sun and solar neighbourhood (Asplund et al. 2009; Kobayashi et al. 2011) we would expect a ratio that is an order of magnitude higher were the disk made of substantial amounts of material from planetary origins. This also supports the findings in Paper I, where the spectral signatures at the companion location seem to indicate that it is accreting material via wind Roche lobe overflow.

6.5.2. Further evolution

Pelupessy & Portegies Zwart (2013) argue that the formation of planets inside evolved disks is unlikely. Only a very narrow combination of physical parameters permit planet formation. The physical properties of the disk of L₂ Pup lie well outside this range. Although no new planets are formed in the disk, the question remains whether the ones present can survive. The discovery of planets around a subdwarf-B star (Silvotti et al. 2014) shows that planets have survived the tumultuous giant phase. In addition, De Marco et al. (2015) deduces that some PN in the Kepler spacecraft field of view may also contain orbiting planets. Evidence thus strongly suggests that planets can persist through the final stages and influence the circumstellar morphology up to the white-dwarf stage.

Detailed studies of post-AGB morphologies Bujarrabal et al.2016 show a vast difference in length scale compared to the L₂ Pup system with disks ranging up to 1000 AU. In addition, their internal density structure is vastly different from the p = −3.1 slope found in this work (Sect. 6.3.2). Not many objects have been studied to a high degree of detail so it is difficult to gauge the significance of the difference between the evolutionary stages. We hypothesise nonetheless, and discuss a few scenarios which could explain the difference.

We speculate that the substantial difference in density distribution may be a consequence of the mass-loss rate of the central star. Low mass-loss rate AGB stars do not provide much material for the planet to interact with. This can be expected to result in a compact and low-density disk with a steep density slope. High mass-loss rates would pump more material into the interaction zone, causing wider and denser disks with flatter density slopes. However, as AGB stars evolve towards the post-AGB phase, it is believed that the mass-loss rate increases. If the circumstellar disk is not dense and massive enough to fully counter weigh this increased outward force, the disk is pushed outwards, slowing down its rotation substantially. This could explain the presence of the giant non-rotating tori observed around many post-AGB stars and planetary nebulae (e.g. Skinner et al. 1998; Sahai et al. 2006, 2011; Lagadec et al. 2011). Rotating disks that are massive enough to survive this “super wind” phase remain as such (e.g. Hillen 2015; Hillen et al. 2016).

The scale differences between AGB and post-AGB disks could also be explained by morphological constraints. As the star transitions to the AGB, equatorial wind material gets injected into a rotating disk by interaction with an orbiting planet. As more material accumulates, the disk becomes denser. At a certain point in time the density of the disk would reach a point where it would substantially obstruct the movement of the AGB wind along the equator. This would be the current stage of the L₂ Pup system, as there is a clear lack of radial gas flow inside the disk. The outflow material left unchanneled by the companion would therefore be forced towards the polar directions (along the path of least resistance), forming a wide bipolar outflow. Akashi & Soker (2008) suggest a mechanism via which the disk can accumulate mass expelled by this outflow, which would cause the disk to grow, and the density slope to flatten.

In any case, the study of the CSE of L₂ Pup strongly suggests that the presence of low mass (planetary) companions play a vital role in the shaping of the CSE of the AGB star, whose properties dictate the further morphological evolution of the system.

7. Summary

The circumstellar environment of L₂ Pup was imaged by ALMA at 15 × 18 mas resolution. The ¹²CO J = 3−2 and ¹³CO J = 3−2 emission show a rotating gas disk in a nearly edge-on orientation. Using the radiative transfer code LIME and the CASA post-processing tools, we constrained the physical properties of the disk:

• The radial temperature of the disk is characterised by three majorregimes. The temperature of the inner disk (<6 AU) is described by a steep drop from around 1000 K to around 500 K. Beyond this radius the temperature levels off at around 500 K, but drops again steeply to very cold temperatures (<100 K) at the outer edge of the disk (23 AU). Vertically the temperature of the inner disk is stratified, with cooler (1000 K) equatorial regions and rather hot inner disk walls (2500 K). We attribute this temperature stratification to complex hydrodynamical interactions at the outflow-disk boundary, which mixes dust formed at larger radii back into the innermost regions of the disk.

• The density of the disk is well described by the solution for a thin, vertically isothermal, non-self-gravitating disk in hydrostatic equilibrium. Its radial decay is expressed by a power law, with a power of −3.1. The flaring rate of the inner disk is quantified by a power of 0.2. Its inner rim is located at around 2 AU, which coincides with the location of the candidate binary companion detected in Paper I. This inner rim has a height of approximately 1.5 AU, indicating that it may be puffed up.

• The equatorial velocity was determined by Kervella et al. (2016), showing a sharp transition from Keplerian to sub-Keplerian at a radius of 6 AU. However, there is also an important vertical dependence, which reduces the expected Keplerian velocities of the inner disk gas to approximately 40% of its expected tangential velocity at high vertical offsets. The Gaussian turbulent velocity field component was deduced to be 0.5 km s^-1.

• The molecular abundance of ¹²CO was assumed to be 10^-4. The molecular abundance of ¹³CO was deduced to be 10^-5, resulting in a ¹²CO over ¹³CO abundance ratio of 10. This is in line with values of other oxygen-rich AGB stars in the literature. This suggests that the bulk of the gas in the disk originates from the star itself.

• The dust in the outer disk (r > 6 AU) seems to play a pivotal role in the thermodynamical and kinematical structure of the disk.

• Constraining the density profile permitted us to deduce the mass of the disk, which was found to be approximately 2 × 10^-4M_⊙. The angular momentum in the disk was found to be over 30 times the expected angular momentum contained within the star L₂ Pup. We calculate that a binary companion of ~1 Jupiter mass would be sufficient to explain the angular momentum contained in the disk.

Acknowledgments

This research made use of the SIMBAD and VIZIER databases (CDS, Strasbourg, France) and NASA’s Astrophysics Data System. We acknowledge financial support from the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, France. W.H. acknowledges support from the Fonds voor Wetenschappelijk Onderzoek Vlaanderen (FWO). L.D. acknowledges support from the ERC consolidator grant 646758 AEROSOL and the FWO Research Project grant G024112N. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.00141.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. KO acknowledges a grant from the Universidad Católica del Norte. We want to thank the referee for the rigorous review of this manuscript and detailed responses.

References

Appendix A: ¹²CO channel maps – ALMA data

Appendix B: ¹²CO channel maps – radiative transfer model

Appendix C: ¹³CO channel maps – ALMA data

Appendix D: ¹³CO channel maps – radiative transfer model

All Tables

All Figures