dc short=DC, long=data centre, dcn short=DCN, long=data centre network, jsd
short=JSD, long=Jensen-Shannon distance,

A Vectorised Packing Algorithm for Efficient Generation of Custom Traffic
Matrices

Christoper W. F. Parsonson,1,* Joshua L. Benjamin,1 and Georgios Zervas1

1 Optical Networks Group, Electronic & Electrical Engineering, UCL

<EMAIL_ADDRESS>
###### Abstract

We propose a new algorithm for generating custom network traffic matrices
which achieves $13\times$, $38\times$, and $70\times$ faster generation times
than prior work on networks with $64$, $256$, and $1024$ nodes respectively.

##  1 Introduction

dc have become critical tools for modern computational tasks. To meet the
ever-increasing demands of dc, recent years have seen a growth in the research
and development of next-generation dc optical systems [1]. However, most
researchers rely on simulations, which require the generation of synthetic
traffic. In doing so, they often make overly simplistic assumptions about the
characteristics of their generated traffic and develop systems which, in
practice, perform poorly under real-world conditions [2]. Furthermore, many
works omit open-accessing their synthetic traffic or even the methodology used
to generate it, bringing problems with reproducibility, benchmarking, and
cross-validation. The lack of a reproducible and high-fidelity synthetic
traffic generation tool has been a long-standing problem in the dc research
community.

Prior works [3, 4] have released traffic generators, but these were either
intended to be unrealistic, were for specific network topologies, required the
cumbersome use of inflexible configuration files, or lacked a reproducibility
guarantee. To address this, recent work presented TrafPy; an open source tool
for generating reproducible dc traffic with custom distributions and
characteristics [2]. However, the authors only demonstrated traffic generation
for $64$ network nodes; far smaller than the $O($$1000$$)$ node dc which are
becoming increasingly common place.

In this work, we first show that the original flow source-destination
assignment algorithm (‘packing’, see Section 2) used in the TrafPy generator
is a major bottleneck because its time complexity scales poorly with the
number of dc nodes $|N|$ for which $|F|$ flows are being generated. This
prevents the generation of traffic for large networks. Next, we propose a
novel vectorised packing algorithm which fits in with the rest of the authors’
traffic generation framework. Finally, we demonstrate our vectorised packer
achieving $13\times$, $38\times$, and $70\times$ faster generation times than
[2] on networks with $64$, $256$, and $1024$ nodes respectively with up to
$\approx 5$M traffic flows, with the speed-up factor increasing with the
network size. We expect this work to unlock a new realm of dc research at
scale and to further facilitate the development of next-generation systems and
common platforms for benchmarking networks. We note that while here we focus
on generating traffic for optical dc, the same traffic generation scheme and
vectorised packing algorithm could be re-purposed and applied to any network
system.

##  2 Custom Traffic Matrix Generation

Problem statement. dc traffic is made up of flows. A flow $f$ is fully
described by its size $f^{s}$ (how much information to send), arrival time
$f^{a}$ (when the flow requests to be transported through the dc, thus giving
rise to the inter-arrival time in a dynamic multi-flow setting), and source-
destination pair $f^{p}$ (which machines in the dc the flow is requesting to
be sent between). In the framework of [2], traffic generation is split into
two stages. In the first stage (‘shaping and sampling’), custom flow size and
inter-arrival time distributions are generated and sampled to attain a set of
sizes $\mathbb{b}^{s}$ and arrival times $\mathbb{b}^{a}$ for $f\in F$ flows
which match the target distributions within some jsd threshold111The jsd
$\in[0,1]$ is a measure of how similar two distributions are to one another
(lower is more similar), and the jsd ‘threshold’ as defined by [2] is a
constraint on how similar the generated traffic characteristics must be to the
target distributions.. In the second stage (‘packing’), given $\mathbb{b}^{s}$
and $\mathbb{b}^{a}$, the task is to assign each flow $f\in F$ to a source-
destination pair such that some target node distribution (a.k.a. traffic
matrix heat map) $\mathbb{P}^{N}$ with nodes $n\in N$ and corresponding
source-destination node pairs $p\in P$ is realised as closely as possible
without exceeding the load capacity limitations of any node. The authors
formulate this task by extracting the fraction of the overall load requested
by each pair $p\in P$ into an array, multiplying each element by the overall
dc’s target load rate to get the per-pair target load rate, and then again
multiplying each element by the simulation duration (the time between the
first and last flows’ arrivals) to get the total amount of information to load
onto each pair, $\mathbb{b}^{p,I}_{target}$, needed in order to achieve the
desired target node distribution $\mathbb{P}^{N}$. The packing task is
therefore reduced to finding the source-destination node pair assignments for
each flow $f\in F$ such that the difference between the actual and the target
per-pair total information loads,
$\mathbb{b}^{p,I}_{target}-\mathbb{b}^{p,I}_{actual}$, is $0$ or, where this
is not possible given any incompatibility between the target node distribution
$\mathbb{P}^{N}$ and the overall dc load rate, to match $\mathbb{P}^{N}$ as
closely as possible (see [2] for further details).

As shown in Figure 1a, as $|N|$ is increased, stage two (packing) becomes a
major bottleneck, taking $\approx$$1\,000\,000$ times longer than stage one
for $|N|=$$1024$. We therefore focus on optimising stage two.

Fig. 1: i) (a) The time for stages one (shaping and sampling) and two
(packing) when generating flows with the original packing algorithm. ii) The
packing (b) time and (c) Jensen-Shannon distance between the target and the
generated node distributions for the original and vectorised packing
algorithms when generating traffic for networks with different numbers of
nodes. (a) shows that the original packing algorithm is the major traffic
generation bottleneck of [2]. (b) shows that as the number of network nodes is
increased, the vectorised packer’s speed-up factor over the original algorithm
increases. (c) shows that both algorithms achieve the exact same resultant
node distribution. Note that the original algorithm’s time results for
$|N|=1024$ are extrapolations since it would have taken $\approx 200$ days to
run the packer.

Original packing algorithm. The original packing algorithm of [2] works by
sequentially iterating through the set of flows and, for each flow, conducting
two passes through the candidate source-destination pairs. In the first pass,
the packer attempts to match the target node distribution by looping through
all pairs, sorted in descending order of the total size of flow information
previously assigned, to find a pair which has not yet met its target
information load given the target node distribution and total flow arrival
duration provided. Failing to find such a node pair, the packer moves to the
second pass, whereby it again loops through each sorted pair but now in search
of a source-destination combination which, if allocated the flow in question,
would not exceed either the source’s or the destination’s maximum load
capacity given any prior flow allocations.

