The Astrophysical Journal, 566:392-398, 2002 February 10
© 2002. The American Astronomical Society. All rights reserved. Printed in U.S.A.

 

Runaway Acceleration of Line-driven Winds: The Role of the Outer Boundary

Achim Feldmeier
Astrophysik, Institut für Physik, Universität Potsdam, Am Neuen Palais 10, 14469 Potsdam, Germany; afeld@astro.physik.uni-potsdam.de
Isaac Shlosman 1,2
JILA, University of Colorado, Box 440, Boulder, CO 80309-0440; shlosman@pa.uky.edu
and
Wolf-Rainer Hamann
Astrophysik, Institut für Physik, Universität Potsdam, Am Neuen Palais 10, 14469 Potsdam, Germany; wrh@astro.physik.uni-potsdam.de

Received 2001 May 4; accepted 2001 October 11

ABSTRACT

Observations and theory suggest that line-driven winds from hot stars and luminous accretion disks adopt a unique, critical solution that corresponds to maximum mass-loss rate. We analyze the numerical stability of the infinite family of shallow wind solutions, which resemble solar wind breezes, and their transition to the critical wind. Shallow solutions are subcritical with respect to radiative (or Abbott) waves. These waves can propagate upstream through shallow winds at high speeds. If the waves are not accounted for in the Courant time step, numerical runaway results. The outer boundary condition is equally important for wind stability. Pure outflow conditions, as assumed in the literature, trigger runaway of shallow winds to the critical solution or to accretion flow.

Subject headings: accretion, accretion disks; hydrodynamics; instabilities; stars: mass loss; waves

     1 JILA Visiting Fellow.
     2 Permanent address: Department of Physics and Astronomy, University of Kentucky, Lexington, KY 40506-0055.

1. INTRODUCTION

     Line-driven winds (LDWs) occur in various astronomical objects, such as OB and Wolf-Rayet stars, in accretion disks in cataclysmic variables, and, probably, in active galactic nuclei and luminous young stellar objects. These winds are powered by absorption and re-emission of photospheric continuum flux in numerous spectral transitions of C, N, O, and Fe, among other ions.

     Castor, Abbott, & Klein (1975, hereafter CAK) analyzed the steady state Euler equation for LDWs. They found an infinite family of mathematical solutions, but only one, hereafter the "critical solution," that extends from the photosphere to arbitrarily large radii. Other solutions either do not reach infinity or do not reach the photosphere. The former solutions are called shallow and the latter ones steep. The unique, critical wind starts on the fastest shallow solution and switches smoothly to the slowest steep solution at the critical point.

     Observational support that LDWs adopt the critical solution comes from measured terminal speeds (Abbott 1982). Furthermore, mass-loss rates of supergiant winds are in general agreement with modified CAK theory (Lamers & Leitherer 1993; Puls et al. 1996). These measurements were recently extended to include Galactic and extragalactic OB and A stars and central stars of planetary nebula (Kudritzki et al. 1999).

     Abbott (1980) interpreted CAK theory as a complete analogy to the solar wind and nozzle flows. The existence of a sonic point defines the unique, transonic solutions for these flows, whereas the existence of a critical point for Abbott waves defines the unique, CAK solution for LDWs. Only from below this critical point can Abbott waves propagate upstream toward the photosphere. Above the critical point, they are advected outward. Because Abbott waves generally propagate highly supersonically, the critical point of LDWs lies at much higher speeds than the sonic point.

     Abbott's (1980) analysis was challenged by Owocki & Rybicki (1986), who derived the Green's function for a pure absorption LDW. The Green's function yields correct signal speeds in the presence of hydrodynamic instabilities. The inward signal speed in a pure absorption line wind is the sound speed, and not the much higher Abbott speed, because photons propagate outward only. Owocki & Rybicki (1986) showed that a fiducial upstream signal, which still propagates inward at Abbott speed, must be interpreted as purely local Taylor series reconstruction. For a flow driven by scattering lines, however, Owocki & Puls (1999) find physically relevant Abbott waves for a numerical Green's function.

     In the present paper, we further analyze the properties of Abbott waves. We show that they are crucial for our understanding of stability of LDWs and must be included in the Courant time step. So far, time-dependent numerical simulations of LDWs from stars and accretion disks have ignored the ability of Abbott waves to communicate in the supersonic regime, which results in a numerical runaway. In particular, this runaway can lift the wind to the critical solution.

     The critical solution is also enforced by applying pure outflow boundary conditions. It is often argued that outflow boundary conditions are appropriate because LDWs are highly supersonic. Instead, they have to be superabbottic. We show that shallow wind solutions, which correspond to solar wind breezes, are everywhere subabbottic. Hence, these solutions are numerically destabilized by applying outflow boundary conditions.

     We formulate boundary conditions that render shallow solutions numerically stable. Those include nonreflecting Riemann conditions for Abbott waves. By allowing for kinks in the velocity law, shallow solutions can be made globally admissible.

