Cylindrical Planetary Nebulae. I. Flow from an Irradiated Ring

: Many bipolar nebulae with a pronounced cylindrical shape, such as Henize 3-401, show no indication whatsoever of interaction between a disk and a stellar wind, or a jet on the nebular axis. I propose that the disk that is observed at the base of the bipolar is itself the source of the outﬂow. In particular, I assume that irradiation from the central star causes the disk to evaporate. I have performed numerical hydrodynamical calculations of outﬂows driven by evaporation of a pseudo-barotropic ring around a hot central star. The ﬁrst results show that the outﬂow shapes are cylindrical, and the internal structures are similar to what is observed in some of these nebulae. Since shape is only the ﬁrst step in the assessment of a model, synthetic observations should be made. For the moment I merely verify that the scalar quantities observed in the archetypical cylindrical nebula Hen 3-401 can be accommodated in my models.


Introduction
Most (pre)planetary nebulae are very aspherical, and many have a distinct bipolar shape. In the past decades, many authors have proposed mechanisms for the formation of such nebulae: interaction between a stellar wind and a circumstellar disk, supersonic jets emanating from a compact disk, and magnetic fields of various configurations. However, some bipolar nebulae with a pronounced cylindrical shape [1], such as in Figure 1, show no indication whatsoever of interaction between a disk and a stellar wind, or a jet on the nebular axis. Moreover, in these cases there is clear evidence of a thick disk at the base of the outflow, and it would be quite remarkable if the outflow did not actually start there. Since the disk has no internal source of energy, the power must be supplied from elsewhere, and the radiation from the central star is the only candidate. Therefore, the cyclindrical shape is very suggestive of evaporative outflow. This prompted me to compute a series of hydrodynamical models in which a dense gaseous ring around a hot luminous star is irradiated so fiercely that the disk surface boils off into space. The disk is presumably approximately barotropic, so that its cross section is quite thick; not really a torus, but more akin to a hollow ring. I do not specify an origin for the disk-ring; a plausible possibility is that the central object is a late post-AGB star, on its way to becoming a white dwarf, surrounded by the remnants of its planetary system (or a super-jupiter binary companion), and/or a gas-dust mixture slowly ejected in an earlier phase of the star's evolution.