Algorithm 1 Vectorised packing algorithm pseudocode.

Input: $F$, $P$, $\mathbb{b}^{p,I}_{target}$

Output: $\mathbb{b}^{p,I}_{actual}$

Initialise: $\mathbb{b}^{p,I}_{actual}=0(|P|)$,
$\mathbb{b}^{p,c}=\frac{\text{node capacity}}{2}(|P|)$

for $f$ in $F$ do

$\mathbb{b}^{p,m}=$where$(\mathbb{b}^{p,c}-f^{s}<0,0,1)$ // Generate boolean
mask

$\mathbb{b}^{p,I,m}_{target}=\mathbb{b}^{p,I}_{target}[\mathbb{b}^{p,m}],\mathbb{b}^{p,I,m}_{actual}=\mathbb{b}^{p,I}_{actual}[\mathbb{b}^{p,m}]$
// Mask invalid pairs

$\mathbb{p}^{max}=$
argmax$(2\cdot\mathbb{b}^{p,I,m}_{target}-\mathbb{b}^{p.I,m}_{actual})$ // Get
furthest pairs

$p^{chosen}=$ random_choice$(\mathbb{p}^{max})$ // Randomly choose pair

update_trackers$(f,p^{chosen})$ // Assign flow to pair

end for

Vectorised packing algorithm. We negate the need for separate first and second
passes and for nested pair for loops by using vector array operations. We
begin by initialising the per-pair remaining capacity vector as the maximum
port capacity (half the per-node capacity, since it is split between the
source and destination ports) $\mathbb{b}^{p,c}$. We then sequentially iterate
through $f\in F$ and, for each flow $f$, we generate a boolean vector pairs
mask $\mathbb{b}^{p,m}$ which masks out any pair indices $i\in[0,...,|P|]$
which would exceed their load capacity were they to be allocated the flow in
question:

$\mathbb{b}_{i}^{p,m}=\begin{cases}0&\text{if }\mathbb{b}_{i}^{p,c}-f^{s}<0\\\
1,&\text{otherwise}\end{cases}$ (1)

We then apply this pairs mask to filter out any invalid pairs, thus ensuring
that any pair chosen from here on would meet the requirements of the second
pass of [2] and also reducing the time complexity of the argmax operation
below in Equation 2 (since the number of candidate pairs is now reduced).
Next, we take the masked candidate pairs’ current distances from the target
information loads, $\mathbb{b}^{p,I,m}_{target}-\mathbb{b}^{p,I,m}_{actual}$,
shift them by $\mathbb{b}^{p,I,m}_{target}$ in order to retain any skewness in
$\mathbb{P}^{N}$ for as long as possible given the overall dc load specified,
and find the pairs in this masked subset which are furthest from their target
information loads, $\mathbb{p}^{max}$:

$\mathbb{p}^{max}=\text{argmax}\bigg{(}2\cdot\mathbb{b}^{p,I,m}_{target}-\mathbb{b}^{p,I,m}_{actual}\bigg{)}$
(2)

In order to avoid any bias towards smaller pair indices and create the fade
phenomenon in the resultant traffic heat map [2], we randomly choose a pair
$p_{chosen}\in\mathbb{p}^{max}$ to which to allocate the flow $f$, thus
meeting the requirements of the first pass of [2]. Finally, we update the
current total information vector’s element for the chosen pair,
$\mathbb{b}^{p_{chosen},I}_{actual}$, and the remaining capacity vector
elements $\mathbb{b}^{p,c}$ for any pairs $p\in P$ which share either a source
or a destination with $p_{chosen}$. The pseudocode for this vectorised packer
is summarised in Algorithm 1.

##  3 Simulation Setup

Fig. 2: Custom traffic matrix distributions generated with $8$, $16$, $32$,
$64$, $128$, $256$, $512$, and $1024$ nodes, where the colour of each source-
destination pair corresponds to the fraction of the overall network load it
requests.

TrafPy enables the production of custom traffic distributions through the
control of a handful of parameters. These include the flow size and inter-
arrival time distribution parameters, the node distribution’s inter- vs.
intra-rack and skew node fractions, and the overall network load. To measure
the packing times for the original and vectorised packing algorithms, we
generated an assortment of custom traffic patterns typical for a ‘university’
dc222University dc service applications such as database backups, distributed
file system hosting, and multicast video streaming, with $\approx 70\%$ of
traffic being inter-rack and $\approx 20\%$ of nodes requesting $\approx 55\%$
of the traffic load. as detailed by [2] for networks with $4$ racks and
$|N|=\\{8,16,32,64,128,256,512,1024\\}$ nodes (see Figure 2). We assumed an
optical dc network with with an overall network load rate of $50\%$. For each
traffic matrix, we generated $|F|=5\cdot|N|^{2}$ flows to ensure non-sparse
packing. Each packing algorithm was ran on a shared cluster with an Intel Xeon
ES-2660 CPU across $4$ seeds to ensure reliable packing times given the
variance in use of the shared cluster, with the 95% confidence interval bands
plotted for any metrics recorded.

##  4 Results & Discussion

Figure 1b shows the packing times taken by the original and the vectorised
packers when generating the distributions shown in Figure 2. The vectorised
packer achieved a $\approx 38\times$ speed-up over the original packer on the
$|N|=264$ traffic matrix and $\approx 70\times$ on the $|N|=$$1024$ matrix.
Although the vectorised algorithm was slightly slower than the original packer
on the smallest $|N|=8$ network due to performing a where vector operation on
all pairs, the absolute generation time was still $O(s)$ and this additional
overhead quickly becomes negligible across $|N|>8$ networks.

To verify that our proposed vectorised packing algorithm was generating the
same node distributions as the original packer used by [2], we measured the
jsd between the target and the generated node distributions for each algorithm
(see Figure 1c). As expected, both packers deterministically reach the same
solution, but the Jensen-Shannon distance will not be exactly $0$ for either
due to the incompatability between the $50\%$ network load and the skewed
target distribution (see [2] for more details).

In conclusion, we have proposed a flow source-destination pair assignment
algorithm which makes novel use of vector array operations to achieve orders
of magnitude faster traffic generation times than the original algorithm used
by [2] when generating custom traffic matrices. This work significantly
improves the utility of an open source traffic generation framework in order
to aid the production of high-fidelity traffic patterns and to test and
develop network systems at scale.

Acknowledgements: EPSRC Distributed Quantum Computing and Applications
EP/W032643/1; the Innovate UK Project on Quantum Data Centres and the Future
10004793; OptoCloud EP/T026081/1; TRANSNET EP/R035342/1; the Engineering and
Physical Sciences Research Council EP/R041792/1 and EP/L015455/1; the Alan
Turing Institute; and Horizon Europe Dynamos.