2. THE WIND MODEL

     In the CAK model for LDWs, both gravity and line force scale with r-2. If the sound speed and, hence, the pressure forces are set to zero, this leads to a degeneracy of the critical-point condition, which is satisfied formally at every radius (Poe, Owocki, & Castor 1990). Thus, in this case, Abbott waves cannot propagate inward from any location in the wind. For finite sound speed, they creep inward at low speed. Inclusion of the finite disk correction factor is much more relevant for LDWs than inclusion of pressure forces. With the finite disk included, the inward speed of Abbott waves below the critical point is significantly higher than the wind speed. Unfortunately, the finite-disk correction factor depends on the (unknown) velocity law of the wind, preventing a simple analysis of the wind dynamics.

     We consider, therefore, a wind model that is analytically feasible and yet prevents the (near-)degeneracy of the CAK point star wind. (Especially, this degeneracy leads to poor convergence of time-dependent numerical schemes.) As a prototype, a vertical LDW from an isothermal, geometrically thin, non-self–gravitating accretion disk is assumed. The sound speed is set to zero. Keplerian rotation is assumed within the disk and angular momentum conservation above the disk. This reduces the flow problem to a one-dimensional, planar one. The radiative flux above an isothermal disk is roughly constant at small heights. On the other hand, the vertical gravity component along the wind cylinder is zero in the disk midplane, grows linearly with z if z ≪ R (with R the footpoint radius in the disk), reaches a maximum, and drops off at large z. To model the launch region of the wind and the gravity maximum, we choose g(z) = z/(1 + z2), with normalization GM = 1 and R = 1, where G is gravitational constant and M is the mass of the central object. The differing spatial dependences of flux and gravity result in a well-defined critical point in the flow.

     For constant radiative flux, the CAK line force becomes gl = C(ρ-1dv/dz)α, where ρ and v are the density and velocity and α and C are constants. We choose α = 1/2. Typical values of α from non–LTE calculations (Pauldrach 1987) and Kramers' law (Puls, Springmann, & Lennon 2000) are α = 0.4 to 0.7. The constant C can be expressed in terms of the unique, critical CAK mass-loss rate. This is the maximum allowed mass-loss rate for a stationary wind. We introduce new fluid variables m and w,



where subscript c refers to the critical CAK solution. For stationary, plane-parallel winds, the continuity equation becomes ρv = const, and m is a constant. In the following, the quantity w′ = dw/dz = v dv/dz will play a central role. The velocity law, v(z), is obtained by quadrature of the wind "solution" w′(z). For α = 1/2, the stationary Euler equation becomes



where we multiplied the nominator and denominator in the line force by v, and introduced gc = g(zc), the gravity at the location of the CAK critical point. The dependence of w′ on z is suppressed in the equation. The constant C is determined as follows. Equation (2) is a quadratic equation in . Consider the E-plane. For arbitrary z and sufficiently small m, the crossings E = 0 of the parabola E() with the abscissa determine two wind solutions w′(z), termed shallow and steep. When m increases, the crossings approach each other, and merge at some maximum m ≥ 1, beyond which no further solutions exist. Hence, the singularity condition ∂E/∂w′ = 0 holds when two solutions merge, as is the case at the critical point. Together with E = 0, this yields w′c = gc and C = 2(gcρcvc)1/2.

     We add a comment that is important for what follows. At every location z ≠ zc, m > 1 for the parabola that meets the abscissa in just one point. The minimum, m = 1, is reached at the unique critical point, zc, which is the flow nozzle and determines the maximum mass loss that can be driven from the photosphere to infinity. Solutions (actually, pairs of shallow and steep solutions) with m > 1 are called overloaded. They correspond to a choked Laval nozzle flow. These solutions become imaginary in the neighborhood of the critical point. We have shown elsewhere (Feldmeier & Shlosman 2000) how overloaded LDW solutions can be defined globally, by letting them decelerate, dv/dz < 0, around the critical point.

     The position of the critical point is determined by the regularity condition, which is derived next. Because E = 0 everywhere, dE/dz = 0 holds, too. Because ∂E/∂w = 0, it follows that ∂E/∂zc = 0, if |w′′c| < ∞. This latter condition singles out the unique critical point, because overloaded winds have kinks at the merging point of two solutions, |w′′| ≡ ∞. Hence, the critical point lies at the gravity maximum, dg/dzc = 0. For g(z) = z/(1 + z2), zc = 1 and gc = 1/2. The line force in equation (2) becomes (2w′/m)1/2. This holds for any solution, whether shallow or steep. Solving equation (2) for ,