Evaporation Outflow Model
In previous studies [3], I have considered the formation of bipolar nebulae under the assumption that the shaping of the nebula is due to a spherical stellar wind impinging on a disk around the star. I found that this mechanism makes nice bipolars, and that it is possible thereby to explain quite a range of observed nebulae. Others have presented computations based on jet-driven or magnetically confined outflows (e.g., Ref. [4] and references therein).
However, there are many bipolar nebulae that do not show the figure-eight shape typical of such interaction models [5]. In many of these objects the circumstellar disk is quite prominent. Considering this, I decided to look at the launching problem 'the other way around' and assumed that the flow originates directly on the disk, instead of on (or very near to) the star. When a spherical outflow is confined by a disk, it is quite hard to make tall cylindrical outflow shapes, due to the contrast between the Mach number of the stellar wind and that of the circumstellar disk. Therefore, I decided to investigate what the opposite viewpoint offers: namely, that the disk itself is doing the launching.
Suppose that the shape of the ring-shaped disk around the star is given in cylindrical coordinates (r, φ, z), with a cross section z(r). Consider an infinitesimal annulus with a surface area ∆A on the ring. Suppose that this annulus is evaporating due to irradiation from the star. Let the velocity v be perpendicular to the ring surface, with gas density ρ. The amount of energy that flows per second through that surface, with speed v, is 1 2 ρv 3 ∆A. This energy is delivered by the radiation from the star which has luminosity L. If the surface ∆A subtends a solid angle ∆Ω as seen from the star, then the energy flux that corresponds to the above kinetic energy is L ∆Ω/4π. Equating the two energy fluxes we obtain The annulus has an inclination ψ − θ with respect to the direction towards the star, where ψ is the inclination of the ring surface with respect to the symmetry axis, and θ is the latitude in spherical coordinates (R, φ, θ). Therefore, in which R = √ r 2 + z 2 is the radial distance from the star to the ring surface (I assume that the radius of the star is negligible). Putting all of this together, we obtain with the stipulation that v = 0 when ψ < θ, because then the annulus is in the shadow of the ring itself. This equation is valid under the assumption that all of the radiative energy is converted into kinetic energy, hereinafter called 'Case A'. However, some of the radiation will be used to heat the surface of the ring. The temperature T(R) of the outflowing gas is determined by the stellar irradiation. I assume that the star and the ring radiate as black bodies, so that the temperature is found from the classical relation T 4 R 2 = const ∝ L. When setting up the outflow boundary conditions, the temperature is not used as such, but enters by way of the gas pressure P. Because the gas density is assumed to be constant on the ring surface, the black-body condition implies that P = P 0 √ R 0 /R in which P 0 and R 0 are fiducial quantities that depend on the local geometry. As an aside, note that this means that the speed of sound on the evaporating surface scales as s ∝ R −1/4 . The thermal energy density is the gas pressure P. Therefore, Equation (3) has to be replaced by This description of the launching conditions on the ring will be called 'case B'. In order to prepare this for the numerical computations, I convert the equation to dimensionless form. Each variable V is replaced by V 0 V, in which V 0 is a fixed fiducial dimensional quantity, and V is now dimensionless. This ploy will show more clearly on what quantities the simulation depends. Writing v 0 v, ρ 0 ρ, P 0 P and R 0 R, I obtain Factoring out the fiducial quantities gives This equation may be rewritten in a more transparent form by using in which s 0 is the fiducial speed of sound, γ the adiabatic index of the gas, M 0 the Mach number, and L 0 the energy flux. Equation (6) then becomes This is an implicit equation for the velocity v(R) in Case B. Using Equation (3) it follows that closer to the star (smaller values of the dimensionless radius R) the kinetic energy density E = 1 2 ρv 2 dominates over the thermal energy density P. The dimensionless equation that corresponds to Case A, Equation (3), is When setting up the initial conditions, I prescribe the shape of the ring cross section, from which I compute ψ(R) and θ(R). In Case A, I use Equation (3) to compute v(R) for plausible values of ρ, M 0 and L/L 0 . All of these dimensionless quantities are of the order of unity. Numerical experience shows, for example in plots such as Figure 2, that in the inner regions of an evaporating disk the gas pressure does not play a significant role. This is a consolation for the fact that in the present preliminary investigation I have not computed the actual ring structure in any detail. That is why I have used Equation (3) to find the launching speed, and have proceeeded from there.
The initial conditions in the gas that surrounds the ring are even more uncertain. It is often assumed that the gas density obeys the proportionality ρ ∝ R −2 , as is the case in a spherical outflow with constant speed. If the temperature of the gas is due to the radiation from the star, the pressure would behave as it does on the surface of the ring. In that case, the speed of sound would scale as s ∝ P/ρ ∝ R 3/4 . I will return to this point in my conclusions.

Shape of the Ring
Determination of realistic disk-wind launching conditions is a very complicated exercise in radiative transfer. In order to make a start, I have concocted a plausible configuration that I prefer to call a 'ring' rather than a disk, because it is quite thick (as may be seen in the observed images of relevant cases, e.g., Figure 1). Assuming that the configuration is cylindrically symmetric, and taking (r, z) to be the radial and axial cylindrical coordinates, I propose the function in which r − and r + are the inner and outer equatorial radii of the ring, respectively; the number n governs the concavity of the ring, typically n ≈ 2; and b is an arbitrary constant of order unity. This shape may be seen as a plausible approximation to the outer equidensity surface of a barotropic disk. The surface of the ring is irradiated by the star. The resulting evaporation flow starts with a velocity that is perpendicular to the ring surface, and has a magnitude such that the flux received from the star equals the outflow energy flux (case B) or the outflow energy minus the thermal energy due to the radiative heating (case A). The inclination ψ(r) of the ring and the latitude θ(r) are determined by The speed is found from Equation (3), and because the outflow is assumed to be perpendicular to the ring surface we obtain the launching velocity in cylindrical coordinates (r, φ, z) (the azimuthal velocity is zero).
In the numerical solutions to be presented below I have adjusted the values of b and n, as well as the ratio r + /r − , until the ring shape looked more or less plausible; that is to say, having a cross section that is sufficiently concave (such that the evaporation velocity is not directed inward too much) but not too thin. Typical shapes and velocity fields are shown in Figure 3. At this stage of the investigation, I have no better argument for adopting these values than noting that b ≈ 1 and n ≈ 2 produce cylindrical chimney-type outflows that resemble the observations rather well.