## References

  * [1] M. Khani et al., “Sip-ml: High-bandwidth optical network interconnects for machine learning training,” SIGCOMM (2021).
  * [2] C. Parsonson et al., “Traffic generation for benchmarking data centre networks,” Optical Switching Networking (2022).
  * [3] M. Alizadeh et al., “pfabric: Minimal near-optimal datacenter transport,” SIGCOMM (2013).
  * [4] W. Bai et al., “Enabling ecn in multi-service multi-queue data centers,” NSDI (2016).


Vol.0 (20xx) No.0, 000–000

11institutetext: Netherlands Institute for Radio Astronomy ASTRON, 79901 PD
Dwingeloo, the Netherlands<EMAIL_ADDRESS>
22institutetext: Center for General Education, Institute for Comprehensive
Education, Kagoshima University, 1-21-30 Korimoto, Kagoshima 890-0065, Japan
33institutetext: Xinjiang Astronomical Observatory, 150 Science 1-Street,
Urumqi, Xinjiang 830011, P.R. China
44institutetext: Amanogawa Galaxy Astronomy Research Center, Graduate School
of Science and Engineering, Kagoshima University, 1-21-35 Korimoto, Kagoshima
890-0065, Japan
55institutetext: School of Natural Sciences, University of Tasmania, Private
Bag 37, Hobart, Tasmania 7001, Australia
Received 20xx month day; accepted 20xx month day

# Fallback in bipolar Planetary Nebulae ?

Willem A. Baan 112233 Hiroshi Imai 2244 Gabor Orosz 223355

###### Abstract

The subclass of bipolar Planetary Nebulae (PN) exhibits well-defined low-power
outflows and some show shock-related equatorial spiderweb structures and
hourglass structures surrounding these outflows. These structures are
distinctly different from the phenomena associated with spherical and
elliptical PN and suggest a non-standard way to simultaneously energise both
kind of structures. This paper presents evidence from the published literature
on bipolar PN Hb 12 and other sources in support of an alternative scenario
for energising these structures by means of accretion from material shells
deposited during earlier post-AGB and pre-PN evolutionary stages. In addition
to energising the bipolar outflow, a sub-Eddington accretion scenario could
hydrodynamically explain the spiderweb and outer hourglass structures as
oblique shockwaves for guiding the accreting material into the equatorial
region of the source. Estimates of the accretion rate resulting from fallback-
related spherical accretion could indeed help to drive a low-power outflow and
contribute to the total luminosity of these sources.

###### keywords:

ISM: planetary nebulae: general – ISM: planetary nebulae: individual Hb 12,
Hen 2-104, MyCn 18 – ISM: jets and outflows – Stars: evolution

##  1 Introduction

Planetary Nebulae (PNe) represent the post-evolution stages of medium mass
stars leaving the Asymptotic Giant Branch (AGB) before its stellar core
eventually turns into a White Dwarf (WD) (Balick & Frank 2002). At the end of
AGB evolution a superwind depletes the AGB envelope until the stellar
structure changes, the photospheric radius shrinks, and the stellar effective
temperature rises (de Marco 2009). The further evolution of PNe is then
affected by the occurrence of multiple ‘late and very-late thermal pulse
cycles’ at the stellar core (Balick & Frank 2002), the fallback of AGB shells
(Perinotto et al. 2004; Chen et al. 2016), and the ionisation of surrounding
material by the central star, which make the PN move back and forth across the
Hertzsprung-Russell (HR) diagram (Werner & Herwig 2006). Radiation
hydrodynamical modelling of AGB-to-PN evolution shows how the central-star
mass and the (post-)AGB mass-loss history determine the shape and the
kinematics of the resulting planetary nebula (Perinotto et al. 2004;
Schonberner et al. 2005).

The complexity of the evolutionary processes of PNe can provide for the wide
variety in colour and shape of spherical and elliptical PNe. However, the
existence of spherical shells and ionization fronts does not explain any of
the structures observed in bipolar PN (Balick & Frank 2002). The small number
of bipolar (elongated) PN structures indicate the existence of a directed
outflow that has been described as a stellar wind emanating from the inner
layers of the stellar structure and blowing out along the symmetry axis
(Mellema & Frank 1995; Balick & Frank 2002). A few sources exhibit a well-
defined ‘hourglass’ axi-symmetric structure that defines the interaction
boundaries of an expanding outflow with the surrounding circumstellar material
(Sahai et al. 1999). Alternatively, bipolar sources have been associated with
a companion star in order to create a system asymmetry or have been explained
by accretion in a symbiotic binary system (SySt), as in the case of the Twin
Jet Nebula (Corradi et al. 2011). One out of five PNe have indeed been
confirmed to host a binary companion (Miszalski et al. 2015) but the evidence
for a ‘binary hypothesis’ is not yet very strong as compared with the ‘single
star hypothesis’ for PNe (de Marco 2009).

Two of the prominent hourglass sources, Hb 12 and MyCn 18, show a Spiderweb
arc structure in the equatorial plane of the system (perpendicular to the
hourglass) that has (following tradition) been explained by repeated stellar
wind outflows (Clark et al. 2014; Sahai et al. 1999), by sequential ejections
resulting from a thermal nuclear runaway (Kwok & Hsia 2007), or by episodes of
activity during earlier evolutionary stages (Vaytet et al. 2009). However,
there is no clear explanation for a simultaneous outflow along the polar axis
and the structure in the equatorial plane. Similarly, intermittent jets or
winds blown by an accreting binary system have also been invoked to explain
both the inner arc structure of the Red Rectangle (HD 44179) and its
(perpendicular) outer structure (Soker 2005; Bujarrabal et al. 2016). Since
the bipolar PNe do not show any ionised shells that are characteristic of
evolved PN, they seem to be either at a very early evolutionary stage or at
the end of PN evolution. The suggestion has indeed been made for these to be
low-excitation and young sources that only recently started ionising their
surrounding shells and represent a short transition stage toward full
development (Clark et al. 2014).

The evolution of PN is known to be strongly affected by fallback by material
deposited in the surrounding interstellar environment by AGB and PN
evolutionary activity. However, no clear scenario has been established to
recognise when fallback actually occurs. While fallback may happen unnoticed
observationally when the accretion rate is very low and streams of fallback
material may leak into the atmosphere of the central star of the PN, a
definite accretion signature should be established when the accretion rate
become high. When the accretion rate becomes higher because of the
availability of accretion material, sub-Eddington accretion may happen where
the low luminosity of the PN source cannot stagnate the accretion flow. The
power generated during this (possibly short-lived) accretion could still be
sufficient to energise activity in the PN. Although existing explanations for
bipolar PN sources have traditionally relied on nuclear activity of the source
itself, an accretion scenario should also be considered to drive an outflow
along the polar axis and also explain any structures close to the equatorial
plane.