where we introduced a new variable, q. The signs refer to steep (plus) and shallow (minus) solutions. For m = 1, the square root in equation (3) vanishes at the critical point. For m > 1, the root becomes imaginary in the neighborhood of the critical point. The solutions w′(z) of equation (3) have a saddle topology in the zw′-plane when m is varied from m < 1 to m > 1. The critical point—zc = 1, w′c = 1/2—is the saddle point. A variable radiative flux above a nonisothermal accretion disk leads instead to two different values of zc, above and below the gravity maximum (Feldmeier & Shlosman 1999). The lower is a saddle, the upper an extremum.

     Note that the CAK critical point of LDWs is a saddle in the z-vv′ plane, whereas the sonic point of the solar wind is a saddle in the zv-plane. In all other respects, we find a deep analogy between LDWs and the solar wind. Only the notion of a critical speed is replaced by a critical acceleration. Shallow solutions are the analogs to solar wind breezes.

     According to equation (3), shallow and steep solutions with m ≤ 1 extend from z = 0 to ∞ for the present model with zero sound speed. CAK found that steep solutions start supersonically at the wind base z = 0, which seems unphysical. In a spherically symmetric wind, shallow solutions do not extend to arbitrarily large z, because they cannot provide the required expansion work. These two results and the requirement for a continuous and differentiable solution imply that the wind has to pass through the critical point. It starts on the fastest shallow solution and ends on the slowest steep one.

     Figure 1 (which is adapted from Fig. 4 in Abbott 1980 and Fig. 8.13 in Lamers & Cassinelli 1999) shows the LDW solution topology in the rv-plane for a finite sound speed. In region II (regions defined in CAK), which is spatially the most extended in real LDWs, shallow and steep solutions exist at each point (r, v). If z ≠ zc, these solution pairs may be overloaded. In the subsonic region I, only shallow solutions (possibly overloaded) exist; steep solutions are necessarily supersonic. Correspondingly, only steep solutions (possibly overloaded) exist in region III. No solutions exist in regions IV and V. Note that the gap in overloaded solutions around zc can be bridged by an extended decelerating region, dv/dr < 0.


Fig. 1   Solution topology of line-driven winds, showing shallow, steep, and overloaded solutions. The circle marks the CAK critical point. Roman numerals refer to cases I–V of CAK. Regions I–III and V are mutually separated by straight lines (sound speed v = a; Parker point r = rP); regions II and IV are separated by the single locus.

     The argument by CAK ruling out global, steep solutions is stringent, but the exclusion of shallow solutions seems too restrictive. The pressure force due to spherical expansion becomes important only at very large distances from the wind base. According to CAK, shallow solutions break down around 300 stellar radii for realistic O star wind parameters. Most shallow solutions have by then overcome the local escape speed. Shallow solutions can be globally defined by a slight generalization of the CAK approach, allowing for negative v′ = dv/dz in the line force ∼. This can be done by replacing v′ with |v′| or with max(v′, 0). The former expresses the "blindness" of purely local line driving to flow acceleration or deceleration. The latter accounts for the nonlocal effect of shadowing in nonmonotonic velocity laws, creating a resonance coupling. All radiation is assumed to be absorbed at the first resonance location. Realistically, the value of the line force lies between these extremes. The decelerating wind solution, w′ < 0, is unique and does not split into shallow and steep branches.

     Hyperbolic differential equations allow for weak discontinuities, which are represented in the zw′-plane by jumps between different solution branches. Hence, a shallow solution may jump onto the decelerating branch. Even an overloaded solution can avoid becoming imaginary by jumping onto the decelerating branch in the vicinity of the critical point. This opens the possibility that global solutions can be composed in a piecewise manner.

     Jumps in w′ introduce kinks in v(z). Such kinks have been found on various occasions in time-dependent LDW simulations. They may play a central role in the observed discrete absorption components in nonsaturated P Cygni line profiles (Cranmer & Owocki 1996). In LDWs from accretion disks in cataclysmic variables, a jump onto the decelerating branch even occurs for the critical solution, which would elsewhere not extend to infinity (Feldmeier & Shlosman 1999). Hence, the smoothness argument of CAK cannot be used to disqualify shallow solutions.

