Models of the Hydrodynamic Histories of Post-AGB Stars. I. Multiflow Shaping of OH 231.8+04.2

We present a detailed hydrodynamic model that matches the present structure of the well-observed preplanetary nebula (“pPN”) OH 231.8+04.2 (“OH231”). The purpose of the model is to present a physically justified and coherent picture of its evolutionary history from about 100 years from the start of the formation of its complex outer structures to the present. We have adopted a set of initial conditions that are heavily constrained by high-quality observations of its present structure and kinematics. The shaping of the nebula occurs while the densities of the flows are “light,” i.e., less than the surrounding AGB-wind environment. The simulations show that pairs of essentially coeval clumps and sprays of the same extent and density, but different outflow speeds, sculpted both the pair of thin axial flow “or spine” and the bulbs. The total ejected mass and momentum in the best-fit model are surprisingly large—3 M⊙ and 2.2 × 1041 gm cm s−1, respectively—however, these values are reduced by up to a factor of 10 in other models that fit the data almost as well. Our ultimate goal is to combine the present model results of masses, momenta, flow speeds, and flow geometries for OH231 with those of other models to be published in the future in order to find common attributes of their ejection histories.


Introduction
The traditional picture of AGB winds and the formation of protoplanetary nebulae ("pPNe") is that they are driven by isotropic radiation pressure absorbed on dust particles that condensed in expanding and cooling quasi-isotropic winds. It is clear that in most cases reality is far more interesting, albeit still only partially understood. Dozens of optical images of pPNe from the Hubble Space Telescope 4 ("HST"), as well as images in molecular emission lines and dust-emitted continuum light from millimeter-wave interferometers such as ALMA, show that the wind streamlines are far from isotropic. Other studies (e.g., Bujarrabal et al. 2001, hereafter "B+01") have demonstrated that the momentum in the outflows is often far too large to be explained by the winds of a single AGB star. Moreover, the predominant kinematic pattern of pPNe is one in which Doppler shifts increase steadily with radius from the star-in direct conflict with expectations of steady jets from stellar outflows.
Observing the active flow acceleration and collimation process presently lies far beyond our grasp. The challenges are that we are generally unable to catch the brief ejection event at its start or to monitor its subsequent evolution using telescopes that spatially resolve the structure of the outflow in its infancy as it grows into its final form. So, the nature of the ejection and collimation processes must be inferred. The most plausible general explanations for the various types of bipolar pPNe-and those for which there is slowly mounting evidence -is that pPNe are formed by mass ejected in magnetically shaped winds from the active surfaces of AGB stars or from disks or mass overflows associated with mass transfers in close binary systems (Balick & Frank 2002;Nordhaus et al. 2007;De Marco 2009;Soker & Kashi 2012;Chen et al. 2016;Staff et al. 2016;Iaconi et al. 2017).
In this paper, we present hydrodynamical models that successfully fit extant observations of one pPN, the appropriately well-observed OH 231.8+04.2 (hereafter "OH231"). However, as the models will show, the observably expanding structures are presently coasting long after the completion of their active shaping process. Thus, in order to surmise its evolution backwards from the present, we must use hydrodynamical models of two pairs of axisymmetric flows to provide insight into the early ejection process. OH231 is one of the most advantageous of its class for modeling owing to its large angular size, well-resolved and regular shape, and a rich literature of observational data.
In particular, we will show that the models yield a simple and robust set of tightly constrained initial values for the hydrodynamic state variables that describe the evolution of OH231 and the mass and momentum of the outflows that have shaped it. The data that we used to constrain our hydrodynamic simulations will be discussed in Section 2. In Section 3, we describe our methodology and the starting point for the simulations. The results of the models are presented in Section 4. In Section 5, we summarize the results and their relevance for understanding the histories of pPNe. Riera et al. (2005, hereafter "RRA05") constructed a hydro model for OH231 consisting of a pair of non-axisymmetric and brief columnar jets, or "jet-like pulses," injected at 200 km s −1 . The northern pulse passed close to a companion AGB star with a continuous asymmetric wind that is located near the end of the lobe. The interaction of the pulse with the outer parts of the AGB wind impeded its progress. We will compare the two sets of results later in the paper.