In this paper, we seek qualitative evidence from the published literature that
would support the possibility of semi-spherical sub-Eddington accretion unto a
bipolar source such as PN Hb 12. The origin of the accreting material may be
the ejecta during earlier AGB stages that have been deposited in the
surrounding medium of a post-AGB environment of a bipolar PN. Alternatively,
accretion could also occur during or at the end of the PN evolutionary path
with diminished/ceased stellar activity but the accretion rate from the ISM
may not be sufficient during these stages. The evidence for an accretion
scenario in a bipolar source such as PN Hb 12 will set the stage for further
research of bipolar PN. After describing the structural evidence of bipolar PN
Hb 12 in Section 2, available observational evidence is presented in Section 3
that would directly support or be consistent with an accretion-outflow
scenario.

Figure 1: (Left) A composite optical/infrared HST image shows the inner
structure of the Matryoshka Nebula Hubble Hb 12 with the [N II]($\lambda$6584
$\mu$m) emission and the NICMOS F160W image showing shocked [Fe
II]($\lambda$1.6435 $\mu$m) emission in the boundary of the bipolar Outflow.
Photo: ESA/NASA/J. Schmidt. (Right) HST NICMOS F160W image from the HST Legacy
Archive showing the detailed [Fe II] emission structures of the Spiderweb. The
arrows point to the Spiderweb arcs detected in NICMOS F220W H2
($\lambda$2.1214 $\mu$m) data. The equatorial arcs are integrated with an
outer and larger bipolar Hourglass structure. The field of view is
40$\arcsec$x40$\arcsec$. Image from (Clark et al. 2014). Figure 2: Schematic
of the structural components of PN Hb 12. The red components relate to the
outflow patterns and the blue components relate to inflow patterns. The Inner
Hourglass outflow boundary, the Outer Hourglass and Spiderweb oblique shocks
are clearly seen in the images of Hb 12 in Figure 1. The black arrows in the
upper right quadrant designate the accretion flow pattern through these (blue)
oblique shock structures. A stagnation point of the (red) outflow in Hb 12 has
also been observed at a significant distance of 75$\arcsec$ from the star.

##  2 The bipolar structure of PN Hb 12

The bipolar Matryoshka Nebula PN Hb 12 has been very well studied and provides
specific fundamental insights into the workings of a bipolar PN (Hora & Latter
1996; Welch et al. 1999; Kwok & Hsia 2007; Vaytet et al. 2009; Clark et al.
2014). Figure 1(left) presents a composite optical/infrared Hubble image of
the central structure of Hb 12, which shows a hourglass-shaped ‘bipolar
outflow’ channels seen in narrow band [N II] ($\lambda$ 6584) observations
(Balick 2003; Kwok & Hsia 2007). The axis of the hourglass outflow structure
in PN Hb 12 with concentric ring structures in its boundary is tipped
clockwise by 10$\deg$ and with the upper part tilted away from the plane of
the sky by 25$\deg$ (Vaytet et al. 2009). In addition, prominent infrared [Fe
II] ($\lambda$1.6435 $\mu$m) line emission shows the two edge section of an
equatorial and axi-symmetric Spiderweb structure together with the inner
sections of a larger scale ‘Outer Hourglass’ structure (Clark et al. 2014).
Figure 1 (right) again shows the detailed [Fe II] equatorial Spiderweb
structure surrounding the star and the continuation towards a bipolar Outer
Hourglass structure that is less collimated than the (inner) outflow hourglass
structure. This Outer Hourglass (shock) structure is expected become less
prominent with distance and will close beyond the stagnation point of the
bipolar outflows; this closure has not yet been observed in Hb 12. The nearly
face-on Hb 12 displays a very detailed east – west symmetry of the axi-
symmetric Spiderweb but also a north – south symmetry at either side. These
structural details of the Spiderweb will be essential for understanding what
is happening in these regions.

The bipolar outflow, the Spiderweb structure, and the outer hourglass
structures found in Hb 12 and other bipolar sources are the defining
structural components to be explained consistently in modelling these sources.
Prominent co-axial structures emanating from the core region are also found in
the Hourglass Nebula PN MyCn 18 (Bryce et al. 1997), as well as in symbiotic
bipolar nebulae such as the Southern Crab Hen 2-104, Menzel’s young PN Ant
Nebula Mz 3, Minkowski’s Twin Jet (or Butterfly) Nebula M 2-9 (Corradi et al.
2001; Clyne et al. 2015, and references therein), and the post-AGB pre-PN
(Rotten) Egg Nebula OH 231.8+04.2 (see Balick et al. 2020).

In this paper we adopt the distance to Hb 12 to be 2.24 kpc, which renders the
15$\arcsec$ extent of these structures to be 0.16 pc (34000 AU). In reality
the distances to Hb 12 (estimates from 2.24 to 14.25 kpc) and other PNe are
highly uncertain, which makes size determinations difficult (Clark et al.
2014). Considering the extreme detail observed in the Spiderweb and the
bipolar outflow in PN Hb 12, the distance may actually be less than the 2.24
kpc estimate. Further hydrodynamic modelling of the Spiderweb structures may
facilitate a more accurate distance determination.

##  3 Evidence for outflow and inflow

The presence of outflows in bipolar PN has been long recognised in the
literature but the simultaneous occurrence of a bipolar outflow combined with
Spiderweb and Outer Hourglass structures have not led to satisfactory
explanations. In the following sections, we consider the existing literature
for any evidence that accretion inflow may be driving the outflows in bipolar
sources and that an accretion scenario consistently explains the nature of all
observed structural components in Hb 12.

###  3.1 A Semi-Spherical Accretion and Outflow Scenario