3. ABBOTT WAVES

3.1. Linear Dispersion Analysis

     The time-dependent continuity and Euler equations for the present wind model are







In the Euler equation, v′ > 0 is assumed. The case v′ < 0 will be discussed below. To derive Abbott waves, small perturbations are applied, v = v0 + v1exp[i(kz - ωt)] and ρ = ρ0 + ρ1exp [i(kz - ωt)]. The subscript 0 refers to a stationary solution, subscript 1 to the perturbation. Linearizing yields







with nondimensional χ ≡ k v0/v and q0 = 1/2. Setting the determinant of the system to zero yields the dispersion relation ω(k). In a WKB approximation, χ ≫ 1, one finds for the phase speed A0 and the growth rate of the inward (minus) mode,



in the observer's frame, and for the outward (plus) mode,



with -ε a small damping term of no further consequence. Along a shallow solution q0 < 1, and along a steep solution q0 > 1. For shallow solutions with small m0 and w, q0 is small, and the inward mode can propagate at arbitrarily large phase speed from every location z. This result is new, because so far Abbott waves along shallow solutions were not considered ("a LDW has no analog to a solar breeze"; Abbott 1980). Because for the CAK point star model every point is critical, q0 = 1, and the inward mode is stagnating everywhere. Stability of non–WKB waves, χ = O(1), is found from a numerical solution.

     Integrating equation (3), the velocity of the critical solution at zc = 1 (the "critical speed") is v ≈ 0.61. On the other hand, a shallow solution with m0 = 0.8 has v > 0.61 for z > 1.51. This is no contradiction to the fact that the shallow solution is subabbottic.

3.2. Characteristic Analysis

     The above analysis holds for linear waves, and the wave speed is determined by the underlying stationary solution. General, nonlinear waves can be derived from a characteristic analysis. The characteristic form of the above equations is, without further approximations,







(Feldmeier & Shlosman 2000), where the Abbott speed, A, in the observer's frame is



     The interpretation of equation (11) is unconventional, in that v′ should now be considered the fundamental hydrodynamic field, instead of v. This happens because of the nonlinear dependence of E on v′ (Courant & Hilbert 1968).

     Hence, -ρv′ in the continuity equation and -g′/ρ in the Euler equation are inhomogenous terms, because they do not contain derivatives of the hydrodynamic fields ρ and v′. The advection operators are ∂/∂t + v∂/∂z and ∂/∂t + A∂/∂z, with characteristic speeds v and A. The former should be read as v + a in zero sound speed approximation. The Riemann invariants ρ and v′/ρ correspond to wave amplitudes. In a frame moving at v, the amplitude ρ is constant, except for sources and sinks -ρv′. In an isothermal gas p ∼ ρ, and the wave amplitude is indeed that of a pressure wave. In a frame moving at A, the amplitude v′/ρ is constant, except for sources and sinks -g′/ρ. The Sobolev optical depth is ∼(v′/ρ)-1, indicating that the inward mode is a radiative wave.

     Even for waves of infinitesimal amplitude, v′ may deviate strongly from v, if the wavelength λ is short. Because A ∼ (v′)-1/2, short-scale Abbott waves are highly dispersive. For sufficiently small λ, v′ becomes negative at certain wave phases. The characteristic equations (10) and (11) hold also for v′ < 0, but A has different meaning then. For a line force ∼|v′|1/2, the general A is