Hydrodynamical Results
In order to test whether star-driven evaporation produces cylindrical nebulae, I have performed numerical hydrodynamical calculations of outflows driven by evaporation of a pseudo-barotropic ring around a hot central star. My hydrodynamical codes have been presented in [6]. Further improvements that I have made since then do not change the key features of my FCTLCD2D hydro solver, and are almost all concerned with upgrading data handling and computing speed (notably vectorization and parallel loop execution). These are too straightforward to merit separate publication.
To obtain quick results that may indicate whether the chosen approach is fruitful, I have implemented the above initial conditions in my hydrocode FCTLCD2D. I have simplified the setup to a constant gas pressure instead of the dependence in Equation (4). The gas surrounding the ring is assumed to be homogeneous, as opposed to the dependence ρ ∝ R −2 mentioned in Section 2. The boundaries of the computational grid were fitted with reflectionless outflow conditions.
The images below show a very small sample of the resulting output. All of the results display a pronounced cylindrical outflow pattern. Many of these show a striking similarity to the shape of Hen 3-401 and its kin ( Figure 4). It came as a very welcome surprise that even the first runs produced shapes that match the observations rather well. In all the images the [RGB] colour coding is what I habitually use in my work: R = gas density ρ, G = pressure P, B = speed v. In some cases I have increased the contrast and range of the lower values, to make the details of the outflow more prominent. In all of my simulations, the internal configuration of sequential X-shocks on the axis is beautifully steady, just as in the case of the classical rocket exhaust nozzle. Under astrophysical circumstances, the surrounding gas (the stellar wind deposited before the evaporation flow started) is much less dense than that of the outflow, whereas in the laboratory setup the atmospheric pressure and density are comparatively high. Therefore, the outer walls of the evaporation flow show vortex structures where the outflow mixes with the ambient gas. Even though the shapes of my outflows seem to be suitably similar to what is observed, there is much more involved than just shape.
A typical example of the outflow velocity field is shown in Figure 5. The turbulent eddies on the periphery should be observable. The streaming velocity in the axial direction (middle panel) appears to be constant after the initial acceleration away from the ring. It would seem, however, that in many cases observations show precisely the opposite, namely, a velocity that increases approximately linearly with distance from the star. Unless there is an extra source of acceleration (such as the stellar radiation driving the gas after it has left the ring surface), I am afraid that these observations are in tension with the present model. Other flow features that should be compared with observations are identified in Figure 6. After completion of this series, I have investigated the consequences of the difference between Equation (8) (Case B) and Equation (9) (Case A) for the shaping of the nebula. In Case B, I assume that the temperature, and therefore the pressure, of the evaporating surface is set by the irradiation from the star. Therefore, the energy available is divided between kinetic and thermal energy. In the latter case I assume that the pressure is given somehow, for example by the hydrostatic conditions in the ring, so that all of the stellar energy is turned into the kinetic energy of the evaporating flow.
A typical example is shown in Figure 7. The 'kinetic-only' Case A (Equation (9), top panel) produces a cylindrical nebula. The nebula in the 'kinetic-plus-thermal' Case B (Equation (8), bottom panel) is much more complex. The reason for this is that close to the star, the kinetic energy of the outflow always dominates (compare Figure 2); however, near the place where the stellar light grazes the ring, the speed drops to zero, whereas the pressure does not.  Consequently, the place near the terminator on the ring (beyond which the star cannot be seen) always produces a dense slow outflow. Further inward, this gas meets the evaporating flow closer to the star and is dragged away from the ring. That interaction makes the flow that intersects on the axis very complicated. Further out, this entrained gas stays in the outer shell of the nebula, leaving a fast low-density stream near the axis. Likewise, the secondary shock that arises when the inward curving shell bends back towards the axis is very different from the more or less ordinary diamond-shock of the case described by Equation (9) (cf. Figures 4 and 5).
The influence of the ring concavity on the outflow conditions may be seen in Figure 3. The hollower the ring, the more radiation is intercepted by the outer regions. Closer to the star, the light impinges in a more grazing direction, wherefore it is effectively less intense; on the other hand, the surface curvature focuses the initial outflow more in the direction of the symmetry axis. I have investigated what consequences the ring concavity has for the resulting flow pattern. Some results are shown in Figure 8. A bowl-shaped ring produces the widest outflow cylinder, whereas the ring with a conical inner surface makes the flow narrower and almost jet-like.