Relevant Observations
We summarize the observations of OH231 that are used to contract and constrain the present model simulations, to wit, the current spatial structure of the hydro state variables density, vector momentum, and temperature. The associated proxy observables are the brightness distribution of the images, patterns of Doppler shifts (Figure 1), and the 12 CO kinetic temperature. In the discussion that follows, and unless otherwise noted, we summarize the characteristics of this iconic and relatively massive pPN compiled from papers by Kastner et al. (1992, hereafter "S+00A"), Sánchez Contreras et al. (1997, 2000, hereafter "S+00B"), Alcolea et al. (2001, hereafter "A+01"), and Bujarrabal et al. (2002, hereafter "B+02"). 5 Geometry. OH231 consists of two opposite flame-shaped "lobes" of similar shapes (Figure 1, yellow and white colors) that are enclosed within bulbous structures called "bulbs." Of vital importance is the high-resolution optical and near-IR images of OH231 obtained using HST. 6 All of the extended structures are outlined in shock-excited lines (Figure 1, blue colors; S+00A&B). The leading tip of the northern (southern) lobe is 18″ (38″) from the dark lane that separates them. The lengths of the 12 CO outflow lobes are slightly smaller. The distance to OH231 is 1.5 kpc, so the tip of the north (south) lobe lies at 32,000 (66,000) au from the nebular core after correcting for the nebular inclination of 35°(a factor of 1.22).
The lobes and bulbs extend from a dense, compact, highly structured, slowly expanding, and only partially resolved central core. All of them share a common symmetry axis. In addition, A+01 make a compelling argument that the core observable features of OH231 are surrounded by an isotropic, stratified, and low-surface-brightness density distribution consistent with a wind from a central AGB star whose origin preceded the ejection of the lobes and bulbs by thousands of years. See A+01 (their Section 3.2.1) for a discussion of the history of the AGB winds. It is commonly assumed that the AGB wind is steady. If so, then its radial density profile declines as r −2 . Figure 1. Images of OH231 in optical light (left) and 12 CO (right) taken from the literature. Left: an overlay of optical and near-IR images of OH231 obtained from the Hubble Legacy Archive. North is up. The white circle provides an image scale. Right: spatial distribution and axial Doppler speeds of 12 CO lines from OH231 obtained at the Plateau de Bure interferometer and republished with permission (CO panel source: a Web document entitled "From the AGB to the PN phase with the Plateau de Bure Interferometer" posted as a PDF file by A. Castro-Carrizo with J.M. Winters and R. Neri, c. 1998, with no public URL. These CO results were subsequently reformatted and republished by A+01). The upper (lower) row shows the images and kinematic results for the 12 CO(1-0) ( 12 CO(2-1)) line.
Aggregate properties. The total mass of OH231 inferred from 12 CO observations is 1 M e (A+01), one-third of which is in the two spines (the fast flows) that lie along the nebular symmetry axis ( Figure 1) and rest is slow gas near the core. The remainder lies in the slow AGB wind and core. The shorter and slower northern spine has slightly more mass and momentum (<15%) than its southern counterpart (A+01). Given their similar masses but differences in their speeds, the slower, brighter northern spine presumably contains about twice the momentum as its southern sibling.
The total linear momentum of OH231 seen in 12 CO, ≈27 M e km s −1 , is about 10 5 times larger than the linear momentum that the central star of OH231 can supply (B+01). Thus, radiation pressure supplies a trivial fraction of the momenta that formed the surrounding structures. Such large momentum "surpluses" are rare among pPNe, but not unprecedented.
Bulbs. Bulbs surround both the north and south lobes whose best simulations' dimensions are about the same size as the lengths of the lobes. The southern lobe is by far the most conspicuous, having dimensions of roughly 27″×30″ (40,000×45,000 au). The vertex at the base of the bulbs is about 4000 au in width. The inner half of the southern bulb has straight edges that open at angles of ∼±50°before they reconverge at the protruding tip of the southern lobe. This lobe shows signs of ripple-like instabilities along its forward edge. Hα lines inside the lobe exhibit a trend of monotonically increasing Doppler shift with distance from the core with speeds as large as 400 km s −1 . A high-contrast rendition of the Hα image shows that the northern bulb is thinner than the southern bulb. Its aspect ratio is about 1:2 and its opening angle is about 40°at the base.
Lobes. The lobes are flame-shaped in outline. (This morphology is common among the lobes of pPNe.) Each lobe is comprised of a spine plus a skirt ( Figure 1). The brighter northern lobe has a rounded leading edge within which a knot of HCO + (a tracer of high-density gas 10 6 cm −3 ) is located. In contrast, there is no prominent peak in any molecular tracer or starlight-scattering dust near the end of the southern lobe. Rather, the spine slowly fades and becomes irregular starting about halfway from the core. Forde & Gledhill (2012) found patchy areas of H 2 emission near the leading edge of the northern lobe. H 2 is typically emitted in shocks of speeds of a few tens of km s −1 . Thus, the H 2 lines may arise behind the lobe tips where putative outflows overtake and shock the leading edge of the bulb at modest speeds (20-50 km s −1 ).
Spines. Thin, bright spines run down the nebular symmetry axis in both lobes ( Figure 1). Spines have a central relevance here, since of all of the features of OH231, they alone have clearly defined kinematic patterns. The overall shapes of both spines are similar. 12 CO within each lobe is intimately associated with its spine (but not the skirts). 12 CO observations show that the tip of the southern (northern) spine has an inclinationcorrected outflow speed of about −400 (+200) km s −1 . The outflow Doppler shift of the spines rises linearly with distance. This kinematic pattern is sometimes called a "Hubble-like flow pattern" that has been interpreted by A+01 and RRA05 as a sign of a ballistic jet-like outflow. (We will offer a very different alternative interpretation later.) The spines share the same speed gradient and, thus, the same kinematic ages, ≈770 yr.
Southern Skirt. The skirt ( Figure 1) is seen only in reflected starlight. Thus, its density, kinematics, and mass are not directly measurable. The wedge-shaped lateral edges of the skirt suggest that the skirt may be starlight that is scattered from dust lying within a stellar illumination cone of opening angle±20°. More speculatively, the skirt may be a dusty conical outflow; however, no swept-up gas is found along its outer edge.
In addition, the central stars of OH231 and their possible roles in the ejection of the nebular features of OH231 are discussed by Sánchez Contreras et al. (2004). Of particular relevance here is that the AGB star, QX Pup, has an A-type companion. The binary originally has a total mass in excess of 4.5 M e . This is an upper limit on the nebular mass today. They also argue that the momentum of the nebula is supplied by winds from an accretion disk of the companion with an orbital radius 50 au. The ejection of the gas took place within 100 yr of its onset. The mass-loss rate needed to supply the fast outflows was of order 10 −3 M e yr −1 .