where the lower sign applies for v′ < 0 and in q- one has to use -v′. Trivially, Abbott waves are advected outward in superabbottic flow. But that Abbott waves propagate downstream as an outward mode in regions where v′ < 0 is probably their single most unique property.

     For a line force ∼[max(0, v′) ]1/2, then, A = v if v′ < 0. This is because the line force drops out of the Euler equation, and both wave modes become ordinary sound. (Actually, A = v - a for upstream-propagating sound.) In this case one has, along a characteristic curve with coordinate s,



as the WKB wave solution of equations (10) and (11). This describes gasdynamical simple waves. Because v′(0) < 0, they steepen and break (v′ = -∞). The derivation of simple Abbott waves is beyond the scope of this paper.

4. BOUNDARY CONDITIONS FOR LINE-DRIVEN WINDS

     We discuss now the outer boundary conditions for shallow solutions. Abbott waves can enter through the outer boundary; hence, pure outflow boundary conditions are not appropriate for shallow solutions.

     As hydrocode we use a standard Eulerian scheme on a staggered mesh. Advection is solved in the conservative control volume approach, using monotonic van Leer (1977) interpolation. The sound speed is set exactly to zero, and a small amount of artificial viscosity is added to handle occasional shock fronts. Details of the scheme can be found in Reile & Gehren (1991) and Feldmeier (1995).

     In a first step, we consider simple boundary conditions for ρ and v that still preserve essential features of Abbott waves. These boundary conditions render shallow solutions stable and allow small perturbations to leave the grid. They are



where 0 and I are the mesh indices of the inner and outer boundary, respectively. The conditions (eq. [15]) are motivated as follows. (The reader may jump to the end of this paragraph.) On a characteristic curve [z(s), t(s)] that leaves the mesh through a boundary, the characteristic equation is solved on the boundary via extrapolation, using one-sided, interior differentials. Such an upwinded scheme is stable (Steger & Warming 1981). For shallow winds, the ρ wave leaves through the outer boundary if v > 0; hence, ρ is extrapolated from the interior in (15). Zeroth-order extrapolation is used for stability reasons. Abbott waves enter the mesh through the outer boundary; hence, a boundary condition must be applied. We choose v(I) = const. [To prevent standing waves, v(I) = m0/ρ(I) is often more appropriate, but is also more susceptible to runaway than v(I) = const.] At the inner boundary, the argument proceeds correspondingly, with interchanged roles of the waves. We find that using ρ(0) = ρ(1) and v(0) = const instead destabilizes shallow solutions.

4.1. The Courant Time Step

     In time-dependent hydrodynamic simulations published so far, it is customary to insert the sound speed in the Courant time step,



with Courant number σ ≤ 1, and the minimum has to be taken over the mesh. We see now that Abbott waves, as the characteristic inward wave mode, have to be included in the Courant time step,



where A is in the fluid frame now. Note that A changes sign with v′. For sufficiently small v′ > 0 and for arbitrary v′ < 0, Abbott waves determine the time step. At velocity plateaus where the wind velocity is more or less constant, A→-∞, from equation (11), and the Courant time step →0. This is an artifact of Sobolev approximation, in which the line radiation force tends to zero in the absence of velocity gradients. In practice, in calculating the Courant time step we do not allow q to drop below a certain minimum value, usually 10-4. This corresponds to σ ≈ 10-4 for an ordinary Courant time step, not including Abbott waves. Our results are not sensitive to the value of this minimum q, if the latter is sufficiently small. Furthermore, velocity plateaus quickly acquire a tilt and propagate at finite speed.

4.2. Nonreflecting Boundary Conditions

     We can then formulate nonreflecting boundary conditions for the Riemann invariants ρ and v′/ρ, which annihilate incoming waves. This prevents boundary reflection of waves that originate on the interior mesh. To annihilate a wave, its amplitude (the Riemann invariant) is kept constant in time at the boundaries, ∂(v′/ρ)/∂t = 0 and ∂ρ/∂t = 0 for Abbott and sound waves, respectively (Hedstrom 1979; Thompson 1987). One may say that nonreflecting boundary conditions drive the numerical solution toward a neighboring, stationary solution without waves. The technical details of implementing these boundary conditions are discussed in Appendix A.

5. STABILITY OF SOLUTIONS

     We discuss here the numerical stability of shallow and steep solutions, depending on the appropriate Courant time step and boundary conditions.