The outflow structure, the Spiderweb, and the outer structures observed in Hb
12 are delineated by forbidden [N II] ($\lambda$6584) and infrared [Fe II]
($\lambda$1.6435 $\mu$m) line emissions that all may be assumed to be related
to shock heating and excitation. While the line emission associated with the
outflow may be easily understood, the other line emission structures would
suggest a non-outflow explanation. The existence of an accretion inflow close
to the stellar surface would naturally result in an axisymmetric system of
oblique shocks that guides the initially radial accretion inflow into the
equatorial region and around the outflow structure as captured in the diagram
of Figure 2). Hydrodynamically an oblique shock stands under an angle with the
flow direction and will thermalise part of the flow velocity component
perpendicular to the shock, while leaving the parallel velocity component
intact. Assuming that the underlying star has ceased its own shell burning
activity and has no substantial stellar wind, the oblique shocks of the Outer
Hourglass serve to decelerate and redirect the in-falling gas at large
distances in order to guide the flow around the bipolar outflow regions and
into the axisymmetric oblique shocks of the Spiderweb structure in the
equatorial region. The very structured shock patterns observed around PN Hb 12
and other sources may thus be consistent with a hydrodynamic inflow pattern
resulting from accretion, either spherical or from a distant companion star.
While full hydrodynamical modelling of this inflow-outflow structure will need
to further confirm this scenario, it will also provide a measure of the
accretion rate, the strength of the oblique shocks, and the actual scale size
of the shock patterns (and indirectly the distance).

###  3.2 Bipolar Outflow in the Inner Hourglass

The bipolar Inner Hourglass outflow structure in Hb 12 represents a prominent
feature in the source with the base seen in forbidden [Fe II] line emission
and the outline of the whole structure seen in the lower density tracing
forbidden [N II] emission (Fig. 1). The observed structure clearly suggests a
low power outflow expanding into a pressure gradient and energising a widening
boundary layer as outlined by the [Fe II] and [N II] line emission (Kwok &
Hsia 2007; Vaytet et al. 2009). The shape of the boundary of this hourglass
outflow channel suggests consistency with a supersonic and accelerating fluid
outflow from a simple DeLaval-like nozzle expanding into a circumstellar
medium with a radial pressure profile p(r) $\propto$ r-2 (using the formalism
from Landau & Lifshits 1959; Choudhuri 1998). The boundaries of this outflow
are self-adjusting and are determined by the equilibrium between the internal
pressure of the accelerating flow and the radially varying ambient pressure in
the circumstellar envelope. The supersonic outflow along the polar axis starts
at a narrowed throat formed in the atmosphere where the gas reaches a proton
escape velocity equaling the local velocity of sound. If the 90 km s-1 base
velocity is indicative of the local escape velocity from a 1 M⊙ star, then
this Mach nozzle may form at a radius of about 1012 cm and require a local
temperature of 105 K, which is higher than the effective surface temperature
of 35,000 K (Hyung & Aller 1996). Our first order modelling of the case with a
quadratic pressure gradient shows that the flow accelerates early in the
nozzle and may reach a final velocity of at least 2.5 times the sound velocity
at the throat of the nozzle. Sustained operation of the nozzle outflow would
require sustained replenishment of heated atmospheric material. Detailed
modelling may provide an estimate of the power and the scalesize of the
outflow.

The boundary layer observed in Hb 12 and other sources will naturally develop
around the outflow where the edges of the flow pattern are slowly decelerating
because of energy and momentum being dissipated. This boundary layer of Hb 12
is very prominent with shock-tracer [Fe II] emission lines at the onset of the
outflow and is seen in [N II] emission away from the star as the interaction
becomes less intense (Clark et al. 2014). The velocities at the base of the Hb
12 outflow are estimated to be 90 $\pm$15 km s-1 in the redshifted (north) and
blueshifted (south) asenses and consistent with the velocity estimate
mentioned above. Similar bipolar structures and outflow velocities have been
detected in the PNe MyCn 18 (V = $\pm$48 km s-1) (Corradi & Schwarz 1993), in
Mz 3 (V = 130 km s-1) (Santander-García et al. 2004), and the Inner Hourglass
structure of the symbiotic Mira Hen 2-104 (V = 50 km s-1) (Corradi et al.
2001; Clyne et al. 2015).

At a certain distance from the central star, bipolar outflows are expected to
stagnate and produce lobes containing accumulated outflow and entrained
material that are analogous to the radio lobes of extragalactic sources. The
stagnation lobes in the low-power outflow of PN Hb 12 appear as multiple knots
at a distance of 75$\arcsec$ (0.84 pc if at 2.24 kpc) that are emitting low-
ionisation [N II] line emission (Vaytet et al. 2009). The velocities of these
knots are at $+$275 (north) and $-$296 km s-1, which is slightly higher than a
few times the base velocity 90 km s-1. Stagnation points have also been
identified in the sources Hen 2-104 and the Rotten Egg Nebula (Corradi et al.
2001; Kwok & Hsia 2007).

The multiple coaxial rings in the walls of the hourglass outflow structures of
Hb 12 are also seen as low-ionisation [N II] thermal emission (Kwok & Hsia
2007) and would result from energy dissipation in the boundary layer (see Fig.
1). These rings appear roughly equally spaced at approximately 850 AU (if at
2.24 kpc) and may result from some axi-symmetric Kelvin-Helmholtz-like
instability resulting from the velocity gradients across the boundary layer.
Alternatively, these rings could represent decelerating shocks only in the
boundary layer as a result of the decreasing ambient pressure shaping the
outflow. Similar rings are also prominent in the sources MyCn 18 (Sahai et al.
1999), M 2-9 (Kwok & Hsia 2007), and Hen 2-104 (Corradi et al. 2001).

###  3.3 Evidence for inflow at the Outer Hourglass

The Outer Hourglass oblique shock structure appears as two axisymmetric
elliptical surfaces that reach around the star and interconnects with Inner
Spiderweb at its outer edges (see Fig. 2). In PN Hb 12 only the base of the
shock-related [Fe II] Outer Hourglass emission structures can be seen in (Fig.
1) but they are expected to extend towards the stagnation points in the
outflows (Vaytet et al. 2009). The axisymmetric region between the boundary
layer of the Inner Hourglass outflow and the Outer Hourglass shocks would also
be filled with in-falling material that has been shock-heated and may show H2
and [N II] tracer emission. All shock-related [N II] and H2 emissions within
the Outer Hourglass structure are expected to decrease in intensity away from
the star because of decreasing densities and decreasing perpendicular
velocities in the oblique shocks. While Hb 12 does not (yet) exhibit such
emissions, the filamentary [N II] $\lambda$6583 emitting cloudlets in PN Hen
2-104 trace a very similar Outer Hourglass structure far beyond the Inner
Hourglass outflow (Corradi et al. 2001; Corradi 2003; Kwok & Hsia 2007),
Similar emissions are found in the Twin Jet Nebula M 2-9, and PN Mz 3 (see
Kwok & Hsia 2007; Clyne et al. 2015).