Methodology
The code AstroBEAR (Cunningham et al. 2005(Cunningham et al. , 2009) was used to construct hydrodynamic models of OH231. Astro-BEAR solves the Eulerian equations of fluid dynamics on a 3D grid using high-performance adaptive mesh refinement as needed to resolve shapes and shocks. The kernel of AstroBEAR is a versatile 3D hydrodynamic solver that can be applied in a wide variety of circumstances, including magnetized/cooling flows and rotating coordinate systems (Huarte-Espinosa et al. 2012). The code does not contain facilities for simulating emission-line images, predicting an emergent spectrum, or for following the transfer of radiation.
The numerical code accommodates three types, or paradigms, of mass injection from a central source into the ambient gas: clumps (a spherical knot of a single ejection age), cylindrical jets (hereafter "jets": a thin flow consisting of parallel streamlines ejected at constant speed), and steady tapered conical flows or "sprays" (diverging streamlines as described later). Each of these can be launched into the AGB environment in any order. In all models, we assume that the ambient gas consists of AGB ejecta and that the flows are axisymmetric and share the same symmetry axis and a common origin ( Figure 1).
We performed the computations in one quadrant of the 64×256 2D model grid (32,000×128,000 au) in which the injected gas is launched into a stratified ambient medium along the y-axis at time t=0. However, computational speed is important since multiple free parameters must be adjusted for a successful fit. Thus, we did not attempt 3D modeling because the amount of computing time per simulation is prohibitive. The computational complexity of 3D models is only needed to follow the growth patterns of small-scale instabilities at shock interfaces. This is not our primary interest.
The basic grid size, 500 au, was decreased by the adaptive mesh algorithm by up to five factors of powers of two in regions where local pressures have steep gradients. Given the number of free parameters ,we ran over 60 different models at modest resolution to find a good low-resolution outcome. Subsequently, we reran many of the most successful simulations, enhancing the resolution by one or two additional factors of 2. Even so, these "full" models often required 24 hr to execute. Thus, after the flow configuration had grown substantially and attained its final form (after a few hundred years), the resolution was relaxed by a factor of two with no noticeable loss of overall morphology.

Flow Concepts
We consider the most appropriate flow concepts for forming the bulbs and the spines seen in HST images. The bulbs are formed by a wide-angle, geometrically tapered polar flow called a "spray" with an opening angle consistent with the HST images (see below). Narrower flows, such as jets and clumps, were considered for the formation of the spines. The age of OH231 is rounded to 800 yr.
Ambient Gas. We assumed that the ambient AGB mass deposited on the grid prior to t=0 is isotropic with a density distribution n amb (r)=n amb,0 (r 0 /r) 2 . Here, n amb,0 is the density of the ambient medium at a fiducial radius r 0 . We adopted n amb,0 =2×10 4 cm −3 at r 0 =3500 au. The total mass of the ambient gas beyond 1000 au is ∼0.7 M e , in accordance with the estimate of A+01. Our simulations are restricted to a 2D slice through this density distribution.
Bulbs. The bulbs are presumably the outcomes of conical nuclear sprays that have displaced the ambient gas within their volumes. The spray exits a nozzle at an opening angle Θ spr with density n spr , speed v spr , and kinetic temperature T spr . The density and speed are modulated in polar angle by a Gaussian of 1/e width f spr , where f spr ≈Θ spr . This modulation of the flow momentum "softens" the edges of the flow. (In the absence of this modulation, a conical flow forms an open wedge-shaped lobe with straight sides when injected into an ambient medium with the same radial density distribution.) Lee & Sahai (2003, hereafter "LS03") showed that some sprays form closed flame-shaped lobes.
Spines. The spines are relatively thin axial features whose shapes and kinematics are defined by 12 CO observations that were presented by A+01. The spines exhibit a linear kinematic gradient seen that is fundamentally inconsistent with any type of steady collimated flow such as a jet. For this reason, we consider that the spines are similar to the wakes of clumps ejected in the same event as the sprays.
The clumps are characterized by their properties at the nozzle upon ejection at t=0: an initial radius r cl , initial density n cl , ejection speed v cl , and temperature T cl . These are free parameters though their starting values are extracted from the observations in Section 2. The mass of the clump at launch is 1.54×10 −8 (n cl /cm −3 ) (r cl /10 3 au) 3 . A+01 estimated that the total mass of high-speed gas (i.e., the spines) is on the order of 0.2 M e . We also adopt T cl =10 K; this value is of no consequence to the outcomes so long as T cl <10 5 K.
Placement of the nozzles. The bulb must be launched into the lee of the spray in order for the wake that follows the clump to remain intact. If their locations are reversed, then a spray of appropriate total momentum flux to account for the bulbs will destroy the structures in the lee of the clump and eradicate the spine. The details are discussed within the context of specific models.
A note about alternate paradigms: A+01, Meakin et al. (2003), RRA05, and others have speculated that cylindrical jets are the most viable paradigm, provided that the jets are accelerated in situ in such a way that they emulate a ballistic jet. The in situ acceleration mechanism is unspecified and we are not able to imagine what it might be.

Results
The best-fit simulation of the flow characteristics of OH231 is shown in Figure 2 panel A. Table 1 lists the initial conditions that were used for the model. The major features of the model nicely match the observations, including the speed gradient observed along both spines.
The densities of the clumps and sprays at the launch surface n cl,0 and n spr,0 ; their launch speeds v cl,0 and v spr,0 ; the radii of the clumps r cl ; and the nozzle widths, geometries, and durations of the sprays, r spr , Θ spr , and f spr , are, in principle, free parameters. Their initial values were estimated as follows: v cl,0 and v spr,0 were set equal to the observed speeds at the end of the spines: 225 (400) km s −1 to the north (south). The clump radius was set to the size of the circle in Figure 1; r cl =3500 au Θ spr and f spr were directly measured from the Hα images of OH231 at the base of the flow and set equal to one another. The densities n cl,0 and n spr,0 were found by trial and error for each lobe.
The northern and southern sprays and clumps are launched at t=0 from two different spherical nozzles that are aligned on the y-axis. The nozzle for each spray has a radius of 5000 au so that the spray initially envelops the clump on the grid (see Section 3.2) and precedes it into the ambient gas (until the clump overtakes the leading edge of the spray after several hundred years). After 800 yr, we estimate from the initial conditions that the total mass (momentum) of gas injected into the spray and clump is 2.8 M e (∼1000 M e km s −1 =2×10 41 g cm s −1 ).
The "best model" in Figure 2 panel A is simply one of a family of very similar simulations that nicely match the shapes and kinematics of OH231 (Figure 1). Other members of the family, or "isoflows," share the same initial parameter values except for their densities. The ratios of all densities, n amb , n spr , and n cl , are the same for all of the isoflows. The "best model" is simply the isoflow that satisfies the 12 CO mass requirements of OH231: mass of slow (fast) gas ≈0.7 M e (0.2 M e ). The bulb is not detectable in 12 CO.
However, the masses and momenta of the flows that shape the spine and bulbs in the best-fit model seem excessive (∼factor of 10-100) relative to their values in other pPNe (B+01, Section 2). Specifically, B+01 estimate a combined momentum for the spines of OH231 of 4×10 39 g cm s −1 (with very large uncertainties). Our result for the spines where the 12 CO lines arise agrees reasonably well with this value. The momentum of the bulbs dominates the combined total in the model. This momentum is plausible because the volumes of the spray-formed bulbs are about 100 times larger than the clump-formed spines. Thus, the sprays have displaced up to 100 times more ambient mass (and momentum) than the spines.
For comparison, the total mass and momentum of the isoflow shown in Figure 2 Panel B is 10 times smaller than the best-fit model even though both simulations have the essentially the same shapes, dimensions, and kinematics. However, this simulation is not quite as good a fit to the morphology of OH231.

Comparison to Key Observations
The outcome of a successful hydro model must replicate the fundamental features of OH231, including the shapes, dimensions, and overall geometry of the spines and bulbs (Figure 1). The speeds of the gas must exhibit matching linear gradients from the core to the lobes' tips with a terminal speed that matches the 12 CO observations. A dense (10 6 cm −3 ) knot-like feature that accounts for the presence of HCO + -a tracer of high-density cold gas-must lie at the leading edge of the northern lobe but not at the southern.
The spray-formed bulbs. The leading edge of ambient gas that is displaced by the spray forms the visible, shocked outer edges of the bulb. The leading edge of this surface quickly displaces about its own mass of ambient gas and slows down. As a result, the clump eventually overtakes the bulb edge and forms polar protrusions in each lobe.
After 800 yr, nearly all of the mass in the original spray has been incorporated into the high-latitude walls of the bulb (Figure 2). At midlatitudes, a fraction of the initial spray still resides in a thin outer crescent. The lobes become bulbous since the walls at these lower latitudes are slower and contain relatively less mass than those at higher latitudes. Moreover, the streamlines strike the reverse shock at midlatitudes obliquely, so the shock speed is about half as large. For these reasons, the walls at midlatitudes are barely visible in shockexcited tracers. Similarly, the Doppler speeds of the mass within the walls will decrease with latitude, as observed. Other regions inside the lobes and outside the spine are essentially empty.
Note that an alternative way to form a bulb is by thermal rather than ram pressure. We explored the injection of a hot, thermally over-pressured low-mass wind (10 6 K) of injection density 10 5 cm −3 and small flow speed <100 km s −1 . However, all such hot-wind simulations failed to recreate the smooth bulb edges. The thermal overpressure behind the walls creates unstable, and persistently and strongly serrated walls as soon as the hot gas is suddenly released into the stationary ambient gas.
The clump-formed spines. Each clump has an initial density of 10 6 cm −3 , a radius of 2000 au, and a mass of 0.1 M e . The clumps are launched into the north (south) lobe at a speed of 225 km s −1 (400 km s −1 ), as implied by the deprojected  Table 1. The inset shows the trail of gas being shed into the spine as the clump is crushed. It applies also to the southern clump at t≈200 yr. Note that the color bars for the speed panels differ for the north and south lobes. Panel B: Like panel A, but for the density distribution of a simulation where all of the initial densities (and the total mass) are reduced by a factor of 10. It demonstrates that the outcomes of isoflow models can have similar large-scale outcomes.
Doppler speeds of the 12 CO at the flow tips. The combined momenta of the clumps are 60 M e km s −1 , or twice the uncertain estimate of A+01.
We emphasize that the clumps are identical aside from their launch speeds. Although the clumps slow down slightly in transit, their terminal speeds remain on the order of their launch speeds as they flow through the co-moving streamlines of the spray. At first, the clumps only gently shock the almost comoving gas through which they flow. By t=800 yr, each clump has overtaken the leading edge of its bulb so the speed of the leading shock is that of the clump or its remnants.
At t=800 yr, the remnant of the southern clump consists of a thin crescent-shaped leading rim of compressed gas with density n cl ≈4×10 6 cm −3 . Within its trailing spinal column, the speed of the gas rises linearly from zero at the clumpʼs launch origin to 400 km s −1 at the flow tip at all times (until the clump is crushed shortly after t=800 yr). The spineʼs density, n spine , rises from 50 cm −3 at y=5000 au to 3000 cm −3 behind the remnant of the clump. The trailing half of the northern clump is still intact at t=800 yr. The thin rim of compressed clump gas at its leading edge has a density n cl ≈4×10 6 cm −3 .
A clump is crushed when it transfers enough momentum to the ambient gas that it significantly slows down and, in doing, so, initiates internal supersonic shear and turbulence. The crushing time t crush =(2r 0 /v cl )·(n cl /n amb ) 1/2 is presumed to be an ambient medium of constant density. In our best simulations, neither the south nor the north clump is crushed in 800 yr, though the south clump is on the verge of crushing when the simulation ends. Note that the onset of crushing of the clumps-especially the northern one-is appreciably delayed since it spends most of its time immersed in the co-moving spray. Note also that the crushing time is the same for all isoflows.
The spines develop immediately in the stagnation column behind the clump. The dominant important source of gas in the spine is gas that was stripped from the outer edges of the clump in transit and falls behind (inset, Figure 2). The speed gradient within the spine develops very quickly and persists through 800 yr (before the clump is crushed).
Thin layers of hot, shocked ambient gas precede the leading edge of both clumps' remnants into the ambient gas after the clumps penetrate their respective bubble edges. This rim hot, shocked gas expands thermally and, for the most part, laterally. It is left behind as the remnants of the clump move onward. Appropriate emission lines will be excited behind the attendant shocks. Although ambient gas does not flow through a stationary contact discontinuity ("CD"), the CD moves steadily into the ambient medium. Under these conditions, slow ambient and fast injected gas will mix to produce a largescale speed gradient in the walls.
At the end of the simulation, the dense gas at the tip of the northern clump can be seen in high-density tracers such as HCO + (see Section 2). The southern clump is almost crushed and its remaining mass and average density are dropping. It might not be detectable in the isoflow with 10 times lower density mentioned above.
Skirts. The simulations fail to reproduce the skirts in the southern lobe that are readily seen in all broadband HST images of OH231. We speculate that the skirts consist of forwardscattered light from dust in the bulb that lie within a stellar illumination cone of opening angle ≈±20°.

Parameter Ranges and Model Robustness
The model discussed above is the "best" of about five dozen models. In the process of exploring the values of the initial conditions, we were able to understand how changes in the starting parameters affect the final outcomes. Here we summarize our experiences.
In our parameter variation studies, we almost always adopted the clump flow speeds and the static ambient density distribution and varied other parameters. We find that changes in n amb (r o ) by factors of about three (with no other changes) clearly had an impact on the evolution of the spray since the spray displaces a considerable volume of static ambient gas and slows down. In short, for a fixed speed of the spray, bulbs of the appropriate dimensions are formed only when that ratio of n spr (r o ) and n amb (r o ) is within 10% of 35. This is a robust result that applies to both lobes.
The outcomes are similar for nearly all values of the clump to spray density ratio so long as the clump is at least as dense as the spray at launch. (However, if n cl (r o )/n spr (r o )3, then the clump or its crushed remains cannot overtake the outer edge of the expanding bulb within 800 yr.) Given the adopted ambient medium, the launch speeds of both flows could not be varied by more than a few percent without a very clear impact on the lengths of the lobes. The requirement that the speed gradients are the same in the spines of both lobes very firmly requires that the launch speeds of the clumps must be in exact proportion to their final offsets from the nucleus; in this case, ∼2:1. The opening angle of the sprays determines the aspect ratios of the bulbs at 800 yr. Departures from our adopted values by 10°or more have noticeable effects on the final dimensions of the bulbs. The same applies to the polar width of the Gaussian taper that we applied to the density and speed of the conical flow to give the spray its overall form.
The large initial radius of the clumps that we have adopted, 2000 au, is obviously an issue given that the dimensions of a binary star system or the dominance of a strong surface magnetic field, 100 μG, will be far below 100 au. We have adopted this clump radius for a practical reason-the clump results in an axial protrusion of the proper total width, ∼7000 au-and an expedient one-we easily resolve it on the model grid. We speculate that, with suitably adjusted initial conditions, the "magnetic bomb" model of Huarte-Espinosa et al. (2012) might produce clumps of the appropriate size and speed within the first hundred years after launch. In any event, readers are cautioned that the clump might simply be a region born with complex internal structure and dynamics.

Comparison with the RRA05Model
RRA05 approached the hydrodynamical modeling of OH231 differently than we did. They launched columnar pulses of various properties from the central star into asymmetric ambient media in order to account for the geometric properties of the lobes. The dense axial cores of the pulses have a radius of 3300 au (about the same as the size of the clump in our bestfit simulations) and a core pulse density of 10 4 cm −3 (100 times sparser than the best-fit clump in our simulations). The ejection speed within each pulse decreases linearly with distance from the source from 200 km s −1 at their leading edge, so the Hubble-like flow pattern that develops as each pulse pushes forward is forced by the initial conditions. (In our simulations, this gradient is a natural, robust, and persistent consequence of the clump−AGB-wind interaction during the early light-flow phase.) One pulse flows through the winds of an AGB star located near the end of the northern lobe and is offset by various amounts from the flow axis. (NB: No AGB star or evidence of a reddened disk appears in about the proper location in an HST image at 1.6 μm.) The AGB wind has an equatorial density contrast ξ that varies from 1 to 50 in their five simulations, and its equator is orthogonal to the pulse flow axis. In models with ξ25, the pulse is impeded as it travels through the edge of the wind disk and slows down. (In effect, the disk of the offset AGB star serves as a sheet-like barrier through which the pulse penetrates.) The lobes created by pulses of disparate initial characteristics successfully develop the observed ratio of the lengths of both lobes and reproduce their speed gradients. However, the ratio of the Doppler shifts at the end of the lobes is close to unity and the lobes have no spines or skirts.

Summary and Conclusions
We have constructed a successful hydrodynamical simulation of the lobes of OH231 using a combination of pairs of clumps that are embedded within sprays of brief duration and launched into static ancient AGB winds. (This is not unlike the ejection of a cork and a spray of high-pressure fluid when a bottle of champagne is opened, except that in our case the spray is ejected shortly before the cork.) After 800 yr, each of the clumps leaves behind a tail, or spine, in which the steady kinematic gradient seen in the Doppler shifts of 12 CO is replicated. The brief spray with an opening angle of±50°e xpands to form the wide southern bulbs. Its northern counterpart is smaller, denser, and thinner since it is ejected at about half the speed and opening angle.
The values of the initial conditions adopted for our best simulation are almost fully constrained with a wide variety of observations of the structure, kinematics, and molecular masses of OH231 (Figure 1). In addition to matching the shapes of the lobes, one feature of the models is especially conspicuous: the linear speed gradients along the outflows. The injected gas is characterized by a single speed in each lobe of OH231 and lies at the lobe perimeter as the simulation evolves. So, the speed gradient is not imposed by initial conditions. It is a natural and robust result of a combination of swept-up ambient gas and clump material that has been shed from the edges of the clumps. Steady jets or sprays of any geometry are incapable of explaining these speed gradients. In addition, they would not replicate the morphology of the bulbs while also producing the highly conspicuous spines.
We also showed that the same outcomes of morphologies and kinematics are expected over a range of injected masses, provided that the ratios of the densities of the ambient gas, the clumps, and the spray are fixed. That is, members of the family of simulations with these ratios but with different masses and momenta are compatible with changes in the total masses of the injected gas within factors of about 10 of the highly uncertain mass and momentum estimates for 12 CO in OH231 (e.g., A +01, B+02).
The initial conditions that we used for the simulations ( Table 1) provide important glimpses into all but perhaps the first century of the nebulaʼs evolution. The similar initial densities and kinematics of the pairs of lobes and sprays suggest that they were formed in a common event. This result is a natural outcome; it was never forced on the models (except as a first guess for the initial conditions at the start of this program). The shaping of the present structures of OH231 occurred at early times when the density of the flows was less than or comparable to that of the ambient gas. Once the ensemble of flows emerge into the outer, low-density zones of the ambient AGB winds, they expand ballistically, as expected for "heavy" or "fast-propagating" flows (Soker 2002).
While the construction of a successful hydro model for OH231 is of interest, the lasting value of this work is really the constraints that the model imposes on the history of the original mass ejection process. Before going further, let us note that the ejection speeds are fixed by measurements, so flow momenta scale with flow masses. Then, from Table 1, we find that the masses and momenta of the southern spine (bulb) derived from the best-fit model are 0.1 M e and 8×10 39 gm cm s −1 (2.3 M e and 1.8×10 41 gm cm s −1 ), respectively. The corresponding values for the northern lobe are 0.1 M e and 4×10 39 gm cm s −1 (0.5 M e and 0.2×10 40 gm cm s −1 ), respectively. Note that the bulb (which is not directly detectable in 12 CO) dominates these estimates. The total ejected mass and momentum are thus 3 M e and 2.1×10 41 gm cm s −1 , respectively. (In effect, they are anchored to the directly measured value of the mass of the slow flow (i.e., the 12 CO(1-0) flux of the narrow line), about which there is considerable observational uncertainty (Section 2).) These values are exceptionally large among all other pPNe. This paper is the first of several like it, in which hydrodynamic models similar to the one here are developed in order to fit detailed observations of other symmetric PNe. Our primary scientific goal in this series is to characterize the aspects of the histories that pPNe may share in common. To that end, a detailed model for one pPN (and a fairly unique one at that) is simply a first step. Hence, our hydrodynamic studies of pPNe will continue on a case-by-case basis into the future.
B.B. wishes to recognize the 50 years of boundless mentoring and loving friendship of his thesis adviser, Prof. Yervant Terzian.
We are grateful to J. Alcolea for useful early discussions. Nearly all of the hydrodynamic models in this paper were run in computers maintained by the Center for Integrated Research Computing of the University of Rochester. (A few of the very early models used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575.) Many of the images used in this project were made possible by HST GO grant 11580. Support for GO11580 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. Additional data that were used in this paper were obtained from the Multimission Archive (MAST) at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts.
Note added in proof. We inadvertently and unfortunately overlooked and neglected to cite a very relevant paper by Akashi & Soker (2008). The paper includes an alternate hydrodynamic model of OH231. However, none of the principal findings of this paper are affected.