5.1. Stability of Steep Solutions

     We start with steep solutions and show that they are unconditionally stable, even when simplified boundary conditions analogous to equation (15) are used. According to equation (3), q > 1 everywhere for steep solutions, and Abbott waves can propagate only toward larger z. At the outer boundary, the extrapolations ρ(I) = ρ(I-1) and v(I) = v(I-1) are therefore appropriate. At the inner boundary, two conditions have to be specified. The obvious choices are either (ρ = const, v = const) or (ρ = const, v′ = const). The first set fixes the mass-loss rate, whereas the second establishes nonreflecting boundary conditions by keeping the wave amplitudes ρ and v′/ρ constant in time, which means that there are no waves.

     We find that steep solutions are stable for both types of boundary conditions. Even when the initial conditions differ profoundly from those of a steep wind, the numerical code converges to a steep solution. Furthermore, we performed tests with an explicit, harmonic perturbation at a fixed Eulerian position in the wind. Even if the perturbation amplitude reached 100%, the wind remained on average on a steep solution, which can therefore be considered as unconditionally stable.

5.2. Stability and Runaway of Shallow Solutions

     A physically more relevant question is the stability of shallow solutions. In a first step, we use an analytic, shallow wind as initial conditions. Table 1 shows the runaway of this shallow solution to the critical velocity law. Abbott waves are not accounted for in the Courant time step, and pure outflow boundary conditions, ρ(I) = ρ(I - 1) and v(I) = v(I - 1), are applied. This corresponds to the typical procedure adopted in the literature. The runaway starts at the outer boundary and generates inward-propagating Abbott waves. The velocity law evolves toward higher speeds, until the critical solution is reached.

Table 1   Numerical Runaway of a Shallow Velocity Law

     The mechanism of the runaway seems to be as follows. The Courant condition for Abbott waves is violated in the outer wind, where the velocity law is flat and the Abbott speed is high. This causes numerical runaway. In contrast to standard gas dynamics, the LDW runaway does not crash the simulation. We speculate that exponential growth in v creates high |v′|, which implies low Abbott speeds. The Courant condition is no longer violated, and runaway stops. Abbott waves are also excited at the outer boundary, by switching from the interior to the boundary scheme. Upstream-propagating Abbott waves are inconsistent with the assumed outflow boundary conditions. The waves propagate to the wind base and drive the inner wind toward the critical solution, which is indeed consistent with outflow boundary conditions.

     Note in Table 1 that Abbott waves can propagate inward from the outer boundary, zM = 4, until the critical velocity law is reached (indicated by bold numbers in the table). This is no contradiction to the critical point being located at zc = 1: the outer, shallow portion of the velocity law (the numbers in italics) is subabbottic, even when the inner, steep velocity law is already superabbottic.

     To prevent numerical runaway of shallow solutions, we apply nonreflecting Riemann boundary conditions and include Abbott waves in the Courant time step. As initial conditions, an arbitrary (for example, linear) velocity law is used. Abbott waves are again excited at the outer boundary and propagate inward. However, the simulation relaxes quickly to a shallow solution. The m-value depends on the initial data. Even for an initial model with height-dependent mass flux, we find quick convergence to a shallow solution.

     Further details on the stability of shallow solutions are given in Appendix B. To summarize this section, both improper outer boundary conditions and Courant conditions that do not account for Abbott waves are responsible for numerical runaway.

6. SUMMARY

     We present a simplified model for line-driven winds from stars and accretion disks that avoids the r-2 degeneracy of the CAK model. Radiative, or Abbott, waves are derived from both a dispersion and a characteristic analysis and for all possible solutions to the Euler equation, i.e., shallow, critical, steep, and overloaded ones.

     We find that Abbott waves can propagate upstream, toward the photosphere, from any position along a shallow wind solution. Hence, for shallow winds, Abbott waves can enter the calculational grid through the outer boundary. The wave propagation speed depends on the velocity slope of the wind. Abbott's (1980) result that these waves are creeping, i.e., have only slightly negative inward speeds, holds globally only for the almost degenerate CAK point star model. In more realistic wind models, Abbott waves can limit the Courant time step.

     The neglect of Abbott waves in either the Courant time step or the boundary conditions leads to numerical runaway toward the critical wind solution of maximum mass-loss rate or to accretion flow. If, instead, incoming waves are annihilated at the outer boundary and the correct Courant time step is used, shallow solutions are stable.

     It is a pleasure to thank Janet Drew, Leon Lucy, Stan Owocki, and Joachim Puls for frequent discussions. This work was supported in part by NASA grants NAG 5-10823, NAG 5-3841, WKU-522762-98-6, and HST GO-08123.01-97A to I. S., which are gratefully acknowledged.