A long-slit velocity-position diagram taken along the major axis of the Hen
2-104 system shows that the bipolar outflow has (accelerated) positive
velocities in the north up to projected $+$90 km s-1, and negative velocities
in the south (Corradi & Schwarz 1993; Corradi et al. 2001). However, the
diagram also shows opposite velocities in the north as well as in the south
that increase up to projected of 30 km s-1at the edge of the Outer Hourglass
structure. While the bipolar Inner Hourglass shows an outflow pattern, the
filaments inside the Outer Hourglass of Hen 2-104 would indeed represent the
expected inflowing accretion component.

###  3.4 Inflow through the Spiderweb Structures

The apparent shock structures in both the equatorial Spiderweb and the Outer
Hourglass structures are consistent with a large scale hydrodynamic flow
pattern leading towards the stellar core in Hb 12. The series of oblique
shocks redirect and decelerate the accretion flow, raise the density and the
temperature, and funnel the flow into the equatorial region. The shock nature
of these structures may be supported by the detection of the shock tracer [Fe
II] ($\lambda$1.6436 $\mu$m) emission, continuum ($\lambda$1.6$\mu$m), and the
H2 ($\lambda$2.1214 $\mu$m) line emission (Hora & Latter 1996; Clark et al.
2014). The two tangential edge sections of this axisymmetric system of
standing oblique shocks show in great detail how the inflow is being guided
onto an equatorial surface of the star (Figure 2). The nearly face-on position
of Hb 12 also shows the north-south symmetry in the Spiderweb, which indicates
that the flow pattern near the star is the same above and below the equatorial
plane. Incidentally, these oblique shocks work on the same principles as a
diffuser (air intake) of a supersonic airplane. The inner Spiderweb and the
Outer Hourglass form an interconnected oblique shock network that takes the
initially spherical accretion flow, bypasses the region of the bipolar
outflow, and brings it into the equatorial region of the star.

Analysis of the molecular hydrogen emission in Hb 12 suggests that collisional
processes would dominate over UV absorption for exciting the H2 emission when
the kinetic temperature Tk $>$ 1000 K (Hora & Latter 1996). The temperature of
free-falling material at a distance R from a star would be considerably high
at Tff = 1.08 x 105 M⊙ R${}_{AU}^{-1}$ K. This suggests that shock heating of
the accreting material would completely dominate the excitation of H2 inside a
radius of about 100 AU from the star.

In addition to PN Hb 12, evidence for similar Spiderweb structures has been
found in bipolar PNe such as MyCn 18 (Bryce et al. 1997), Henize 2-104
(Corradi et al. 2001; Clyne et al. 2015), R Aqr (Corradi & Schwarz 1995), and
Mz 3 and M 2-9 (Clyne et al. 2015).

###  3.5 Infall in the vicinity of the stellar surface

A final standing shock must exist above the equatorial stellar surface in
order to further decelerate the (then) nearly radial flow pattern to a
subsonic velocity and to let this inflow integrate with the stellar atmosphere
and a possible circumstellar disk. Existing observational data of the H2 and
[Fe II] emissions in Hb 12 taken just north of the core, at the core, and just
south of the core provide the evidence for this inflow inside the Spiderweb of
Hb 12 (Figure 2 in Clark et al. 2014). At the ‘northern edge’ of the nuclear
region and covering the edge of the Inner Hourglass, the shock tracer [Fe II]
emission shows a redshifted northern outflow with a (projected) velocity in
the range of $+$30 to $+$60 km s-1. However, there are also positive
velocities in the [Fe II] and H2 data south of the core at the edge of the
Spiderweb in the range 0 to $+$30 km s-1. Similarly, the blueshifted
‘southern’ outflow is seen south of the core in the [Fe II] data at $-$60 to
$-$30 km s-1, while the H2 data also shows emission north of the core at the
Spiderweb structure in the range $-$30 to 0 km s-1. While the strongly
positive and negative velocity components in the [Fe II] data can be explained
as being part of the bipolar outflow, the smaller velocity components in the
same direction at the opposite sides of the core do not fit the picture. A
simple interpretation of these components would be that velocities in the same
direction (positive or negative) at either (north and south) side of the
tilted source indicate that one can represent an outflow and then the other is
an inflow. This suggests that the observational data presented by Clark et al.
(2014) show a outflow-inflow velocity pattern that is consistent with an
accretion scenario.

###  3.6 Evidence for a disk around the central star

The accretion flow into an equatorial footprint may result in the formation of
a thick disk/torus around the central star as a form of temporary mass storage
before material transfers into the atmosphere of the star. Evidence for such
an equatorial torus has indeed been found in Hb 12 using the bright He I
(2.0585 $\mu$m) and Br$\gamma$ (2.1655 $\mu$m) recombination line emission
showing a central cavity with a 0.2$\arcsec$ x 0.41$\arcsec$ (450 $\times$ 930
AU if at 2.24 kpc) inner torus structure (Clark et al. 2014). However, the
observed velocity profile perpendicular to the symmetry axis indicates a
rotation velocity of 90 km s-1, which would suggest an orbital radius in the
disk of only 0.22 AU around a one M⊙ star (and an uncomfortable distance of
only 2.2 pc to the observer).

The presence of a disk in a bipolar PNe may also be confirmed by any maser
activity to be detected at their central regions. OH (1665 MHz) and H2O (22.12
GHz) maser features at the core region of the bipolar PN K 3-35 suggest a disk
with an extent of 98 AU (if at 5 kpc) (Miranda et al. 2001). In addition, the
OH and H2O features at the core K 3-35 are redshifted relative to the OH (1667
MHz) maser components found in the bipolar outflows, which would further
confirm inflow into the central region. OH and H2O masers are also found in
the core regions of other bipolar PNe (Gómez et al. 2008; Uscanga et al. 2014;
Gómez et al. 2015). Similarly, the in-falling SiO (43 GHz) masers in the
bipolar post-AGB pre-PN (Rotten Egg) OH 231.8+04.2 have been interpreted as
being part of a torus (Sánchez-Contreras et al. 2002).

The shocked-heated material in the disc around the stellar core will
contribute to emission of the source in addition to ongoing or enhanced shell
burning possibly triggered by the accreted material. The estimated effective
surface temperature of Hb 12 is 35,000 K based on modelling of the observed
emission lines (Hyung & Aller 1996), which suggests an early O-star
classification and would indicate an early stage in PNe evolutionary models
(Perinotto et al. 2004; Werner & Herwig 2006).

###  3.7 The possibility of sub-Eddington Accretion