Application to Henize 3-401
To connect with the available observations, I have to map the computed quantities onto the properties that are observed in actual nebulae. A case in point is Henize 3-401, because the shape of that nebula is very suggestive of evaporative outflow. The images show no indication of interaction between the disk and a stellar wind, or a jet. A stellar wind, confined by the disk, would produce a bipolar nebula, but with the shape of the classical Figure 8, and not the elongated cylinder observed in Hen 3-401 and similar morphologies. A jet would be preceded by a bow shock, which is not observed.
I have assembled the most recent values (as far as I can determine) of the observed properties of Hen 3-401 from [7][8][9][10]. Unless indicated otherwise, the calculations below all refer to the data listed in Table 1.

Gauging of Dimensionless Variables
My computations of evaporation flows from circumstellar rings were performed in dimensionless units. Each variable a is scaled by writing A = a/a 0 , where a 0 is a fiducial quantity to be determined for each specific case. The equations to be solved are cast in terms of the resulting dimensionless variables A. The primary scaling quantities are the mass density ρ 0 , disk radius R 0 , flow speed v 0 and luminosity L 0 . These must be related, because the only independent quantities are mass, length and time. Using Equation (7), we may compute the value of the mechanical luminosity L 0 for some appropriate values of the fiducial quantities ρ 0 , R 0 and v 0 .
Unless the distance to Hen 3-401 has been grossly overestimated, we may consider R 0 as fairly safe, but it is not easy to assess the mass density ρ 0 . The data from [7] suggest n H = 3 × 10 15 m −3 . Another way of estimating ρ 0 is by assuming that the nebular mass 0.01 M = 1.99 × 10 28 kg is distributed uniformly in a cylinder with the size of the nebula. For the above value of R 0 , the volume is πR 2 0 H = 1.9 × 10 46 m 3 , in which H ≈ 28 R 0 is the total length of the nebula. In that case, ρ 0 = 10 −18 kg m −3 , corresponding more or less to the lowest value reported [7], namely 10 9 m −3 , that is ρ 0 = 2 × 10 −18 kg m −3 . With these numbers, we obtain a particle density of 6.3 × 10 8 m −3 and L 0 = 4.5 × 10 27 W = 12 L .
Browsing through the observations, it seems clear that the determination of the length scale R 0 is probably the most robust. The speed v 0 is next, although its value varies a lot in space, and the inclination of the axis of Hen 3-401 with respect to the plane of the sky is somewhat uncertain. Because v 0 appears in the expression with a third power, a relatively small adjustment would produce a noticeable change in L 0 . However, ρ 0 is the least well known, and it seems reasonable to compute L 0 for a range of densities to see if there are values of ρ 0 that yield acceptable values of the mechanical luminosity. The result is that densities in the neighbourhood of 10 −15 kg m −3 (or 6 × 10 11 m −3 ) produce mechanical luminosities that more or less correspond to what the central star of Hen 3-401 might deliver. This falls comfortably in the range [7] from 10 9 m −3 up to 10 15 m −3 .
There is one more constraint on the mass density ρ 0 to be considered. Because the outflow is presumably driven by the irradiation from the central star, the density must be such that the distance R 0 corresponds at least to an optical depth of approximately unity. This places an upper limit on the mass density. That is to say, if the mean scattering cross section of the nebular matter is σ, then R 0 ≈ 1/n 0 σ so that ρ 0 ≈ m/σ R 0 where n 0 is the particle density corresponding to ρ 0 , and m is the mean particle mass. For example, using the Thomson cross section σ T , and taking R 0 = 6 × 10 14 m from the above, we obtain ρ 0 = m H /σ T R 0 = 4.2 × 10 −14 kg m −3 and n 0 = 2.5 × 10 13 m −3 . This upper limit is fairly consistent with the data in Table 1. All of this considered, I think that ρ 0 = 10 −15 kg m −3 = 6 × 10 11 m −3 and L 0 = 4 × 10 30 W = 10 4 L are plausible values.
One more parameter is essential in the hydro simulations, namely the initial Mach number M of the evaporating outflow. Because the circumstellar ring is supposed to have no energy sources other than the stellar irradiation, we may assume that M ≈ 1, so that v 0 is approximately equal to the local sound speed s. To determine s, we must know the local gas pressure P, which may be derived from the local temperature T by means of the ideal gas law: where k is Boltzmann's constant, m is the mean particle mass, A is the stellar surface, and σ is the Stefan-Boltzmann constant. In my evaporation model, T 0 is due to the central star.
If the central star is B1, then T = 2.6 × 10 4 K and L = 1.35 × 10 4 L = 5.2 × 10 30 W, and using r = 10 14 m we obtain which compares fairly well with the IRAS values [10]. Taking the IRAS temperature they report, 350 K, I obtain s = γkT/m H = 2.2 km/s which should be seen as a lower limit, because the gas temperature is expected to be larger than the dust temperature.
Finally, there is one more quantity that should be compared with the local speed of sound, namely the orbital speed v d of the disk at radius R 0 : v d = √ GM * /R 0 = 4 km/s if M * = 11 M and r 0 = 10 14 m (B1 star). If the disk is in hydrostatic equilibrium, the local sound speed should be about an order of magnitude less. We may conclude that the Mach number of the initial outflow would be in the range 1-5. Most of the simulations I have run have a Mach number in this range.