APPENDIX A

NONREFLECTING BOUNDARY CONDITIONS

     We adopt the following procedure to calculate ρ and v on a boundary. (1) If the ρ wave leaves the grid, = -(ρv)′ is calculated using one-sided, interior differentials. (2) If the ρ wave enters the grid, = 0 is set (wave annihilation; see § 4.2). (3) If the v′/ρ wave leaves the grid, ′/v′ is calculated from the Euler equation (eq. [11]) in characteristic form, using one-sided, interior differentials. Note that v′′ appears here. (4) If the v′/ρ wave enters the grid, ′/v′ = /ρ is set, assuming WKB waves. (5) The new values for ρ and v on the boundary are found by integrating and ′ using a time-explicit scheme.

     In step 5, one actually needs , not ′. To integrate from time step n - 1 to n at the outer boundary using ′, we write



Using v′n = (v - v)/Δz, and similarly for v′n-1, this becomes



whence v depends on v. Because of frequent variable updating in an operator splitting scheme during the time step, v would have to be stored as an extra variable. This seems undesirable, and indeed it causes numerical runaway of shallow solutions. Fortunately, a simple approximation to equation (A2) yields satisfactory results, namely, setting v ≡ v, or



APPENDIX B

NUMERICAL STABILITY OF SHALLOW SOLUTIONS

     For numerical stability tests, we assume the following parameters, if not otherwise stated: 200 equidistant grid points from z = 0.1 to z = 4. Below z = 0.1, the stationary wind speed is very small. Because δv/A0 = δρ/ρ 0, one has v0 ≪ A0 for an inner boundary zm→0. Even small perturbations, δρ, then cause negative speeds at zm. This alters the direction of the ρ characteristic and often causes numerical problems. A shallow solution with m = 0.8 serves as start model. The sound speed is set to zero. The flow time from z = 0.1 to z = 4 is 12.4 for a shallow solution with m = 0.8, and 9.1 for the critical solution (in the units specified in the main text). For m = 0.8, linear Abbott waves propagate from z = 4 to z = 0.1 in a time of 6.0. All simulations are evolved to time of 50. Table 2 lists results from simulations using different outer boundary conditions and Courant time steps. On the inner boundary, conditions (eq. [15]) were used.

Table 2   Runaway in Simulations Applying Different Boundary Conditions and Courant Time Steps

     Consider first the line force ∼|v′|1/2. The table shows that inclusion of Abbott waves in the Courant time is mandatory to prevent runaway (runs 1–5; "shallow" in the column for the outer boundary conditions refers to eq. [15]). For small Courant numbers, Abbott waves are, to a degree, already accounted for by an ordinary time step. Runaway may then be prevented, as in run 2, though the solution does oscillate at all times. If the outer boundary is shifted outward, the runaway occurs again (run 5). If Abbott waves are accounted for in the time step, outflow boundary conditions (extrapolation of ρ and v) do not necessarily lead to runaway. Instead, a stable, shallow solution may be maintained (runs 6 and 9; also run 2). However, the solution oscillates at all times. Shallow solutions become more unstable when the outer boundary is shifted outward, into the shallow part of the velocity curve where Abbott waves are easily excited. Finally, the control run 10 shows that the initial shallow solution is numerically stable if Abbott waves are included in the time step and boundary conditions of equation (15) are used.

     For a line force ∼ 1/2, the scheme is even more susceptible to numerical runaway. All cases with outflow boundary conditions undergo runaway. However, for the boundary condition v(I) = v(I - 1), the solution no longer tends toward the critical solution. Instead, v(I) drops to zero and with it the interior wind velocity. The outcome of the simulation is not clear, because we have not formulated appropriate line-driven accretion boundary conditions. The reason for the drop in v(I) is that an accidental v′ < 0 at the outer boundary implies zero line force. This causes further velocity drop, and v′ becomes more negative. Abbott waves carry this wind breakdown to the interior mesh. It is not clear whether this type of line force should actually be applied in the proximity of the outer boundary.

REFERENCES