Considering that bipolar PNe do not (yet) exhibit the shell-like
characteristics of spherical and elliptical PN, the possibility should be
considered that the bipolar activity is driven by accretion in sources such as
Hb 12. For early-type PN the material may be accreted from a circumstellar
cocoon/envelope or a large circumbinary disk that was formed during the later
stages of post-AGB evolution (Perinotto et al. 2004; Herwig 2005; Akashi &
Soker 2008). Simulations show that depending on the mass of the underlying
star, such envelopes/shells may exist at distances ranging from a few to 10
$\times$ 1017 cm and have particle densities ranging from 20 to a few hundred
cm-3 (Perinotto et al. 2004; Chen et al. 2016). During the early evolutionary
PN phases lasting a few 104 yr, the free-fall timescale
$t_{ff}=R^{3/2}/(2GM)^{-1/2}$ suggests that material within a radius of $R$ =
0.22 $\times$ 1017 cm can fall back unto the star.

An accretion model that most adequately describes sub-Eddington accretion from
a (distant) surrounding medium is the Bondi-Hoyle spherical accretion model
(Bondi 1952; Mestel 1954). While the luminosity of the source is not
sufficient to halt the spherical accretion, a hemispherical balance may be
established between the polar outflows and equatorial accretion flow regions,
while the accreting material away from the star does not yet know about any
complicated flow pattern close to the star. Assuming a uniform surrounding
medium, a spherical accretion rate may be expressed for the temperature
limited case with $\gamma$ = 3/2 as (Bondi 1952):

${\dot{M}_{acc}}=2\pi\alpha(GM_{st}^{2}(V^{2}+c_{s}^{2})^{-3/2}\rho_{ism})=3.5\times
10^{-9}M_{st}^{2}\,n_{100}\,T^{-3/2}_{10}\,\hbox{M${}_{\odot}$}/yr,$

where $M_{st}$ is the stellar mass assumed to be 1 M⊙, $V$ is the spatial
velocity of the star and is assumed to be zero, $c_{s}$ is the speed of sound
that depends on the temperature of the surrounding medium expressed in units
of $T_{10}$ = 10K, and $\rho_{ism}$ is the particle mass density of the medium
far from the star expressed in units of $n_{100}$ = 100 cm-3. The factor
$\alpha$ is between 1 and 2 depending on environmental conditions and a value
of 1.5 is adopted here. The accretion rate favours higher stellar mass, higher
gas densities, and lower temperatures. The luminosity resulting from the
estimated accretion rate may be estimated as:

$L_{acc}=\dot{M}_{acc}GM_{st}/R_{st}=2.9\times
10^{33}M_{1}^{3}\,n_{100}\,T^{-3/2}_{10}\,R_{10}^{-1}\,erg/s,$

where the stellar radius $R_{st}$ is expressed in units of 1010 cm.

An estimate of the luminosity of Hb 12 may follow from the estimated effective
surface temperature of 35,000 K based on modelling of the observed emission
lines (Hyung & Aller 1996), which results in a (minimum) blackbody luminosity
for a stellar core with an effective radius $10^{11}$ cm to be
$L_{bb}=8.8\times 10^{36}$ erg/sec. This value is less than the Eddington
Luminosity LEdd = 1.2 $\times 10^{38}M_{1}$ erg/sec, which would verify that
sub-Eddington accretion is possible for PN Hb 12. Assuming some reasonable
parameters for the surrounding environment, a scenario with spherical
accretion from fallback of circumstellar material could account for part of
the apparent luminosity of the star and result in driving the outflow activity
of Hb 12.

##  4 Conclusions

The presence of bipolar outflows seen in a subset of the PNe population may
naturally result from accretion unto the star and the existing observational
evidence of bipolar PN Hb 12 and similar sources would (already) support such
a scenario. Published literature provides the evidence that fallback of
gaseous shells deposited during the earlier stages of post-AGB and pre-PN
evolution may indeed provide sufficient material for sub-Eddington accretion.
The shape of the boundaries of the low-power outflow in Hb 12 suggests that it
is an accelerating and supersonic flow into an environment with a quadratic
pressure gradient. A DeLaval-like nozzle may form in the stellar atmosphere
where the gas reaches escape velocity equalling the sound speed from where the
flow expands into the pressure gradient. The supplement of the accretion
energy to the stellar atmosphere would be a trigger for establishing a bipolar
outflow. The low luminosity of the stellar PN core is unable halt the
accretion flow and establish an omni-directional stellar wind.

The unusual Spiderweb and Outer Hourglass structures observed in PN Hb 12 (and
in similar sources) may be interpreted hydrodynamically as oblique shocks that
decelerate the initially spherically accretion flow and redirect the flow
pattern around the hemispherical sections of the bipolar outflows and into the
equatorial region of the star. Evidence for infall follows from a detailed
inspection of the imaging data of the [Fe II] and H2 line emission of HB 12,
where inflow velocity components are indeed co-located with outflow velocities
in the core region. Similarly, the axi-symmetric region between the Inner
Hourglass boundary layer of the bipolar outflow and the oblique shocks of the
Outer Hourglass would be filled with inflowing gas, as seen in emission in the
similar Spiderweb source Hen 2-104. Also the presence of an accretion disk
surrounding Hb 12 with a clear velocity gradient confirms that infalling
material is collecting around the star. The estimate of the mass accretion
rate using the Bondi-Hoyle model for spherical accretion shows that for
reasonable parameters for the surrounding ISM environment, accretion of
fallback shells could add significantly to the luminosity of the pre-PN source
and provide enough energy to drive the bipolar outflows.

An episode with sub-Eddington accretion undoubtedly defines a separate
evolutionary stage in the life of a PN. If the accreting material results from
post-AGB activity, the star will acquire new material in its atmosphere,
although some will be bipolar-ejected again and may fall back at a later time.
After this accretion stage is finished, the star will initiate a further path
of PN evolution with a new atmosphere. If this type of accretion happens
during a later evolutionary stage of the PN, this will mark a temporary
interruption of the PN activity to be followed by renewed PN activity with a
replenished atmosphere. Since the availability of material to be accreted is
an important factor, episodes of accretion depend very much on the details of
AGB to PN and it is likely that stages of accretion will occupy only part of
the evolutionary lifetime of a PN.

Further verification of the hydrodynamic interpretation of the accretion
scenario will follow from (ongoing) modelling of the spherical Bondi-Hoyle
accretion flow and the oblique shocks of the Spiderweb and the Outer Hourglass
structures. Hydrodynamic modelling of the accretion flow and the supersonic
DeLaval-nozzle bipolar outflow will provide an indication of the actual sizes
of the oblique shock components and the outflow, and provide an independent
estimate for the (uncertain) distance of these sources. Modelling efforts of
the bipolar outflow and of the accretion flow have been initiated. A general
verification of the accretion scenario would also result from (planned) radio
observations of the shocks of the Spiderweb and the Outer Hourglass structure
in Hb 12 and the (potential) detection of maser action in the shocked dense
material of the disk in the system.