Discussion and Conclusions
My original aim was a preliminary investigation to find out if a circumstellar ring, irradiated by the central star, would produce a cylindrically collimated outflow. The conclusion is that it does, for a plausible family of ring shapes, and under quite a wide variety of launching conditions. The generic flow features, described in Figure 6, are very suggestive of features found in many pseudocylindrical bipolar nebulae, as far as the observational resolution allows. However, a morphological match is merely a minimum requirement of any theoretical model. This should be followed up by the construction of observables, which I will reserve for a future publication.
Two more avenues for future research are: the inclusion of rotation about the symmetry axis, and fully three-dimensional hydrodynamics. The azimuthal speed estimated above might not be very much smaller than the axial speed on launch. It is to be expected that the conservation of angular momentum would prevent the flow from focusing too strongly on the axis, and might contribute to the increase of the density on the cylinder walls. Simulations of rotating flows are readily possible in two dimensions, and I will work on that in the future. The 3D avenue is more technically demanding and will have to wait.
The main item of concern is the velocity field. In the majority of the computational runs, the outflow is only weakly accelerated. This is at variance with the observation that in many such nebulae, the speed of the outer walls increases approximately linearly with the distance along the cylinder. This shortcoming may be due to my initial choice for the distribution of the gas surrounding the ring, which I have given a constant and rather low mass density. As mentioned in Section 2, the gas density in a spherical stellar atmosphere expanding at constant speed obeys ρ ∝ R −2 . If the temperature of the atmosphere is due to the radiation from the star, the speed of sound would scale as s ∝ P/ρ ∝ R 3/4 , so that the Mach number would decrease away from the star: M ∝ R −3/4 . If the density of the gas into which the outflow propagates is much larger, and if indeed ρ ∝ R −2 , the outer walls would be much denser, and would couple to the ambient gas through turbulent viscosity (cf. the eddies in Figure 5). These effects together might produce an acceleration with v ∝ z, in accordance with the observations (see also [11]). To make matters even more complicated, one should bear in mind that the late stages of stellar evolution on the red giant branch produce outflows that become faster in the course of time, and may even be episodic. A possible consequence, proposed by an anonymous referee, is that a secular decrease of the stellar luminosity leads to a decrease of the launching speed.
Published observations do not readily provide the main physical quantities needed for theoretical modeling. In particular, it is not clear to what part of the nebula the various quantities refer. Moreover, the error margins are difficult to assess. I am not qualified to judge the quantitative aspects of the observations, but purely from the nebular shapes I would say that, for example, Hen 3-401, He2-437, OH231.8 and V445Pup are fitted very well by my models, without undue effort in adjusting the outflow parameters.
My model has at least this virtue, that the source of the driving energy (the light from the central star) and the cause of the collimation (the hollow inner surface of the ring) are evident (as opposed to, e.g., jet-driven models, where the mechanism that causes the jet in the first place is often unspecified). This is obviously not enough to be an incontrovertible explanation for nebulae of the Hen 3-401 type, but it is a start.