## Acknowledgements

WAB thanks the Graduate School of Science and Engineering of Kagoshima
University and the Department of English Communication of Toyo University for
their hospitality during the JSPS Invitation Program for Foreign Researchers
(S15126). WAB also thanks the XinJiang Astronomical Observatory for their
hospitality during the Program of High-end Foreign Experts of the XinJiang
Uygur Autonomous Region (grants 20166500004 and 20176500001) and under the
PIFI grant No. 2019VMA0040 from the Chinese Academy of Sciences. HI was
supported by the JSPS KAKENHI program JP16H02167. GO acknowledges support from
the ARC Discovery project DP180101061 of the Australian Government, the CAS
LCWR Program 2018-XBQNXZ-B-021 of China, the Japanese MEXT scholarship, the
Leids Kerkhoven-Bosscha Fonds (LKBF17.0.002) and also thanks the support and
hospitality of the Joint Institute for VLBI ERIC (JIVE).

Images made with the NASA/ESA Hubble Space Telescope are from the Hubble
Legacy Archive, a collaboration between the Space Telescope Science Institute
(STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA),
and the Canadian Astronomy Data Centre (CADC/NRC/CSA).

## References

  * Akashi & Soker [2008] Akashi, M. & Soker, N., 2008, New Astronomy, 13, 157
  * Balick [2003] Balick, B., 2003, Symbiotic Stars Probing Stellar Evolution, eds. R.L.M. Corradi et al., ASP Conf. Ser., 303, 407
  * Balick & Frank [2002] Balick, B. & Frank, A., 2002, ARA&A, 40, 439
  * Balick et al. [2020] Balick, B., Frank, A., Liu, B, Huarte-Espinosa, M., 2020, ApJ, 843, 108
  * Bondi [1952] Bondi, H., 1952, MNRAS, 112, 195
  * Bryce et al. [1997] Bryce, M., Lopez, J.A., Hollowat, A.J., & Meaburn, J., 1997, ApJ, 478, 161
  * Bujarrabal et al. [2016] Bujarrabal, V., Castro-Carrizo, A., Alcolea, M, Santander-Garcia, H., et al., 2016, A&A, 593, A92
  * Choudhuri [1998] Choudhuri, A. R., 1998, The physics of fluids and plasmas : an introduction for astrophysicist
  * Chen et al. [2016] Chen, Z., Frank, A., Blackman, E.G., & Nordhaus, J., 2016, MNRAS, 457, 3219
  * Clark et al. [2014] Clark, D.M., Lopez, J.A., Edwards, M.L., et al., 2014, ApJ, 148, 98
  * Clyne et al. [2015] Clyne, N., Akras, S., Steffen, W., Redman, M.P., et al., 2015, A&A, 582, A60
  * Corradi [2003] Corradi, R. L. M., 2003, Symbiotic Stars Probing Stellar Evolution, eds. R. L. M. Corradi et al., ASP Conf., 303, 393
  * Corradi & Schwarz [1993] Corradi, R.L.M., & Schwarz, H.E., 1993, A&A, 268, 714
  * Corradi & Schwarz [1995] Corradi, R.L.M., & Schwarz, H.E., 1995, A&A, 293, 871
  * Corradi et al. [2001] Corradi, R.L.M., Livio, M., Balick, B., et al., 2001, ApJ, 553, 211
  * Corradi et al. [2011] Corradi, R. L. M., Balick, B., & Santander-García, M., 2011, A&A, 529, A43
  * de Marco [2009] De Marco, O., 2009, PASP, 121, 316
  * Gómez et al. [2008] Gómez, J.F., Suárez, O., Gómez, Y., et al.. 2008, AJ, 135, 2074
  * Gómez et al. [2015] Gómez, J., Suárez, O., Benjoya, P., et al.. 2015, ApJS, 799, 186
  * Hyung & Aller [1996] Hyung, S. & Aller, L.H.. 1996, MNRAS, 278, 551
  * Herwig [2005] Herwig, F.. 2005, ARA&A, 43, 435
  * Hora & Latter [1996] Hora, J.L., & Latter, W.B.. 1996, ApJ, 461, 288
  * Kwok & Hsia [2007] Kwok, S., & Hsia C.H.. 2007, ApJ, 660, 341
  * Landau & Lifshits [1959] Landau, L.D. & Lifshitz, E.M.,1959, Pergamon Press, London
  * Mellema & Frank [1995] Mellema, G. & Frank, A., 1995, MNRAS, 273, 401
  * Mestel [1954] Mestel, L.,1954, MNRAS, 114, 437
  * Miranda et al. [2001] Miranda, L.F., Gómez, Y., Anglada, G., & Torrelles, J.M., 2001, Nature, 114, 284
  * Miszalski et al. [2015] Miszalski, B.; Manick, R.; McBride, V., 2015, The Physics of Evolved Stars, eds. E. Lagadec et al., EAS PS, 71-72, 117
  * Perinotto et al. [2004] Perinotto , M., Schönberner, D., Stefffen, M. , & Calonaci, C., 2004, A&A,414, 993
  * Sahai et al. [1999] Sahai, R., Dayal, A., Watson, A.M., et al., 1999, AAJ118, 468
  * Sánchez-Contreras et al. [2002] Sánchez-Contreras, C., Desmurs, J.F., Bujarrabal, V., et al.,2002, A&A, 385, 1
  * Santander-García et al. [2004] Santander-García, M., Corradi, R. L. M., Balick, B., Mampaso, A., 2004, A&A, 426, 185
  * Schonberner et al. [2005] Schönberner, D., Jacob, R., Steffen, M., et al., 2005, A&A, 431, 963
  * Soker [2005] Soker, N., 2005, AJ, 129, 947
  * Uscanga et al. [2014] Uscanga, L., Gómez, J.F., Miranda, L.F., et al., 2014, MNRAS, 444, 217
  * Vaytet et al. [2009] Vaytet, N.M.H., Rushton, A.P., Lloyd, M., et al., 2009, MNRAS, 398, 385
  * Welch et al. [1999] Welch, C.A., Frank, A., Pipher, J.L., et al., 1999, ApJ522, L69
  * Werner & Herwig [2006] Werner, K. & Herwig, F., 2006, PASP, 118, 183