Hydrodynamics of astrophysical winds

driven by scattering in spectral lines

Achim Feldmeier

Habilitation thesis in Astrophysics

Potsdam University, June 2001

Meinen Eltern

Instead of a quote

Abstract. Line driven winds are accelerated by the momentum transfer from photons to a plasma, by absorption and scattering in numerous spectral lines. Line driving is most efficient for ultraviolet radiation, and at plasma temperatures from 10⁴ K to 10⁵ K. Astronomical objects which show line driven winds include stars of spectral type O, B, and A, Wolf-Rayet stars, and accretion disks over a wide range of scales, from disks in young stellar objects and cataclysmic variables to quasar disks. It is not yet possible to solve the full wind problem numerically, and treat the combined hydrodynamics, radiative transfer, and statistical equilibrium of these flows. The emphasis in the present writing is on wind hydrodynamics, with severe simplifications in the other two areas. I consider three topics in some detail, for reasons of personal involvement. 1. Wind instability, as caused by Doppler de-shadowing of gas parcels. The instability causes the wind gas to be compressed into dense shells enclosed by strong shocks. Fast clouds occur in the space between shells, and collide with the latter. This leads to X-ray flashes which may explain the observed X-ray emission from hot stars. 2. Wind runaway, as caused by a new type of radiative waves. The runaway may explain why observed line driven winds adopt fast, critical solutions instead of shallow (or breeze) solutions. Under certain conditions the wind settles on overloaded solutions, which show a broad deceleration region and kinks in their velocity law. 3. Magnetized winds, as launched from accretion disks around stars or in active galactic nuclei. Line driving is assisted by centrifugal forces along co-rotating poloidal magnetic field lines, and by Lorentz forces due to toroidal field gradients. A vortex sheet starting at the inner disk rim can lead to highly enhanced mass loss rates.

Zusammenfassung. Liniengetriebene Winde werden durch Impulsübertrag von Photonen auf ein Plasma bei Absorption oder Streuung in zahlreichen Spektrallinien beschleunigt. Dieser Prozess ist besonders effizient für ultraviolette Strahlung und Plasmatemperaturen zwischen 10⁴ K und 10⁵ K. Zu den astronomischen Objekten mit liniengetriebenen Winden gehören Sterne der Spektraltypen O, B und A, Wolf-Rayet-Sterne sowie Akkretionsscheiben verschiedenster Größenordnung, von Scheiben um junge Sterne und in kataklysmischen Veränderlichen bis zu Quasarscheiben. Es ist bislang nicht möglich, das vollständige Windproblem numerisch zu lösen, also die Hydrodynamik, den Strahlungstransport und das statistische Gleichgewicht dieser Strömungen gleichzeitig zu behandeln. Die Betonung liegt in dieser Arbeit auf der Windhydrodynamik, mit starken Vereinfachungen in den beiden anderen Gebieten. Wegen persönlicher Beteiligung betrachte ich drei Themen im Detail. 1. Windinstabilität durch Dopplerde-shadowing des Gases. Die Instabilität bewirkt, dass Windgas in dichte Schalen komprimiert wird, die von starken Stoßfronten begrenzt sind. Schnelle Wolken entstehen im Raum zwischen den Schalen und stoßen mit diesen zusammen. Dies erzeugt Röntgenflashes, die die beobachtete Röntgenstrahlung heißer Sterne erklären können. 2. Wind runway durch radiative Wellen. Der runaway zeigt, warum beobachtete liniengetriebene Winde schnelle, kritische Lösungen anstelle von Brisenlösungen (oder shallow solutions) annehmen. Unter bestimmten Bedingungen stabilisiert der Wind sich auf masseüberladenen Lösungen, mit einem breiten, abbremsenden Bereich und Knicken im Geschwindigkeitsfeld. 3. Magnetische Winde von Akkretionsscheiben um Sterne oder in aktiven Galaxienzentren. Die Linienbeschleunigung wird hier durch die Zentrifugalkraft entlang korotierender poloidaler Magnetfelder und die Lorentzkraft aufgrund von Gradienten im toroidalen Feld unterstützt. Ein Wirbelblatt, das am inneren Scheibenrand beginnt, kann zu stark erhöhten Massenverlustraten führen.

Contents

§ 1 Introduction

Chapter 1. The line force

§ 2 Pure absorption line force. Sobolev approximation

§ 3 Scattering in SSF and EISF approximation

Chapter 2. Unstable winds

§ 4 The de-shadowing instability

§ 5 Evolved wind structure

§ 6 Wind structure and line profiles

§ 7 Including energy transfer

§ 8 X-rays and clouds

Chapter 3. Runaway winds

§ 9 Solution topology

§ 10 Abbott waves

§ 11 Abbott wave runaway

§ 12 Overloaded winds

Chapter 4. Disk winds

§ 13 Analytical model

§ 14 Numerical model

§ 15 Magnetized line driven winds

Outlook and Acknowledgments

Appended Papers

§ 1 Introduction ¹

Why line driven winds? The study of winds from astronomical objects started with Parker's work on the solar wind in 1958, one of the landmark theories in astronomy. At present, four important types of hydrodynamic winds are known in astronomy: the thermal wind from Sun; dust driven winds from red supergiants; line driven winds from blue stars and accretion disks; and magnetocentrifugal winds from accretion disks, either around stars or in quasars.

Thermal winds are accelerated by gas pressure in a hot corona. Parker's breakthrough idea was that the solar wind can be described as a transsonic hydrodynamic flow, instead of a discrete particle flux. Dusty winds are driven by continuum radiation pressure acting on dust grains in a relatively cool environment. Line driven winds are also driven by radiation pressure, yet, in numerous ultraviolet spectral lines. Finally, magnetocentrifugal winds are launched from accretion disks, via centrifugal forces acting along poloidal field lines, or Lorentz forces caused by the toroidal field.

By contrast to geological winds which are essentially horizontal flows caused by pressure gradients and the Coriolis force, the above four astronomical winds directly oppose gravity, and carry away mass and momentum from the central object. For line driven winds, both the mass loss rate, [M\dot], and the momentum rate, [M\dot] v_¥ (with terminal wind speed v_¥), are large. These winds are therefore important in two respects:

Stellar evolution. Hot, massive stars lose a large fraction of their initial mass through winds, and the winds control stellar evolution. Unfortunately, the phases of strongest mass loss, the LBV (luminous blue variable) and Wolf-Rayet phase, are not well understood, and empirical formulae have to be used in evolutionary calculations. The LBV phase may be characterized by the star reaching its Eddington limit (radiation pressure on electrons larger than gravity; Langer et al. 1999).

Star formation. Hot stars often reside in environments rich of gas and dust. The stellar wind enriches the interstellar medium with metals and triggers (bursts of) star formation (Leitherer et al. 1999). Figure 1 shows an aspect of this process which was recently discovered. Shown is a Hubble Space Telescope image, where radiation from an O7 main sequence star may prevent planet formation. The stellar wind is seen blowing off matter from the protoplanetary accretion disk.

Figure 1: NASA press release, April 26, 2001. ``Hubble Watches Planetary Nurseries Being Torched by Radiation from Hot Star. Planet formation is a hazardous process. This snapshot, taken by HST, shows a dust disk around an embryonic star in the Orion Nebula being `blowtorched' by a blistering flood of ultraviolet radiation from the region's brightest star. Within these disks are the seeds of planets. Evidence suggests that dust grains in the disk are already forming larger particles, which range in size from snowflakes to gravel. But these particles may not have time to grow into full-fledged planets because of the relentless `hurricane' of radiation from the nebula's hottest star, Theta 1 Orionis C. In the picture, the disk is oval near the center. Radiation from the hot star is heating up the disk, causing matter to dissipate. A strong `stellar wind' is propelling the material away from the disk.''

We add more reasons why the hydrodynamics of line driven winds is an important and interesting new research area.

Spectroscopy. Fundamental parameters of hot stars like radius and mass can be derived from quantitative spectroscopy of spectral lines forming in their winds (Pauldrach, Hoffmann, & Lennon 2001; Kudritzki & Puls 2000; Hamann & Koesterke 1998a; Hillier & Miller 1998). Most importantly, the star's luminosity follows from the wind-momentum luminosity relation (Kudritzki, Lennon, & Puls 1995). In the near future, hot, massive stars may compete with Cepheids as primary distance indicators. Wind hydrodynamics affects spectral line formation and therefore quantitative spectroscopy in a fundamental way. Observed spectral line features which are indicators for time-dependent flow features are (we mention only the terms, without going into explanations): black troughs, discrete absorption components, enhanced electron scattering wings, bowed variation contours, variable blue edges, and discrete emission components.

New force. Line driving is a new hydrodynamic force, and defines a new class of `radiating fluids'. In astrophysics, radiation hydrodynamics is of similar importance as magnetohydrodynamics. The unique quality about line driving is `Doppler tuning'. Each spectral line can absorb photons in a narrow frequency band only, its width being determined by the thermal speed of ions. After a minute acceleration of the gas by the line force, the Doppler-shifted spectral line can absorb at bluer frequencies, from the `fresh' stellar continuum. The ratio of terminal to thermal speed, which is roughly the maximum Mach number, Ma, of the flow, is therefore of order the width of the UV frequency band divided by the line width: Ma > 100! By contrast, flows driven by continuum absorption of radiation, the solar wind, and MHD flows (if the Mach number with respect to Alfvén waves is considered) have Ma £ 10 in their accelerating regions.

New waves. Actually, Ma £ 10 also holds for line driven winds, with respect to a new type of line driven or Abbott waves (Abbott 1980). These waves are discussed controversially in the literature, and a strict proof of their existence is still missing. They may play a central role in the ubiquitous `discrete absorption components' observed in unsaturated P Cygni line profiles (Cranmer & Owocki 1996); and they may drive line driven winds towards a unique, critical solution.

Range of objects. The importance of line driving is also clear from a list of astronomical objects which share this flow type. Line driven winds occur in O and B stars (Lucy & Solomon 1970), Wolf-Rayet stars (Lucy & Abbott 1993; Gayley, Owocki, & Cranmer 1995; Hamann & Koesterke 2000), central stars of planetary nebula (Koesterke & Hamann 1997), and in A supergiants (Kudritzki et al. 1999). The latter are the optically brightest stars, and are central in extragalactic distance determination. A relatively new idea is that line driven winds are launched from accretion disks. Relevant cases span a huge range, reaching from active galactic nuclei (quasars and Seyfert galaxies; Weymann, Turnshek, & Christiansen 1985) to cataclysmic variables (white dwarfs with a late-type, main-sequence companion; Heap et al. 1978) and young stellar objects (bright B protostars; Drew, Proga, & Stone 1998).

Interaction of line and magnetic driving. Line driven winds from magnetized accretion disks may show an intricate interplay of three driving forces: the radiative, Lorentz, and centrifugal force. The assistance of line driving may help to overcome problems encountered with pure magnetocentrifugal driving. This should be relevant for accretion disks in quasars and young stellar objects, where magnetic fields and radiation fields are strong.

After this general motivation, and before we go over to more detailed discussions in the main text, we give an overview of the results and ideas treated in the following.

Wind instability and X-ray emission. Line driven winds are subject to a new radiation-hydrodynamics instability (Lucy & Solomon 1970). If a fluid parcel experiences a small, positive velocity perturbation, it gets Doppler-shifted out of the absorption shadow of gas lying closer to the radiation source (the photosphere). It sees more light and experiences stronger driving, hence, is further accelerated: an amplification cycle results, termed de-shadowing instability. Since the flow is highly supersonic, perturbations are expected to quickly grow into shocks.

In an important paper, Owocki, Castor, and Rybicki (1988; OCR from now) calculated for the first time the unstable flow structure numerically, along a 1-D, radial ray assuming spherical symmetry. The initially smooth flow is transformed into a sequence of dense shells. Since the shells are accelerated outwards, they are Rayleigh-Taylor unstable and should fragment. In linear approximation, the de-shadowing instability has no lateral component (Rybicki, Owocki, & Castor 1990), and one may expect that the R-T debris maintains a relatively large, lateral scale. However, this is presently mere speculation, since 2-D simulations are lacking due to computational limitations. The fragmented shells are separated radially by broad regions of rarefied, steeply accelerating gas. The thin, fast gas is eventually decelerated in strong reverse shocks on the inner, starward (or disk) facing rim of the shells. The shells propagate outwards at a speed similar to that of smooth, stationary flow. Figure 2 shows the evolved wind structure. The origin of X-ray emission in cloud-shell collisions is also indicated in this sketch.

Figure 2: Expected structure of an unstable, line driven wind. Dense shells fragment via Rayleigh-Taylor instability. Fast clouds collide with the shell fragments, creating X-ray flashes.

The detection of X-rays from early-type stars was one of the major discoveries of the Einstein satellite (Seward et al. 1979; Harnden et al. 1979). For O stars, the ratio of X-ray and total luminosity is L_x/L_bol ~ 10^-7±1. X-ray temperatures are 10⁶ to 10^7.5 Kelvin, two or three orders of magnitude above photospheric temperatures. Early-type stars have no envelope convection zones and are thus quite different from solar-type stars with convection zones, magnetic fields, and hot coronae. It is generally believed that the X-ray emission from hot stars originates in their winds, possibly in shocks which result from de-shadowing instability.

While this idea is consistent with many observational facts, especially the absence of K-shell absorption edges, theoretical modeling encounters severe problems. In a phenomenological model of strong forward shocks in the wind (Lucy 1982b), L_x was found to be a factor of 100 smaller than derived from Einstein data (Cassinelli & Swank 1983). The same problem occured for the hydrodynamic wind models of OCR, which show an X-ray flux deficiency by factors of 10 to 100 (Hillier et al. 1993).

Instead of a continuous X-ray emission from quasi-steady wind shocks, we propose X-rays flashes of short duration. These flashes occur when fast wind clouds collide with a dense, cold shell. The clouds are created via turbulent ablation from a `mildly' dense gas reservoir lying ahead of the overdense, highly compressed shells. At the wind base, roughly one half of the upstreaming gas is quickly compressed into shells. The remaining gas is available for cloud formation at larger heights in the wind. The clouds are accelerated through empty intershell space, until they collide with the next outer shell. Time averaged X-ray spectra synthesised from these models match Rosat observations well. Variations between individual snapshots are large, however, in contrast to the observed constancy of X-ray fluxes.

Since the clouds are of turbulent origin, we speculate that their lateral scale is much smaller than that of shell fragments. In future 2-D modeling, X-ray flux constancy may be achieved via angle averaging over independent, radial wind rays, each with its own cloud-shell collisions taking place.

Wind runaway caused by Abbott waves. After the discussion in Chapter 2 of localized flow features due to de-shadowing instability, we turn in Chapter 3 to the global solution topology of steady, line driven winds. The question we pose is: why do line driven winds from stars and accretion disks adopt a unique, critical solution? In a first attempt to answer this question, we suppress de-shadowing instability by adopting the simplifying Sobolev approximation in the line force. Future work has to unify both aspects, global and local. (A relevant question could then be: can local, unstable flow structure provide the seed perturbations required for global solution transformation?)

In their fundamental work on line driven winds, Castor, Abbott, and Klein (1975; CAK from now) showed that the stationary Euler equation has an infinite number of possible solutions. They come in two classes, `shallow' and `steep'. Steep solutions are supersonic, and cannot connect to the subsonic wind base. Shallow solutions fail to reach large radii, as they cannot perform the required spherical expansion work. CAK concluded that the wind adopts the unique, critical solution which starts shallow and ends steep, switching smoothly between the classes at some critical point. This defines the critical or CAK solution.

Evidence mounts that the argument given by CAK is too restrictive. The question arises whether the CAK solution is still unique if discontinuities are allowed for in derivatives of flow quantities (i.e., kinks). Is the CAK solution a dynamical attractor in the sense of mechanical system theory? And how does the transition between shallow and the critical solution occur?

Radiative or Abbott waves are the key to answer these questions. We show that shallow solutions are sub-abbottic, and Abbott waves propagate inwards towards the photosphere from any point in the wind. Numerical wind simulations published so far assume pure outflow boundary conditions, which apply if all characteristics leave the mesh. The standard argument is that the outer boundary is highly supersonic, at Ma > 100. Instead, we see now that Ma < 1 for Abbott waves along shallow solutions. We show that, by chosing outflow boundary conditions, the numerical scheme is forced to relax to a super-abbottic solution. If the wind is sub-abbottic (yet, supersonic), outflow boundary conditions cause numerical runaway.

Abbott waves define the characteristics of line driven winds, and as such have to be included in the Courant time step, to prevent numerical runaway. This has not be done so far. Accounting for Abbott waves in the Courant time step and applying non-reflecting Riemann conditions at the outer boundary, we find that the numerical scheme converges to shallow solutions from a wide class of initial conditions.

Will line driven winds in nature adopt a shallow solution? We suggest that this is not the case, and identify a new, physical runaway which drives shallow solutions towards the critical one. This runaway depends on new and strange dispersion properties of Abbott waves: negative velocity gradients propagate downstream and outwards, as opposed to upstream propagating, positive velocity gradients. This asymmetry causes systematic evolution of the wind towards larger speeds. The runaway stops when a critical point forms in the flow, and prevents waves to penetrate downwards to the wind base. This is the case for the CAK solution.

A new solution type occurs with generalized critical points. If the wind is perturbed in the sub-abbottic region below the CAK critical point, runaway does not stop on the CAK solution. Instead, the wind becomes overloaded. A new critical point forms, showing up as a kink in the velocity law at which the wind starts to decelerate. This is a direct consequence of overloading, since super-CAK mass loss rates cannot be accelerated through the CAK critical point. Increasing the overloading further, negative flow speeds result, and a steady solution is no longer possible. Shocks and shells form in the wind, and propagate outwards. There is some observational evidence suggesting that the wind of the LBV star P Cygni has indeed a broad, decelerating region.

Winds from accretion disks. The last, forth chapter also deals with a simple Sobolev line force, however, in a more complicated flow geometry. Observations indicate that line driven winds are launched by the ultraviolet radiation field in certain types of accretion disks. For example, P Cygni line profiles from broad absorption line quasars show terminal speeds of 10% the speed of light (Turnshek 1984). The flow is rather massive, with mass fluxes of order one solar mass per year in luminous quasars. Both facts are suggestive of a line driven wind. However, the ionizing radiation from the central source poses a serious problem. Wind gas which is not highly compressed or shielded from this radiation becomes fully ionized, and line driving stalls. Shielding could be provided by a cold disk atmosphere. The wind is launched vertically from the disk, reaches large speeds already within the atmosphere, and escapes on ballistic orbits after being exposed to the central radiation (Shlosman et al. 1985). Alternatively, dense regions of hot, ionized gas (of unspecified origin) may occur between the central source and the wind (Murray et al. 1995). This gas may block ionizing X-rays, but could be transparent to UV radiation. The flow is again launched vertically by local disk radiation, but quickly bends over and makes a shallow angle with the disk when irradiated by the central source. Numerical simulations show that the shielding region may consist of highly-ionized, failed wind (Proga, Stone, & Kallman 2000). Strong observational support that winds from broad absorption line quasars are driven by resonance line scattering comes from ``the ghost of Lya'' (Arav 1996).

While very fascinating, line driven quasar winds are still a matter of debate. On the other hand, observational evidence for line driven disk winds in cataclysmic variables (CVs) is unambiguous (Heap et al. 1978; Krautter et al. 1981). These systems consist of a white dwarf and a late-type, main sequence companion. The latter fills its Roche lobe and feeds an accretion disk onto the white dwarf. Observed P Cygni line profiles clearly indicate a biconical outflow from the disk, not a spherical outflow from the white dwarf itself (Vitello & Shlosman 1993). In Chapter 4, we discuss numerical simulations and a semi-analytical model for CV winds. We encounter a two-dimensional eigenvalue problem, for mass loss rates and wind tilt angles with the disk. The derived mass loss rates are much smaller than was hitherto expected. We discuss why the latter expectations were overestimates. Still, a large discrepancy remains, and it is not yet clear why CV disk winds are so efficient in nature. Possibly, magnetic fields assist line driving. This leads over to the last topic of this writing, magnetized winds.

Magnetocentrifugal winds could occur in young stellar objects and in quasars. In the classical model of Blandford & Payne (1982; see Fig. 3 for a realistic scenario), outflow occurs along poloidal magnetic field lines. If the magnetic pressure in the disk corona is much larger than gas pressure, the field lines co-rotate with the disk. The field lines act as lever arms and, at inclination angles < 60 degrees with the disk, flung the gas outwards. Above the Alfvén point, ram pressure dominates magnetic pressure in the wind. The gas parcels start to conserve angular momentum, and lack behind the disk rotation. The poloidal field gets wound-up into a toroidal field. Hoop stresses of the latter confine the centrifugal outflow to bipolar jets, and jets are indeed often found associated with young stellar objects (Eisloeffel et al. 2000). There is a large body of literature on the theory of magnetocentrifugal disk winds. Some key references include: Pudritz & Norman (1983; 1986), Königl (1989), and Heyvaerts & Norman (1989). An excellent review is Königl & Pudritz (2000).

Figure 3: Realistic scenario for a magnetocentrifugal wind from a quasar disk. Gas clouds are centrifugally driven outwards along sufficiently inclined, co-rotating magnetic field lines. From Emmering, Blandford, & Shlosman (1992).

Recently, Drew, Proga, & Stone (1998) suggested that line driving could assist magnetocentrifugal driving in bright young stellar objects (young B stars). The cold accretion disk is assumed to be illuminated by the hot central star, and either scatters the incident UV radiation or is itself heated to temperatures resulting in strong UV emission. The disk radiation field launches a flow. The terminal speed of line driven winds scales with the escape speed at their base, and is much smaller from outer disk regions than from the central star. This scenario could therefore explain the small observed outflow speeds from certain objects.

At the end of Chapter 4, we present explorative 2-D simulations of magnetized line driven winds. The resulting flow dynamics above accretion disks with small Eddington factor is rather intricate. For realistic magnetic field strengths, mass loss rates may be dramatically increased. We find that dominant magnetic driving is not via the centrifugal force along poloidal field lines, but via the Lorentz force caused by toroidal field gradients. This scenario, complementary to the Blandford & Payne model, was proposed by Contopoulos (1995). Still, we find that the poloidal field is mandatory for increased mass loss. A vortex sheet forms in the poloidal wind velocity and magnetic field, and carries the toroidal field to such heights that it can assist in driving the enhanced mass loss through the critical point (the `bottleneck' of the flow).

Much work remains to be done on magnetized line driven winds, in order to understand the relevant physics and rule out numerical artefacts. This is the reason why no paper or preprint on this subject is appended to the present writing.

CHAPTER 1. THE LINE FORCE

§ 2 Pure absorption line force. Sobolev approximation

For one line. Line driven winds stand and fall with the formulation of the radiative force. The force due to momentum transfer by absorption and re-emission of photons in spectral lines is termed `line force' from now. The basic problem in calculating the line force is the inclusion of scattering in the wind. The absorbed radiation determines together with particle collisions the ionization degree and occupation numbers of the plasma, which in turn determine the emitted radiation field. For a stationary and smooth flow, the feedback between radiative transfer and statistical equilibrium of the gas can be solved using ALI techniques (Cannon 1973; Scharmer 1981; Hamann 1985). However, coupling the radiative transfer and thermal equilibrium to time-dependent hydrodynamics over » 10<sup>5</sup> time steps (instead of `1' for stationary winds) and for » 5000 spatial grid points in a highly non-monotonic velocity law (instead of 50 for smooth flow) is not yet possible. Hydrodynamic simulations have to severely approximate the radiative transfer and thermal equilibrium of the wind. The line force from all driving lines (between 10 and 10⁵) is calculated from the force on a single line via a power law line distribution function with two parameters (Castor, Abbott, & Klein 1975; CAK in the following). For dense O supergiant winds, these parameters are constrained within relatively narrow margins (Gayley 1995; Puls, Springmann, & Lennon 2000). The radiative transfer in the remaining, `generic' spectral line is further simplified. In Chapter 2, we consider pure line absorption and a simplified approach to scattering. In Chapters 3 and 4, Sobolev approximation is used.

The force which acts on a gas absorbing radiation of intensity I is derived in the textbooks by Chandrasekhar (1950) and Mihalas (1978). To avoid angle integrals at the present step, we consider a radial, spherically symmetric flow with coordinate r, which is accelerated by absorption (no scattering) of radiation from a point source located at r=0. The line force per mass, g_l, is

gl(r) = gt(r) ó õ ¥ -¥  dx  f(x-v(r)/vth)e-t(x,r),

(1)

where t is the line optical depth,

t(x,r)= ó õ r 0  dr¢k(r¢) f(x-v(r¢)/vth))

(2)

and

gt(r) = kDnD Fn0(r)/c

(3)

is the radiative acceleration if the line were optically thin, t << 1. A Doppler profile function is assumed, f(x) = p^-1/2 exp(-x²), with normalized frequency variable x=(n-n₀)/Dn_D. Here, n₀ is the line frequency and Dn_D = n₀ v_th/c the Doppler width of the line, with c the speed of light. For simplicity, the thermal speed v_th is assumed to be constant. The radial wind velocity is given by v, and v/v_th in the profile function accounts for Doppler shifts. In an accelerating wind, ions can absorb `blue' photons, n > n₀, in the line transition n₀, since they appear redshifted. Finally, F_n is the radiative flux per frequency interval dn at frequency n, and k is the mass absorption coefficient of the line, in units cm² per gram. For line transitions to the ground state (resonance lines) and to metastable levels, which both dominate the line force, k is to a good approximation constant, as shall be assumed in all the following.

In hydrodynamic wind simulations of the de-shadowing instability, eqs. (1) and (2) have to be calculated as they stand, by explicit quadrature over dr and dx (and over angle; Owocki, Castor, & Rybicki 1988; Owocki 1991; paper [1]). The r and x integrals are time consuming, since high spatial and frequency resolution is required: the thermal band has to be resolved in both r and x. The thermal speed of metal ions is a few km/s in hot star winds, whereas the terminal wind speed is thousand km/s and larger. Hence, thousands of frequency and spatial grid points are required. The angle integral, on the other hand, is cheap for stellar wind simulations with a high degree of symmetry. A one- or two-ray quadrature may be sufficient. This is no longer true for winds above accretion disks, due to the radial temperature stratification of the disk and the complicated flow geometry (Proga, Stone, & Drew 1998; [6]).

Sobolev approximation. In Sobolev approximation, the profile function f is replaced by a Dirac d function: the wind is assumed to accelerate so steeply that a spectral line is Doppler-shifted into resonance with a photon over a narrow spatial range only. This region is called a resonance or Sobolev zone. The quantities r, k and dv/dr are assumed to be constant over the Sobolev zone. By definition, v increases by a few thermal speeds over the zone.

Radiative transfer takes place then on a `microscopic' scale within the zone, hydrodynamics (i.e., changes in r and dv/dr) on a `macroscopic' scale (Rybicki & Hummer 1978). The fact that the wind speed v is neither a truly microscopic nor macroscopic variable causes some difficulties in our understanding of radiative waves and critical points in line driven flows [10]. A characteristic analysis in Chapter 3 shows that in Sobolev approximation, dv/dr and not v is a fundamental hydrodynamic quantity besides r.

To calculate the line force in Sobolev approximation from a point source of radiation, substitute r ® [x\tilde] in the optical depth (2), where [x\tilde] = x-v(r)/v_th. This is allowed for monotonic v(r), and gives

t( ~ x   ) = t0 F( ~ x   ),

(4)

where F([x\tilde])=ò_[x\tilde]^¥ dy f(y), and t₀=rkv_th/(dv/dr) is the total optical depth of the Sobolev resonance zone. The line force at r is caused by lines which absorb at r, hence [x\tilde](r)=0 and t₀=t₀(r). Furthermore, [x\tilde] º ¥ is assumed at the wind base. Substituting x® F in the line force integral (1) gives

gl(r)=gt(r) b(r),         with          b(r)= 1-e-t0 t0

(5)

the so-called `photon escape probability.' For t << 1 and t >> 1, the line force scales as k and r^-1 dv/dr, respectively.

The expression (5) holds also for the Sobolev force when line scattering is included, and the diffuse radiation field is fore-aft symmetric, i.e., the latter does then not contribute to the line force. We will apply a Sobolev force in Chapters 3 and 4. In Chapter 2, we consider the more complicated SSF and EISF forces, which are extensions of the general force (1) to the case of line scattering, using ingredients from Sobolev approximation to calculate the radiative source function.

Force from all lines. To calculate the line force from thousands of spectral lines, the CAK line distribution function is used throughout. By assumption, there is no wind-velocity induced line overlap, and each photon is scattered in one line at most. The CAK line distribution function is given by

N(n, k) = 1 n    1 k0   æ ç è k0 k ö ÷ ø 2-a   ,

(6)

introducing two parameters, k₀ and a, where 0 < a < 1. Integrating (1) over n and k using (6), the total line force becomes (we keep the symbol g_l from the single line force),

gl(r) = G(a)k01-a vth c2   F(r)  ó õ ¥ -¥  dx   f(x-v(r)/vth) ha(x,r) ,

(7)

where G is the Gamma function, F is the frequency integrated flux, and we introduced h º t/k (constant k for each line), with t calculated from (2). Inserting instead for a moment the Sobolev optical depth (4), one finds for the total line force in Sobolev approximation,

gl(r) = G(a)k0 vth (1-a) c2   F(r)  t0(r)-a,

(8)

where we redefined t₀=k₀rv_th/(dv/dr) (in order to avoid adding another subscript 0 to t₀). Precise values for k₀ (or, equivalently, the CAK parameter k) and a must be obtained from a detailed NLTE treatment of the wind. For the present purposes, some universal estimates are sufficient. Puls et al. (2000) derive a = 2/3 from Kramers' opacity law for hydrogen-like ions; this value should apply for dense winds. In thin winds, a £ 1/2 (Pauldrach et al. 1994). The absorption coefficient k₀ corresponds roughly to the strongest line in the flow, k₀=O(10⁸ cm²  g^-1) in dense winds. Alternatively, k₀ can be expressed in terms of an effective oscillator number Q (Gayley 1995),

k0 vth se c = Q G(a)-[1/(1-a)],

(9)

where s_e is the absorption coefficient for Thomson scattering on electrons. In a fully ionized hydrogen plasma, s_e = 0.4 cm²/g. For O supergiant winds, Q » 2000 (Gayley 1995), from which k₀ can be calculated.

§ 3 Scattering in SSF and EISF approximation

First time-dependent hydrodynamic simulations of unstable O star winds by Owocki et al. (1988) assumed a pure absorption line force like that in (1,2) since the de-shadowing instability vanishes in Sobolev approximation. This led to certain unexpected results, most notably, that the time-averaged wind does not adopt the CAK solution, but a steeper solution. (The stationary CAK solution is treated in Chapter 3. Until then, it suffices to know that this solution is unique, has a maximum mass loss rate, and a critical point which is not the sonic point.) The defects of the pure absorption model are still not fully understood (see Owocki & Puls 1999 for new insights), but they vanish when line scattering is included. (The diffuse force scales ~ v_th/v, and vanishes only in Sobolev approximation.) Hence, we turn to a simplified treatment of scattering now.

In the `smooth source function' or SSF method (Owocki 1991), a purely local radiative source function from Sobolev aproximation is assumed. The principal idea goes back to Hamann (1981a), and was adopted in the SEI method of Lamers, Cerruti-Sola, & Perinotto (1987). The remaining `formal solution' of radiative transfer is no more complex than for pure line absorption. Especially, the optical depth for photons which are backscattered to the photosphere is by symmetry related to the optical depth for photospheric photons. No extra integrals are required in the SSF method, and the computational costs are practically the same as for pure absorption. The SSF method accounts for the important line drag effect (Lucy 1984), which stabilizes the flow via the mean diffuse radiation field.

The next step of sophistication beyond pure absorption and SSF is the `escape-integral source function' or EISF method (Owocki & Puls 1996), which accounts for the perturbed diffuse radiation field. Already linear stability theory becomes very complex when perturbations in the diffuse radiation field are included (Owocki & Rybicki 1985). Yet, the basic idea of the EISF method is clear and straightforward: in spherical symmetry, the direct radiative force due to photon absorption in a single line is, including now angle integrals explicitly,

ga(r) ~ ámI*(m) b(m,r) ñ.

(10)

Angle brackets indicate angle averages, m is the cosine of the angle of a photon ray with the radial direction, I_*(m) is the angle dependent, photospheric radiation field, and b(m,r) is the escape probability for direction cosine m (cf. eq. 1),

b(m,r) = ó õ dx  f(x-mv(r)/vth) e-t(x,m,r).

(11)

The diffuse or scattering force, on the other hand, is

gs(r) ~ S(r) ámb(m,r) ñ,

(12)

with isotropic source function S. For pure scattering lines in Sobolev approximation, one derives from (5) and the transfer equation (see Owocki & Rybicki 1985 for details) that

S(r)= áI*(m) b(m,r) ñ áb(m,r)ñ .

(13)

In the SSF method, different expressions are used for the escape probabilities in g_a, g_s and in S: in the g's, b is calculated from the actual, time-dependent flow structure; whereas in S, purely local escape probabilities for a smooth flow are used. By contrast, in the EISF method b from the actual, structured flow is used both in g and S. No quantities are introduced in EISF which did not already occur in SSF. Still, the calculational cost is much larger in EISF since two spatial integrals are required. In the first integral, the source function is calculated over the whole mesh, the second serves for a formal solution. Furthermore, a quadrature over frequency x is required. Owocki & Puls (1996) remark that these 3 integrals can be reduced to 2 again for a single spectral line; but not for a line ensemble.

The EISF approach allows for the first time to study phase reversal between velocity and density fluctuations as a consequence of perturbations of the diffuse radiation field. This subtle yet important effect is discussed further on page . The simulations discussed in detail in the next chapter were calculated using the SSF method.

CHAPTER 2. UNSTABLE WINDS

It is now believed that line driven winds from single stars show structure in all three spatial directions. In the polar (rq) plane, wind-compressed disks may (Bjorkman & Cassinelli 1993; Owocki, Cranmer, & Blondin 1994) or may not (Owocki, Cranmer, & Gayley 1996) form around rapidly rotating stars. In the equatorial (rf) plane, co-rotating interaction regions may cause `discrete absorption components' (Cranmer & Owocki 1996) observed in non-saturated P Cygni line profiles. Most `simply', however, already in 1-D radial flow shocks and dense shells develop due to a new hydrodynamic instability. This instability is the subject of the present chapter.

To give an overview of the field, we start with a bibliography of relevant papers. The instability mechanism is discussed, and results from linear stability analysis are summarized. Because of the complexities of line scattering, linear theory is not complete to the present day. This is not just a mathematical curiosity, but the origin of a fundamental debate. Namely, since the Green's function for the case of pure line scattering is not yet known, the nature of signal propagation in these unstable winds is - mysterious. Signal propagation plays a key role in understanding how the flow adapts to boundary conditions. The main topic of this chapter is the evolved, nonlinear wind structure found from numerical simulations, and its relation to observed X-ray emission from O stars.

§ 4 The de-shadowing instability

History. A new, radiation hydrodynamics instability of line driven flows was first suggested by Lucy & Solomon (1970). The mechanism is similar to one proposed by Milne (1926) for the solar chromosphere. Approximate linear stability analysis was performed by MacGregor, Hartmann, & Raymond (1979), Carlberg (1980), and Abbott (1980). In the former two papers, unstable growth rates were derived, whereas Abbott found a new, marginally stable, radiative-acoustic wave mode. This contradiction was resolved by Owocki & Rybicki (1984), who showed that the opposing results apply on different length scales. Lucy (1984) found that the diffuse radiation field from line scattering causes a drag effect which could prevent the instability. Owocki & Rybicki (1985) derived that complete cancellation occurs only close to the photosphere. A very puzzling result was derived by Owocki & Rybicki (1986), who showed that Abbott waves do occur in absorption flows (they shouldn't), but as a pure mathematical artefact. What are the consequences for Abbott waves if scattering is included? Rybicki, Owocki, & Castor (1990) proved that wind instability occurs in flow or radiative flux direction only. In lateral direction, line drag stabilizes the flow. A linear stability analysis for Wolf-Rayet stars was performed by Owocki & Rybicki (1991) and Gayley & Owocki (1995) in diffusion approximation. They found that unstable growth rates are reduced by the multi-scattering factor. Feldmeier (1998, paper [5]) derived that the instability occurs already in Sobolev approximation, if velocity curvature terms are included. This issue of fore-aft asymmetries when crossing the Sobolev zone is addressed in papers by Lucy (1975; Sobolev vs. Newtonian derivative); Owocki & Zank (1991; radiative viscosity) Gayley & Owocki (1994; radiative heating), and Owocki & Puls (1999; source function depression). First numerical simulations of the evolved wind structure were given in breakthrough papers by Owocki, Castor and Rybicki (1988; OCR) for the case of pure line absorption, and by Owocki (1991, 1992) for line scattering in SSF approximation. Poe, Owocki, & Castor (1990) suggested that the nodal topology of the sonic point in line driven flows causes either solution degeneracy or convergence to a steep, non-CAK solution. The classical paper on steady wind solutions and the critical point topology is Castor, Abbott, & Klein (1975; CAK). The issue of solution topology is again related to the inclusion of line scattering. The perturbed diffuse radiation field was treated in linear analysis by Owocki & Rybicki (1985), and implemented in numerical simulations by Owocki & Puls (1996) via the EISF method. EISF simulations which clarified deeper aspects of the Sobolev approximation were performed by Owocki & Puls (1999). Phenomenological wind shock models of X-ray emission were suggested by Lucy & White (1980), Lucy (1982b), Krolik & Raymond (1985), and MacFarlane & Cassinelli (1989). The energy equation and radiative cooling in post-shock zones was included in hydrodynamic simulations of unstable winds by Feldmeier (1995, [1]), and Feldmeier, Puls, & Pauldrach (1997, [4]) suggested X-ray emission from turbulent cloud collisions. P Cygni lines from structured, unstable winds were calculated by Puls, Owocki, & Fullerton (1993) and Puls et al. (1994).

1 The microscopic mechanism.² The basic mechanism of the instability is as follows. Consider the velocity law v(r) of a line driven wind. Because of its width of a few v_th as caused by thermal motions, v(r) may be called a thermal band. The thickness of the thermal band is characterized by the Sobolev length, L=v_th/dv/dr. This is the natural length scale for line driven winds. Instability occurs for perturbations with wavelength l < L (we shall, however, find below that l > L is also unstable). An arbitrary, positive velocity fluctuation +dv Doppler-shifts a gas parcel out of the absorption shadow of gas lying closer to the star (or the accretion disk). The enhanced radiative flux on the parcel accelerates it to larger speeds and de-shadows it further. Since the single-line force scales as g_l ~ e^-t, the amplification cycle can be written as dv® -dt® dg_l ® dv.

Once the parcel got `kicked out' of the thermal band, further de-shadowing is impossible. Carlberg (1980) concluded that the instability should not cause observable perturbations of the velocity law, but microscopic fluctuations of order v_th only. This is correct for short scale perturbations l < L. We shall see below that long scale perturbations, l >> L, result in Dv >> v_th. It is therefore slightly misleading to say (as is occasionally done) that the instability generally causes short scale structure.

Milne (1926) described a runaway process which he held responsible for the ejection of high-speed atoms from the static solar chromosphere. The process is so similar to the de-shadowing instability that we give a quotation from Milne's paper.

``An atom which, due to some cause or other, begins to move outward from the sun with an appreciable velocity will begin to absorb in the violet wing of the absorption line corresponding to the same atom at rest, owing to the Doppler effect. It will therefore be exposed to more intense radiation, and the atom will be accelerated outwards. It will therefore move still further out into the wing, where it will be exposed to still more intense radiation, and so on, until it eventually climbs out of its absorption line.''

Note that Milne describes a plasma instability for single ions, whereas the de-shadowing instability is a hydrodynamic instability for fluid parcels. The other difference is between a static atmosphere and a wind.

2 Linear instability for absorption lines. Owocki & Rybicki (1984) gave the first, full derivation of instability growth rates. For long scale perturbations l >> L, the imaginary growth rates turn into a real dispersion relation for Abbott waves (or, in the older literature, `radiative-acoustic waves'). The derivation given in equations (1) to (30) of Owocki & Rybicki (1984) is very compact. Hence we refer to this paper, and quote only some results. Starting point of the derivation is the line force from a point radiation source, eqs. (1, 2). Perturbations dv enter in both profile functions f(x). After substitution from spatial to frequency variables; applying Sobolev approximation for the mean flow; and introducing harmonic perturbations dv(r) ~ exp(ikr) obeying WKB approximation, one arrives at the perturbed line force (a subscript 0 refers to the mean flow),

dgl dv = iKw0t0 ó õ ¥ -¥  dx  f(x) e-t0 F(x) ó õ ¥ x  dx¢f(x¢)e-iK(x¢-x).

(14)

Here, w₀=g_t/v_th, t₀=k₀r₀v_th/(dv₀/dr), K=kL (the wavenumber in units of the Sobolev length, L), and F(x) as above. The perturbed line force depends on the profile function f. Using an ingenious integration trick, Owocki & Rybicki (1984) solved the double integral analytically for t₀ >> 1, where it becomes independent of f. The result is

dgl dv = wb ik cb +ik              (t0 >> 1),

(15)

where (introducing an opacity c = rk),

wb=w0  f(xb),             cb = c0  f(xb),

(16)

and the blue edge frequency x_b is defined by

F(xb) º 1/t0.

(17)

For short scale perturbations, k®¥ and dg_l/dv = w_b. This implies instability, since the phase shift between velocity and force perturbations is 0. In the limit k® 0, on the other hand, dg_l/dv=ikw₀/c₀. With 90 degrees phase lag between velocity and force perturbations, this corresponds to marginally stable waves: Abbott waves. For long scale perturbations with finite l, one has (slightly) unstable waves. Again from (15), the growth rate drops as (L/l)² for l >> L.

3 Macroscopic linear instability: 2nd order Sobolev. Instead of going through the above, `exact absorption' analysis, Abbott (1980) applied Sobolev approximation to the mean flow and to velocity perturbations. He finds marginally stable waves which are not affected by the instability. This is odd since we saw above that long scale waves, l >> L, are unstable, if only at reduced growth rates. But for l >> L, Sobolev approximation should apply. Indeed, instability occurs in second order Sobolev approximation, including curvature terms d²v/dr² [5]. The second order Sobolev optical depth t is found to be,

t( ~ x   ,r) = t0(r) é ê ë F( ~ x   ) + vth 2v¢ æ ç è v¢¢ v¢ (r) - r¢ r (r) ö ÷ ø æ è f( ~ x   )-2 ~ x   F( ~ x   ) ö ø ù ú û ,

(18)

with t₀ as before, and primes indicating spatial derivatives. Using (18) to calculate dt in the perturbed line force one finds unstable Abbott waves, with growth rate ~ (L/l)². Despite the different approximations made in the two approaches, the latter growth rates agree to within 20% with those of Owocki & Rybicki (1984).

Equation (18) offers an intuitive understanding of the wind instability at large perturbation wavelengths, which complements the picture of parcels being kicked out of the thermal band described earlier (see Fig. 1 in [5]). When the velocity law experiences a small upward bend, v¢¢ < 0, the optical depth is reduced (we assume that v¢ averaged over the resonance zone is left unchanged). This implies an increase in line force, g_l ~ e^-t. The `elevated' region is accelerated to larger speeds. This means it gets further elevated, and -v¢¢ grows further. An amplification cycle -dv¢¢® -dt®dg_l ® dv ® -dv¢¢ results. The corresponding argument holds for a depression +dv¢¢.

4 Nonlinear steepening: a first look at evolved wind structure. We can already at this stage predict some features of the evolved wind structure, before going into details of hydrodynamic simulations. A depression in v(r) as considered in the last paragraph will eventually cause the mean velocity gradient v₀¢ to become smaller (assumed above to be left unchanged). Since t₀ ~ 1/v₀¢ increases, the blue wing frequency x_b becomes more negative, and the growth rate w_b ~ f(x_b) drops steeply. The depression region will not evolve (`depress') further. Depression regions will always remain close to the stationary flow. On the other hand, elevations of the thermal band will continue to grow until the flow becomes optically thin and no further de-shadowing is possible. In the nonlinear regime, fast gas overtakes slow gas ahead of it, and the velocity law evolves a triangular sawtooth shape. The jumps in the `shark fins' decelerate the gas, hence, are reverse shocks. The robust morphology of the wind velocity law, which shows up in most numerical simulations, is therefore: a sawtooth pattern of steep, accelerating regions and pronounced reverse shocks.

5 How large are the shock jumps? Whereas short scale perturbations l £ L saturate at the microscopic level of a few v_th (Carlberg 1980), large scale perturations l >> L give velocity jumps of order v_¥, the wind terminal speed. Overall, a short scale, noisy structure (v_th) is superimposed on a large scale, coherent tilt of the wind velocity law [5].

6 Scattering: the line drag effect. We consider an effect which proves the fundamental importance of line scattering, and that the pure absorption line case may be a degenerate limit (Owocki & Rybicki 1986). Lucy (1984) found that line scattering, via a so-called line drag effect, may prevent wind instability! The origin of the effect is simple. A wind parcel may again experience a perturbation +dv. The parcel leaves the absorption shadow of gas lying closer to the star, leading to instability. However, the parcel experiences also a stronger, backscattered radiation field from larger heights, because it is Doppler-shifted into resonance with faster gas further out. Hence, the inward push grows. Assuming a plane-parallel atmospheric slab and that the star fills a hemisphere, and furthermore that scattering is fore-aft symmetric, one shows that the direct and diffuse force perturbations in an optically thick line cancel exactly. We leave consideration of the outward pointing diffuse force from gas lying at smaller height than the perturbation site to the reader.

For spherically symmetric flow, the instability growth rate is back to 50% of its pure absorption line value at one stellar radius above the photosphere, and reaches 80% of the absorption value at large radii (Owocki & Rybicki 1985). Still, the line drag effect is of great importance both numerically, preventing fast growth at the inner boundary, and in nature, with regard to the much-speculated `photospheric connection', i.e., whether wind structure grows from photospheric perturbations or is wind intrinsic (Henrichs 1986; Henrichs et al. 1990).

§ 5 Evolved wind structure

After decades of fascinating and frustrating mathematical research on nonlinear growth of fluid instabilities (Landau equation, bifurcation, chaos & catastrophes), computers have shown a pragmatic way into evolved, unstable flow via direct, time-dependent simulations. In this and the next sections, we discuss the structure of fully developed, unstable line driven winds. To model the evolution of the instability, a standard Eulerian grid code is used. For the pure hydrodynamics part, we coded a program following the detailed, technical descriptions given in Hawley, Smarr, & Wilson (1984), Norman & Winkler (1986), Reile & Gehren (1991), Stone & Norman (1992a,b). Some of the techniques are also summarized in [1]. The applied techniques comprise: `consistent' advection (Norman, Wilson, & Barton 1980) using van Leer (1977) or `piecewise parabolic' interpolants (Colella & Woodward 1984) on control volumes of staggered grids; non-reflecting Riemann boundary conditions (Hedstrom 1979; Thompson 1987, 1990) ; tensor artificial viscosity (Schulz 1965); pressure predictor (Norman & Winkler 1986).

As for the radiative line force, I followed Owocki et al. (1988) which includes technical subleties like line-list cutoff; a Schuster-Schwarzschild photospheric layer to prevent unstable growth close at the inner boundary; and the strict one-sidedness of the direct force. The SSF method for treating the diffuse radiation field is described in Owocki (1991, 1992), an unpublished draft from 1990, and in Owocki & Puls (1996). Except for one workshop paper (Owocki 1999), all simulations published so far are in radial direction only, assuming spherically symmetric flow. CPU time requirements for 2-D wind instability simulations are huge. The present work is no exception, and deals with spherical symmetric simulations only.

7 The basic entity of wind structure. Figure 4 (taken from paper [1]) shows the evolved wind structure. The pronounced features are: (1) broad rarefaction regions of accelerating gas; this gas is (2) braked in strong reverse shocks, and fed into (3) narrow, dense shells. (4) The shells propagate into gas which remains close to stationary initial conditions. If forward shocks occur at all at the outer shell edges, they are weak.

As a seed perturbation for wind instability, a coherent photospheric sound wave of period 5,000 sec and pressure amplitude 1 percent is introduced at the inner boundary of the model. The wave period determines the spacing of wind shells. Namely, the perturbation wavelength in the photosphere, l = aT=0.009 R_* (with sound speed a, period T, and stellar radius  R_*), is stretched in the accelerating flow by a factor v_¥/a=90. Indeed, the shell distance far out in the wind is 0.8 R_*, cf the figure.

8 Subshells and overtones. Figure 4 of paper [1] shows that 50 (!) overtones of the photospheric sound wave can be clearly distinguished in the wind. There is no indication of stochasticity at all in the wind, hence the structure is strictly deterministic. This is rather unexpected from glancing at the velocity law between 2 and 5 R_*, in which region the wind appears to be rather chaotic (not meant in the new, technical sense of the word). Up to these heights, collisions occur between multiple shells per photospheric excitation period. As noted in [1], these sub-shells are related to non-linear steepening and harmonic overtones. The dynamical details are still not clear. Figure  shows the mass loss rate in the wind as function of radius and time. After initial transients have died out, the wind settles to a limit cycle. Per photospheric perturbation period, three subshells are created and mutually collide with each other up to around 3 R_*.

Figure 4: Wind mass loss rate as function of radius and time. After initial transients have died out, the wind settles to a limit cycle. Per photospheric excitation period, three subshells are created and mutually collide.

9 Two-stage instability. Even for these periodic models the de-shadowing instability is sufficiently delicate that new details become visible when model parameters are varied. Figure 5 below shows a run which differs from the above one mainly by a shorter perturbation period of 1,000 seconds (besides this, it was performed 7 years later). The difference is marked. In the new model, sub-shell collisions terminate already around 2.5 R_*. Afterwards, the shocks decay quickly, out to 6 R_*. Between 7 and 9 R_*, new shocks occur, which was not the case in the old model.

How can new shocks occur after the instability went into saturation, and left a fully developed flow? What happens here is a second stage of the instability. We noted above that negative velocity perturbations saturate quickly, and the wind remains close to stationary initial conditions in regions ahead of pronounced shells. Looking closely at the region from 4 to 6 R_* in Fig. 5, one sees a growing velocity perturbation: in the evolved wind structure, eventually the pseudo-stationary regions become unstable, too, via secondary perturbations. We expect that the latter are related to overtones. The new perturbations steepen into new shocks. The shocks accelerate, overtake and merge with the `old' shock front from the first growth phase of the instability. This creates a single, strong shock.

Figure 5: Snapshot of an O supergiant wind velocity law and density stratification. A periodic sound wave of 5,000 sec period is applied as photospheric seed perturbation for the de-shadowing instability. Dots show individual mesh points. The numbers in the upper panel refer to the sub-shells per excitation period. Taken from [1].

Figure 6: Snapshot of the wind velocity law, assuming a photospheric sound wave with period 1,000 sec. Two stages of unstable growth are seen. In the second stage, perturbations grow within the quasi-stationary gas ahead of dense shells, and lead to secondary shocks.

10 Lamb ringing. Introducing coherent sound waves as instability seeds allows to study the ultimate fate of linear perturbations from harmonic stability analysis. But how relevant are coherent perturbations physically? Below we shall argue that, to understand the observed X-ray emission from hot star winds, one has to refrain from periodic, coherent perturbations and consider random ones instead. Still, relatively coherent perturbations may be expected from stellar pulsations. They cannot be modeled in 1-D simulations, where only pressure or p modes (sound waves) can be excited. Internal gravity waves or g modes cannot propagate vertically through an atmosphere, since lifting a planar atmospheric layer does not give a buoyancy force.

Interestingly, coherent atmospheric perturbations do actually not require an external piston like stellar pulsations, but could be excited intrinsically in the atmosphere: an atmospheric resonance frequency exists, Lamb's acoustic cutoff.

As the reader may already get weary of deflections from the straight path brought about by `unexpected' hydrodynamic effects shooting in from left and right, we add a historical note, trying to emphasize the importance of the acoustic cutoff. Lamb built on earlier work by Rayleigh (1890), who first derived the dispersion relation for sound waves in an isothermal barometric density stratification. Rayleigh failed to see the physical relevance of the frequency w<sub>a</sub> = a/2H defined by atmospheric parameters (a the sound speed, H the scale height). Only after the detection of a single-frequency, atmospheric response after the big Krakatao volcano eruption, Lamb realized the importance of w<sub>a</sub> as an atmospheric resonance. His argument (Lamb 1908), which is far from trivial, can be found on p. 544 of his `Hydrodynamics' (1932). Especially, he proves that in an isothermal atmosphere, a white-noise sound spectrum evolves into a one-spike spectrum, the spike located at the acoustic cutoff period. The acoustic cutoff is therefore indeed a resonance.

The acoustic cutoff was held responsible for the 5 minute oscillations of Sun (Schmidt & Zirker 1963; Meyer & Schmidt 1967). But the appropriate acoustic cutoff period is 3 to 4 min, not 5 min. Ulrich (1970) developed the alternative and correct theory that the 5 min oscillations correspond to acoustic waves trapped in a resonance cavity reaching from the deep solar interior to the top of the convection zone. This established the field of helioseismology. Deubner (1973) observed a photospheric subsignal of 3 min period, where the wave trains are correlated with the appearance of bright granules. These granules are thought to struck the photosphere from below, and excite Lamb ringing. For a review on solar 5 minute oscillations, see Stein & Leibacher (1974).

Simulations of line driven winds from O stars show a self-excitation of the atmosphere at the acoustic cutoff frequency. Its `perturbing' influence on the outer wind structure was occasionally noted in the literature (e.g., Blondin et al. 1990), but not traced back to its physical cause. We close this excursion by expressing our belief that `Lamb ringing' could be responsible for rather coherent trains of pronounced shells in hot star winds.

11 Reverse shocks. Reverse shocks occur because the instability steepens the velocity law of the wind, and the fast gas is eventually decelerated. As a reminder, we add here the distinction between reverse and forward shocks: forward shocks overtake slow gas and accelerate it. In any reference frame, the shock is faster than both the pre- and postshock gas. Forward shocks occur in explosions. Reverse shocks, on the other hand, decelerate fast gas. The shock is slower than the gas on both sides. A reverse shock occurs if a supersonic stream hits a wall or an obstacle.

12 Rarefaction regions. The rarefaction regions and reverse shocks in numerical wind simulations are quasi-stationary. Seen from a comoving frame, they change only little. Rarefaction regions can therefore, to a good approximation, be identified with (patches of) `steep' wind solutions (Owocki, priv. comm.; Feldmeier et al. 1997c). This type of solution to the stationary Euler wind equation will be discussed in the next chapter.

If a stationary rarefaction region is to feed gas into a shell through a reverse shock, a gas source must exist. Rarefaction regions lie directly above depression regions of the velocity law (negative Abbott half waves). Here, the flow does not evolve, but maintains the initial conditions. The depression region extends inwards, to the outer edge of the next shell.

13 Contact discontinuity. Assuming a Sobolev line force and zero sound speed, a contact discontinuity separates the rarefaction and depression regions. The velocity law is continuous there but has a kink. Using terminology only introduced on page , the wind jumps from the critical, initial conditions to a steep solution. Because of subcritical mass flux in the latter, a density discontinuity occurs at the contact discontinuity. Gas cannot penetrate through a contact discontinuity, hence we are still left without gas source to feed the next outer shell.

Any discontinuity in a derivative of the fluid variables v and r propagates at characteristic speed (Courant & Hilbert 1968). Indeed, the contact discontinuity moves at sound speed. Yet, since a=0, it moves along with the fluid. Allowing for finite sound speed instead, the `gate opens' and the rarefaction region eats slowly into the depression. Constant mass flux is maintained through the rarefaction region. The situation is sketched in Fig. 6, which also shows earlier and later stages in the evolution of the instability. Figure 7 shows results from a numerical simulation.

Figure 7: Evolution of winds structure from the de-shadowing instability (schematic).

Figure 8: How the rarefaction zone eats through the mass reservoir. The plot shows subsequent snapshots of wind density (logarithmic) as function of radius (linear). No scales are given on the axes, since the same structures occur for small or large perturbation periods, at small or large separation between pronounced shells. The dotted lines show the underlying, stationary wind model. The gas reservoir ahead of dense shells is unaffected by instability, and has practically stationary densities.

14 Forward shocks and EISF noise. We come to what may at first seem a rather subtle, technical issue. The central importance of line scattering to line driven wind hydrodynamics will, as we hope, become clearer throughout the subsequent discussions. Owocki & Rybicki (1985) derived from linear stability analysis that the perturbed diffuse radiation field turns anti-correlated density and velocity fluctuations into correlated ones. Anti-correlated fluctuations steepen into reverse shocks, correlated fluctuations steepen into forward shocks. In SSF, only the mean diffuse radiation field is treated, and photospheric perturbations evolve into reverse shocks. If, on the other hand, the perturbed diffuse radiation field is included, the phase lag between velocity and density perturbations may be inverted, and strong forward shocks may occur instead of reverse shocks. This argument was first made by Puls (1994), and stimulated development of the EISF method.

Owocki & Puls (1999) proved from EISF simulations, which include the perturbed scattered radiation field, that forward shocks are not important in the unstable wind. EISF and SSF wind structure are essentially identical, and are both dominated by reverse shocks. The reason is that phase reversal occurs only for short scale fluctuations below the Sobolev length (Owocki & Rybicki 1985). Short scale perturbations saturate at velocity amplitudes of order v_th (see page pageref). Correlated perturbations and forward `shocks' (if at all) appear therefore as short scale, small amplitude noise superimposed on the long scale, large amplitude tilt of the thermal band leading to reverse shocks.

15 Are shells enclosed by forward shocks? Furthermore, forward shocks are not required to enclose the shells on their outer edges, to prevent them from expanding away. The argument is simple (Feldmeier et al. 1997c): shells are geometrically thin and consist of subsonic post-shock gas. Hence, their internal sound crossing time is small compared to the flow time, and hydrostatic equilibrium can be assumed. According to the equivalence principle, the outward directed acceleration of the shell by the line force is indistinguishable from an inward gravity binding a static atmosphere. Hence, the shells are held together without the necessity of outer forward shocks.

16 Shell densities. This same argument can be used to estimate shell overdensities with respect to smooth, stationary gas densities. The inner-shell gas expands via thermal pressure. This thermal pressure and shell acceleration define a (`gravitational') scale height, H, which measures the shell thickness. We approximate the line acceleration by v<sub>¥</sub><sup>2</sup>/n R<sub>*</sub>, where v<sub>¥</sub> is the wind terminal speed and n is `a few'; hence, H=n R_*a²/ v_¥². Assuming that all gas from within a rarefaction region of average thickness l/2 (l being the perturbation wavelength) is fed into the shell; and that l »  R_* for the longest perturbations which can still grow into saturation, the overdensity, o, of shell gas with respect to stationary wind gas becomes,

o » l 2H = 1 2n   æ ç è v¥ a ö ÷ ø 2   .

(19)

For v_¥/a » 100, this is of order ³ 10³, in good agreement with numerical simulations. At large radii, radiative acceleration of the shell ceases, and H grows. The shells should have expanded away, and the wind be homogeneous again by 30 to 100 R_*. This, too, agrees with numerical simulations, see Fig. 6 in [4].

17 Inner shell velocity law. At first surprisingly, the velocity law has a negative gradient inside shells (Owocki 1992; [3]). This is also true for solar wind shells (Simon & Axford 1966). The reason is that the velocity law inside the shell reflects the gas history. Gas lying close to the outer shell edge was shocked earlier than gas lying near the inner edge, close to the reverse shock. The shell velocity law is close to a stationary CAK velocity law, which means that the shell is constantly accelerated on its trajectory. With the gas velocity at the reverse shock increasing in course of time, a negative velocity gradient results inside the shell.

§ 6 Wind structure and line profiles

Motivation. Quantitative spectroscopy of UV and optical lines (especially of Ha, Puls et al. 1996) allows to determine mass loss rates, terminal speeds, and metal abundances of the wind. All this information (plus the stellar radius) flows into the wind-momentum luminosity relation, from which the stellar luminosity is inferred (Kudritzki et al. 1999). The premise is that hot stars will in the near future become a primary distance indicator of the same quality as Cepheids. It is therefore important to understand the influence of wind structure on line formation.

P Cygni lines from 1-D, structured wind models. Puls et al. (1993, 1994) calculated P Cygni line profiles for resonance lines from winds structured by the de-shadowing instability. Radiative transfer in the highly non-monotonic velocity law is solved for by iterating the source function, to account for multiple resonance locations (Rybicki & Hummer 1978). The resulting line profiles agree well with profiles from stationary wind models. This is at first surprising, given the `large amount' of structure in the wind, cf. Figures 4 and 5, and Figure 8 below. The reason for the agreement is that most of the mass in the unstable, structured wind still follows a smooth, CAK-type velocity law. The reasoning is actually more subtle, and we refer to Puls et al. (1993) for an in-depth discussion. Some of the differences between line profiles from stationary and unstable, time-dependent wind models resemble observed variability features like black troughs, narrow absorption components (NACs), and blue edge variability.

Observed line variability, and ideas of its origin. Besides X-ray emission, which is discussed further below, the major observational evidence for pronounced flow structure in winds from hot stars comes from variability in optical and UV spectral line profiles. We give a brief, phenomenological overview of variability features. The books edited by Moffat et al. (1994) and Wolf et al. (1999) are good entry points in the large body of literature.

Of central importance are DACs (discrete absorption components), found in unsaturated P Cygni line profiles of OB stars, and so-called discrete wind emission elements, observed in flat-topped emission line profiles from W-R stars (Moffat 1994; Lépine & Moffat 1999) and O supergiants (Eversberg, Lépine, & Moffat 1998). The DACs come in company of `bowed variation contours' or simply `bananas' (Massa et al. 1995). DACs are pure absorption phenomena, while bananas are modulative. Bananas are explained by dense spiral-arm structures in the wind (Owocki et al. 1995; Fullerton et al. 1997), similar to CIRs (co-rotating interaction regions) in the solar wind. Their period is an even divisor (2 or 4) of the rotation period. The origin of the DACs was first supposed in the de-shadowing instability, but this was later excluded (Owocki 1994). Then CIRs were suspected as their origin, but they can only explain the much faster bananas. Hence, DACs remain enigmatic. A promising idea is that Abbott waves cause velocity plateaus in the wind which lead to enhanced absorption and DACs (Cranmer & Owocki 1996). Interestingly, velocity plateaus and not density enhancements were the first idea to explain DACs (Hamann 1980). Discrete wind emission elements, on the other hand, are thought to be caused by compressible blob turbulence. Whereas turbulence in incompressible fluids leads to an eddy cascade, supersonic or compressible turbulence leads to a shock cascade. The eddy cascade is direct, with big eddies feeding energy into small eddies. The shock cascade, on the other hand, is inverse: strong, fast shocks overtake weak, slow shocks, and merge with them into one strong and fast shock. Hence, energy is transferred from small to large scales. A connection may exist between this compressible turbulence and the wind clouds thought to be responsible for X-ray emission from hot stars. This is discussed in the following.

Wind clumping. Terminal speeds and mass loss rates are the central parameters of line driven winds. They determine the metal entrichment of the interstellar medium, and the star formation rate in starbursts. Furthermore, wind mass loss determines the evolution of hot, massive stars. And finally, knowledge of wind parameters allows to derive stellar luminosities and distances, a matter of prime importance in astronomy. Hence, reliable measurements are needed for v_¥ and [M\dot]. However, the mass loss rates of Wolf-Rayet stars, and probably also of O supergiants, are fundamentally affected by instability-generated flow structure. This was first realized by Hillier (1984), who introduced wind clumpiness to explain observed electron scattering wings of emission line profiles in Wolf-Rayet stars. The line core depends linearly on gas density, whereas the wings develop with density squared. Hence, line wings are pronounced in clumped winds. Quantitative profile fits indicate that the wind gas only fills 10 to 30 percent of the available volume around the star (Hamann & Koesterke 1998b). This leads to a reduction of mass loss rates for W-R stars by factors of 2 to 4, clearly demonstrating the importance of hydrodynamic wind structure for quantitative spectroscopy of hot stars with winds. This reduction in the formerly tremendous mass loss rates of W-R stars also opened the way to new, quantitative modeling of their winds being radiatively driven (Lucy & Abbott 1993; Springmann 1994; Gayley et al. 1995; Gräfener, Hamann, & Koesterke 2000).

§ 7 Including energy transfer

18 Radiative shocks. Wind shocks are radiative shocks, consisting of a narrow viscous layer in which the gas is heated, and a subsequent cooling zone in which the gas cools again by radiative losses. So far we assumed implicitly that radiative cooling is very efficient in the wind, and that shock cooling zones are narrow. Radiative shocks can then be viewed `from far' as isothermal shocks, since both heating and cooling occurs on microscopic, unresolved length scales. In this approximation, a solution of the energy equation in the wind is not required. This is the reasoning which led OCR to undertake isothermal wind calculations.

However, the numerical simulations discussed above make the assumption of isothermality questionable at certain heights above the photosphere, and isothermality seems not justified a posteriori. We found that the wind gas is highly rarefied at the end of a rarefaction region, before it undergoes the reverse shock transition. Efficient radiative cooling cannot be assumed. To find a self-consistent wind structure including the effects of radiative cooling (possibly on long scales), Cooper & Owocki (1992) included radiative cooling in numerical wind simulations for the first time. This led to the strange, unexpected results of unresolved cooling zones, and all shocks were still isothermal. As the numerical mesh was chosen sufficiently fine to resolve cooling zones, the conclusion was that the latter got somehow collapsed.

Advective diffusion. We offer in this and the next paragraph two alternative explanations for cooling zone collapse. The explanations are slightly technical, and the reader primary interested in the main physical argument may want to skip two paragraphs ahead. - Owocki (1993, private communication; see also Cooper 1994) explained shock collapse by advective diffusion. This is a manifestation of Field's (1965) local thermal instability. Consider a propagating temperature jump, i.e., a contact discontinuity. Diffusive errors of the advection scheme spread the jump over a few grid points. The gas at intermediate temperatures cools better than hot gas: in pressure equilibrium, cold gas is denser than hot gas, giving more collisions and stronger cooling. The broadened jump is sharpened again by the different cooling rates, and thereby gets slightly shifted into the hot gas. The jump introduces new diffusive errors, and the cycle repeats. Advective diffusion should occur within the viscous shock layer, which is spread out over » 3 grid points by artificial viscosity. The shock front `eats' then through its own cooling zone. This argument assumes that errors pile up, as they indeed do at a contact discontinuity where always the same gas is located. But this is not the case in a shock transition, where gas passes through. A detailed calculation shows that only a slight modification of the cooling zone results from advective diffusion, but no collapse [1].

19 Oscillatory thermal instability. Besides Field's local thermal instability, a second, global thermal instability occurs in radiative shocks. This instability was found in numerical simulations of accretion columns onto white dwarfs (Langer et al. 1981, 1982) in magnetic (AM Her) cataclysmic variables. The linear stability analysis is due to Chevalier & Imamura (1982), and the instability mechanism is explained in detail in Langer et al. (1982), Gaetz, Edgar, & Chevalier (1988), and Wu, Chanmugam & Shaviv (1992). The instability is of oscillatory type, and causes periodic contraction and expansion of the cooling zone. The contractions are strong, and a fine grid is required to resolve them. On a coarse grid, the contracted zone drops at some point below grid resolution, and an isothermal shock remains [1]. Obviously, there is no thermal instability for an isothermal shock: the shock will not re-expand, and the collapse is permanent. We add as a side remark that oscillating radiative cooling zones show rich dynamics. For example, tiny condensations within the cooling zone can grow into secondary shocks and propagate through the cooling zone (Innes et al. 1987).

Altering the cooling function. We assume in the simulations below that radiative cooling is parameterized in power law form, L = A r² T^d (units erg s^-1 cm^-3). The parameter d is derived from fits to calculated cooling functions (Cox & Tucker 1969). Global thermal instability occurs for d £ 1/2. We find that the above cooling zone collapse is prevented when an artificial, stable exponent d > d_c is used at low temperatures, as is demonstrated in Fig. 6 of [1]. Typically, assuming d = 2 at T < 5×10⁵ K prevents cooling zone collapse in wind simulations. X-ray spectra of O stars, in which we are primarily interested, indicate temperatures between 10⁶ and 10^7.5 K. Hence, X-ray emission should be largely unaffected by the modified cooling function.

20 Shock destruction. How does radiative cooling influence the wind structure? What happens at intermediate and large radii in the wind, when radiative cooling ceases to be efficient in rarefied gas undergoing a reverse shock transition? Using the above method, we find in [1] that isothermality is a good approximation for O supergiant winds out to » 5 R_*. At these radii, reverse shocks are suddenly destroyed, in marked contrast to their gradual decay in isothermal calculations, which continues out to » 20 R_*. A direct comparison is made in Figures 1 and 10 of [1].

21 Shocks becoming adiabatic. The reason for shock destruction is the depletion of intershell gas. Radiative cooling becomes inefficient, and the cooling zones broaden as they become adiabatic. This drives the shock front through the rarefaction region, towards the next inner shell. Since the shock propagates into the preshock gas, the postshock temperature raises, which makes radiative cooling even more inefficient, pushing the shock further into the gas. Actually, this is the very mechanism of the global thermal instability of Langer et al. (1981). Eventually, the shock merges with the next inner shell (this may create forward shocks, cf. page ). Tenuous, hot gas at temperatures > 10⁷ K fills the whole space between shells.

Why does shock destruction occur so suddenly around 5 R_*? At this height, the gas reservoir ahead of a shell, which remains first at stationary conditions, is used up, i.e., was fully fed into the next outer shell. Further sub-shells - we term them `clouds' in the following, for reasons which will become apparent - cannot occur in between shells. These clouds propagated outwards and pushed ahead of them the adiabatic shock front which tried to expand in the opposite, inward direction. Once the clouds cease, the shock front can expand freely through empty intershell space, leaving hot, thin gas.

22 Outer corona. Volume filling factors of hot gas at 10⁶ to 10⁷ K can almost reach unity between 5 R_*, the location of shock destruction, and 20 R_*, where hot gas has significantly cooled by adiabatic expansion [1, 4]. In winds from OB supergiants, X-ray emission from this low density gas is negligible. It is still possible that an X-ray emitting, outer corona occurs in thin winds from B stars near the main sequence. Here, the first strong wind shock can heat large fractions of the gas. UV line profiles show indeed that the cold wind vanishes before it reaches terminal speed (Lucy & White 1980; Hamann 1981b).

§ 8 X-rays and clouds

23 Clouds. The sub-shells or clouds are not quite the accidental, secondary feature as which they were treated so far (cf. secondary shock formation on page pageref.) We will argue in the following that clouds are a primary agent in explaining X-rays from hot, massive stars [4]. We distinguish from now on strictly between shells (or shell fragments) and clouds. Shells (shell fragments) are instability generated, highly overdense as compared to stationary wind densities, and have probably large lateral scale even after Rayleigh-Taylor fragmentation (cf. the sketch in Figure 2). Clouds, on the other hand, are turbulence-induced, have roughly stationary wind densities, and, if turbulence is approximate isotropic, a tiny lateral scale.

Hillier et al. (1993) estimated that OCR wind models should fail by one or two orders of magnitude to produce the observed X-ray emission from hot stars. The reason is the low density of shock heated gas immediately behind the reverse shock. The gas density lies orders of magnitude below stationary wind densities. Note that X-ray emission, as a collisional process, scales with the density squared. We found in [1] that large-amplitude, periodic perturbations in the photosphere lead to chaotic wind structure, whereas strictly deterministic flow results from low-amplitude perturbations (cf. Figures 10 and 12 in [1]). We speculated, therefore, that random or turbulent boundary perturbations, still of small amplitude, could lead to an `active' wind showing enhanced X-ray production. This is indeed the case, and clouds are the means of enhanced X-ray emission. For random boundary perturbations, wind clouds are prevalent, far above the level of a few `overtone' sub-shells in models with coherent base perturbations (Figure 4). Clouds are dense and collide with shells, leading to strong X-ray emission.

24 Langevin boundary conditions. To mimic photospheric turbulence, a velocity perturbation, u, is applied on the inner simulation boundary [4]. The perturbation fulfills the Langevin equation for a continuous Markov chain, du/dt + u/t_c = G(t). This equation is integrated in time-forward manner. The correlation time, t_c, is a free parameter. We choose a value not too far from the acoustic cutoff period. G is the stochastic force, with a white-noise correlation function. The force amplitude is the second free parameter, typically chosen at 30% the speed of sound. This is well below the limit of measured turbulent velocity dispersion in hot star atmospheres (Conti & Ebbets 1977). The power law index of this turbulence model is -2, not too different from the Kolmogorov index -5/3 for eddy turbulence. Some classic papers on the Langevin equation are collected in Wax (1954).

Figure 9: Clouds, marked with filled circles, in a snapshot of the density stratification in an O supergiant wind.

25 Cloud ablation. We expect that the average time interval between the passage of two dense shells is not too different from the acoustic cutoff period. Clouds, on the other hand, are excited at much shorter periods. They form in the dense gas ahead of a shell, i.e., in the decelerating part of a harmonic perturbation where unstable growth quickly saturates. The birth density of clouds is the stationary wind density, and they roughly maintain this density until collision with the next outer shell. In Fig. 8, clouds are marked explicitly in a snapshot taken from a numerical simulation. Figure 9 shows the cloud dynamics, their propagation through empty rarefaction regions between shells, and collision with the latter.

Figure 10: Clouds, marked with +, ×, \vartriangle, \square, and *, are ablated from the gas reservoir ahead of a pronounced shell, propagate through the rarefaction region, and collide with the next outer shell. Compare with Fig. 8 on page pageref.

A model case: z Ori. The O9 supergiant z Ori is a standard object for X-ray observations from O stars. Most recently, it was the first O star for which an X-ray line spectrum was taken with the Chandra satellite. With excellent Rosat observations available at the time (1996), we made z Ori the test case for the turbulent-cloud scenario. Figure 13 in paper [4] shows a snapshot of the wind structure of z Ori. The X-ray spectrum from this snapshot is shown in Fig. 14 of the same paper. The spectrum is calculated using a formal integral approach on the 1-D, time-dependent wind. K-shell opacities are assumed for cold gas, Raymond-Smith emissivities for hot gas. The model flux is only a factor of 2 or 3 below the observed one, but not by factors of 10 to 100 as estimated by Hillier et al. (1993) from OCR models.

26 Cloud-shell collisions. Strangely then, Figure 15 in [4] shows that all the X-ray emission arises from a single shock at 4.5 R_*. The gas density in this shock lies between ±1 dex of the stationary wind density, far higher than in the rarefaction region. What happens at this location? Figure 17 in [4] shows that a collision of a fast cloud with a dense shell takes place at this time, at this radius. The figure shows also that any appreciable X-ray emission from the wind is due to cloud-shell collisions.

27 X-ray flux constancy. From Fig. 17 in [4] it is also evident that no more than one or two cloud collisions take place at any time. This leads to strong variability of the X-ray flux, ±1 dex around mean in spherically symmetric models. By contrast, observed X-ray flux variations for z Ori are a mere few percent. We expect that spherical symmetry is a very poor approximation for turbulent structures. A realistic wind model should instead consist of independent radial cones, each with its own cloud-shell collisions taking place. A few thousand such cones would guarantee the observed flux constancy.

28 X-ray variability. There are only few detections of X-ray variability so far. Berghöfer & Schmitt (1994) report on an episodic raise in the Rosat X-ray count rate of z Ori. This gained wide attention, for example in the science section of the New York Times. The possible origin for the raise in count rates is the breakthrough of a dense shell through an X-ray photosphere. The event was found in one energy band only. Oskinova et al. (2001) propose observations at higher signal-to-noise, to detect different variability levels in soft and hard bands. This would allow to test certain aspects of the shock model.

Berghöfer et al. (1996) find periodic X-ray variability in z Pup, which maybe caused by absorption of X-rays in a co-rotating interaction region (a spiral arm) in the wind. Oskinova (2001, priv. comm.) finds other examples of periodic X-ray variability in ASCA observations of O stars. We recently proposed a long (180 ksec) observation of the O star x Per with the Chandra satellite, to find first signatures of X-ray line variability. The reason for this expectancy is that this star has amongst the shortest recurrence time scale for DACs (Kaper et al. 1999; de Jong et al. 2001).

O stars observed with ROSAT. Dedicated Rosat observations are available for 42 stars in the spectral range from O3 to O9. These data were analyzed assuming a 2-component wind, consisting of cold, X-ray absorbing gas and a random ensemble of X-ray emitting shocks. Density and temperature stratifications of radiative cooling zones behind shocks are included in this approach (Chevalier & Imamura 1982; [3]). These postshock stratifications allow to reduce the number of free fit parameters to the observations from 4 (Hillier et al. 1993) or 3 (Cohen et al. 1997) to 2 (see [3]). Application of a simplified version of the method is described in Kudritzki et al. (1996). We find there that the scatter in the L_x/L_bol relation for O stars is smaller than from Einstein data. We currently try to decide whether this scatter is caused by a dependence of L_x on a second wind parameter, besides L_bol.

29 Adiabatic shells. Paper [3] also contains an analytic treatment of outer, adiabatic shocks, whose cooling length is not small compared to the wind scale. We follow there a classic paper by Simon & Axford (1966) on solar wind shells.

30 Forward shocks, again. We discuss in the appendix of [3] that shock destruction as described above can be mimicked by a simple numerical test. Instead of calculating a fully structured wind model using an unstable line force, we consider multiple adiabatic shells in a flow expanding spherically symmetric, and at constant speed. The merging of reverse shocks with the next inner shell leads there to a train of forward shocks. Hence, it is presently not quite clear whether insufficient radiative cooling causes shock destruction or shock transformation (reverse to forward). This question is of some relevance since, for the O supergiant z Pup, soft X-rays originate from very large radii, far above the location of the final cloud-shell collisions (Hillier et al. 1993; Schulz et al. 2000).

Future work on X-ray emission. What determines the azimuthal scale of shell fragments and clouds? This is presently the most pressing question. The demand on computer time is huge for 2-D hydrodynamic modeling of the wind instability. One has to employ computational tricks, like using a specially designed mesh for a 3-ray radiative transfer (Owocki 1999) or 2nd order Sobolev approximation [5]. Both methods have their caveats: the specially designed mesh stretches too strongly as function of radius, and instability generated structure may `fall between the mesh nodes'. For 2nd order Sobolev approximation, inclusion of nonlocal couplings (shadowing of one shell by another) may be impossible (Wegner 1999).

Another major question is: why are X-rays from O stars universal? Why does the intricate wind hydrodynamics and radiative transfer lead to a simple L_x/L_bol=10^-7? And why does universality break down for thin winds from B stars (Cassinelli et al. 1994) and for Wolf-Rayet stars (Baum et al. 1992; Wessolowski 1996; Ignace, Oskinova, & Foullon 2000)?

31 Symmetric X-ray lines. With the launch of the X-ray satellites Chandra and Xmm, observations of X-ray lines from O stars has become possible. First observations of O stars led to the unexpected result that X-ray emission lines in z Ori und z Pup (the usual suspects) are symmetric and almost not blue-shifted (Schulz et al. 2000; Kahn et al. 2001; Waldron & Cassinelli 2001). The line width is roughly half the terminal wind speed. These results are puzzling if X-rays originate indeed in a dense wind, where one expects different optical depths for photons from the front and back hemisphere of the star, resulting in asymmetry or effective blue-shift of the line profile. However, Owocki & Cohen (2001) have recently shown that even in a homogeneous wind, observed line symmetry can be partially attributed to low instrumental resolution. Even better, recent Chandra observations of z Pup ``show blue-shifted and skewed line profiles, providing the clearest evidence that the X-ray sources are embedded in the stellar wind'' (Cassinelli et al. 2001).

This should be the final `out' for coronal models of X-ray emission from O supergiants (Hearn 1972, 1973; Cassinelli & Hartmann 1977). Coronal model were practically excluded due to missing K shell absorption (Cassinelli & Swank 1983; Cohen et al. 1997) and the missing green coronal iron line (Baade & Lucy 1987). Recently, however, coronal models gained again some attention (Waldron & Cassinelli 2001).

We plan to perform line-synthesis calculations in a clumped wind, where almost all gas is confined to narrow shell fragments. If the lateral scale of the fragments is not too small, photons from the stellar back hemisphere may escape to the observer by passing between neighboring shells (cf. Figure 2). Furthermore, on the back hemisphere one looks onto the X-ray emitting, inner rim of the shell; whereas on the front hemisphere, the shock is hidden behind a dense, absorbing shell. Both effects (lateral escape between neighboring shells; inner shell rim emission only) level off differences between the red and blue line wing.

32 Clouds and EISF structure. We add another technical remark. The turbulent wind clouds described above and in [4] have short length scales of a few mesh points only. They were calculated using the SSF method, which does not account for the perturbed diffuse radiation field. On short scales, the latter may turn anti-correlated fluctuations like clouds into correlated ones. Future EISF simulations are therefore vital to the cloud model. To achieve proper resolution of the thermal band, an adaptive mesh technique will be used.

X-rays in stationary wind models. The focus of the present work is on wind hydrodynamics, with a highly simplified treatment of NLTE (using a power-law line list; and treating one resonance line only) and radiative transfer (using Sobolev approximation, SSF, or EISF). The complementary viewpoint: full NLTE and radiative transfer, and assuming parameterized hydrodynamics is of central importance for quantitative spectroscopy of hot stars. Precise UV fluxes from OB stars are required to calibrate the wind momentum-luminosity relation (Kudritzki et al. 1999); to determine ionizing fluxes in H  ii-regions (Sellmaier et al. 1996); and for population synthesis in starbursts (Leitherer et al. 1999) and in blue galaxies at redshifts z > 3 (Steidel et al. 1996). In synthesising UV fluxes, the full statistical equilibrium problem of the wind is treated, however, assuming stationary flow (Pauldrach et al. 2001, and literature therein). To explain the observed `superionization', inclusion of X-rays is still required (Pauldrach et al. 1994, MacFarlane et al. 1994). The shocks are treated in parameterized form, on a underlying, monotonic velocity law. Furthermore, inclusion of an artifically high microturbulence is necessary to reproduce observed UV line profiles (Hamann 1980). In future calculations on structured wind models, microturbulence and shock strength will no longer be independent, free parameters, but will both appear as consequence of wind hydrodynamics. Lucy (1982a, 1983) and Puls et al. (1994) showed that black troughs in saturated P Cyg profiles, so far explained by microturbulence, indeed result from an instability-generated velocity law.

High mass X-ray binaries. Neutron stars or black holes which orbit an O or B supergiant act as point-like probes in the wind. Bondi-Hoyle accretion of the supersonic wind turns the compact object into a strong X-ray source which ionizes a large volume of the wind gas. Line driving stalls in this overionized region, and a dense wake forms via the Coriolis force. The beautiful and intricate wind dynamics is discussed by Blondin et al. (1990) and Blondin, Stevens, & Kallman (1991). The presence of the dense wake can be infered from asymmetries of X-ray lightcurves at ingress and egress. For the system Vela X-1, the azimuthal extent of the wake is consistent with the observed asymmetry [2]. Kaper, Hammerschlag-Hensberge, & Zuiderwijk (1994) suggest that observed, short-term variability of X-ray lightcurves is due to fluctuations in the wind velocity as caused by the de-shadowing instability. Hence, high-mass X-ray binaries are another test case to diagnose line driven flow structure.

CHAPTER 3: RUNAWAY WINDS

§ 9 Solution topology

We turn now to a simpler level in our description of line driven wind hydrodynamics, and consider the simplest line force possible, in first order Sobolev approximation. The wind is then no longer subject to de-shadowing instability. We consider why line driven winds from stars and accretion disks adopt a unique, critical solution out of an infinite number of possible solutions, the latter falling into two classes: shallow and steep.

This unique solution is defined as follows (Castor, Abbott, & Klein 1975): it starts in the photosphere as the fastest possible, shallow solution, and crosses at some critical point smoothly and differentiably to the slowest, steep solution. From all accelerating wind solutions, this is the one with maximum mass loss rate. Shallow solutions are discarded as global wind solutions because they break down at large radii, becoming imaginary there. Steep solutions, on the other hand, are everywhere supersonic. They cannot connect to the wind base which is assumed to be subsonic, and are also discarded as global wind solutions.

There are indications that this reasoning is too restrictive (see paper [6]): (1) Shallow solutions break down around 300 R_*. Why then do time-dependent hydrodynamic simulations which extend out to 10 R_* always adopt the critical solution, instead of remaining on shallow initial conditions? (2) Spherical expansion work, as a thermodynamic effect, scales with the gas temperature and sound speed. It is not inherent to line driving. (3) For line driven winds from accretion disks, neither shallow nor the critical solution reach infinity, because the disk flux drops off steeper than disk gravity [6]. Unavoidably, the wind has to decelerate at large heights. This is unproblematic, since the wind speed is above the local escape speed, both for the critical and for sufficiently fast, shallow solutions. (4) Once the wind is allowed to decelerate in certain regions, an initially shallow solution could jump to the decelerating branch when the speed is larger than the escape speed, and reach infinity. Wind deceleration is also discussed by Koninx (1992) and Friend & Abbott (1986), the latter authors modeling velocity laws of rotating O star winds.

A simple wind model. We consider now a simple model for line driven winds which allows to discuss the above questions in detail. The model is used in numerical simulations [8, 9, 10] to study wind runaway caused by Abbott waves. In Sobolev approximation, the line force scales as

gl ~ ó õ dw n  In æ è r-1 ¶nvn ö ø a   .

(20)

The integral is over solid angle, n is the unit direction vector, and I is the frequency-integrated intensity of the radiation field. A CAK line distribution function is assumed, and we adopt a = 1/2 for simplicity in the following. Furthermore, the velocity gradient is taken out of the integral, assuming that the gradient of the flow speed in flow direction gives the main contribution. This is equivalent to the CAK point star approximation. Assuming 1-D, planar geometry,

gl(z) ~ F Ö   r-1 ¶v/¶z   ,

(21)

with frequency-integrated, radiative flux F in z direction. We assume that F is constant, but g(z) may be an arbitrary function. Constancy of both F and g leads to solution degeneracy for zero sound speed (Poe et al. 1990). If the disk is extended and isothermal (the latter is a bad approximation), F is indeed constant at sufficiently small z. For g we choose

g(z)= GM q2 z/q 1+(z/q)2 ,

(22)

which roughly resembles effective gravity (gravity minus centrifugal force) close to a thin accretion disk. M is the mass of the central object, q is the footpoint radius in the disk. Note that g(0)=0. Furthermore, the sound speed is set to zero, which turns the Euler equation into an algebraic equation in vv¢. Note that the Sobolev line force is independent of v_th, and therefore of a. The continuity and Euler equation become,

C º . r   + (rv)¢=0, (23) E º . v   +vv¢+g - ~ K   Ö   v¢/r   =0. (24)

Dots indicate temporal, primes spatial differentiation. The flux was absorbed into [K\tilde]. The Euler equation was divided through by GM/q², which means that speed is measured in units of the local Kepler speed, Ö[(GM/q)], height in units of q, and time in units of a Kepler flow time through a distance q. We keep the symbols v, z, and t also for the normalized quantities. For a stationary wind, the continuity equation becomes rv=const. We introduce m=rv/r_c v_c and w¢=vv¢. Subscripts c refer to the critical or CAK solution, to be introduced shortly. The Euler equation turns into,

w¢+ g - K   ____ Öw¢/m   =0,

(25)

For stationary winds, m and w¢ replace r and v as fundamental hydrodynamic variables.

Shallow and steep solutions. At each z, (25) is a quadratic equation in Ö[(w¢)], with solutions

__ Öw¢   = K ±   ______ ÖK2-4gm   2Öm .

(26)

Solutions with `-' are termed shallow, solutions with `+' are termed steep. For sufficiently small m, shallow and steep solutions are globally defined. At the critical point of the critical solution, termed z_c from now, the square root vanishes, and two globally defined shallow and steep solutions merge. By definition, m=1 for the critical solution. For m > 1, shallow and steep solutions become imaginary in a neigborhood of z_c. In this region, gravity g overcomes line driving ~ K. These failing winds are termed overloaded [8]. We may also allow for v¢ < 0, by introducing Ö[(|v¢|)] in the line force. This expresses that the Sobolev force is blind to the sign of v¢. A single decelerating branch results from this generalization. Figure  shows the `solution topology' of the quadratic equation (25) in the zw¢ plane (see also Fig. 3 in CAK; Fig. 3 in Cassinelli 1979; Fig. 4 in Abbott 1980; Fig. 1 in Bjorkman 1995; and Fig. 6 in [6]).

Figure 11: Solution topology of the equation of motion (25) for stationary winds. The right panel shows acceleration solutions, the left panel decelerating ones. For the latter, Ö[(|w¢|/m)] was assumed in the line force. Note the saddle point at x=1, w¢=1/2, which is the CAK critical point for Abbott waves.

From (26), the critical point is a saddle point of E in the zw¢ plane (Bjorkman 1995). Note the important difference that the sonic points for the solar wind and the Lavalle nozzle are saddle points in the zv plane instead. The critical point being a saddle, the shallow solution with m=1 has a discontinuity in v¢¢ at z_c. This discontinuity is avoided by switching to the steep solution. Discontinuities in derivatives of the fundamental hydrodynamic variables lie on characteristics (Courant & Hilbert 1968) . Characteristics are space-time curves along which Riemann invariants or wave amplitudes are constant, or change according to an ordinary, not a partial differential equation. Indeed, the flow becomes super-abbottic at the CAK critical point (Abbott 1980).

The constant K can be expressed in terms of flow quantities at the critical point. Setting the square root in (26) to zero gives, for m=1,

K = 2   __ Ögc   ,

(27)

where g_c=g(z_c). Since the critical point is a saddle of E in the zw¢ plane, ¶E/¶z_c = ¶E/¶w¢_c=0 holds. The latter (`singularity') condition leads to w¢<sub>c</sub>=g<sub>c</sub>. The former (`regularity') condition leads to dg/dz_c=0, hence the critical point coincides with the gravity maximum. This reveals the role of the critical point as bottleneck of the flow. For height-dependent F, g/F² determines the nozzle function instead. If F and g are constant, the critical point degenerates, and each location z becomes critical (for zero sound speed). Inserting k and w¢_c, the general, stationary wind acceleration becomes,

w¢= gc m (1±   _______ Ö1-mg/gc   )2.

(28)

The velocity law v(z) is found by (analytic or numerical) quadrature of w¢=vv¢.

§ 10 Abbott waves

Since Abbott waves are of central importance in the following, we derive them in three different ways: by dispersion analysis [9]; characteristic analysis [10]; and Green's function analysis [10]. The first and last are necessarily linear, but the characteristic analysis holds for arbitrary wave amplitudes. This should not mask the fundamental limitation of our approach: a derivation of Abbott waves including line scattering in any approximation (SSF, EISF) going beyond Sobolev approximation was not given so far, due to mathematical complexities.

Dispersion relation. For all wind solutions, K=2Ö[(g_c)] holds. Inserting this into the time-dependent Euler equation,

. v   + vv¢ = -g + Ö   2rcvc  v¢/r   .

(29)

We consider small harmonic perturbations on a stationary wind solution r₀, v₀ (not necessarily the critical solution), r = r₀+r₁ exp[i(kz-wt)], v=v₀+v₁ exp[i(kz-wt)]. Linearizing the continuity and Euler equation gives,

é ê ë -i w v0¢ +ic+1 ù ú û r1 r0 +(ic-1)    v1 v0 =0, (30) 1 q0 r1 r0 + é ê ë -i w v0¢ +1+ic æ ç è 1- 1 q0 ö ÷ ø ù ú û v1 v0 =0, (31)

where c º k  v₀/v₀¢ and q₀=Ö[(2 m₀ w₀¢)] were introduced. Setting the determinant of the system to 0 gives the dispersion relation w(k). In the WKB approximation, c >> 1, the phase speed and growth rate of the downstream mode (subscript +) become, in the observers frame,

vf = Re(w+)/k = v0,       Im(w+)=-e.

(32)

The small damping term -e is of no further consequence. In the comoving frame, the wave speed is 0 in this zero sound speed limit: the downstream mode consists of sound waves. For the upstream or (-) mode,

vf = Re(w-)/k = v0 æ è 1-1/q0 ö ø ,       Im(w-)=0.

(33)

This new wave type is caused by radiation pressure, and is termed an Abbott wave (after Abbott 1980). For shallow solutions, q₀ < 1, and Abbott waves propagate towards smaller z: shallow solutions are the analog to solar wind breezes. For steep solutions, q₀ > 1, and the waves propagate towards larger z. At the critical point, m=1, w_c¢=g_c, and q₀=1, hence the critical point is a stagnation point for Abbott waves. At smaller (larger) radii, the waves propagate towards smaller (larger) z. For shallow solutions with m₀ << 1 and w¢₀ << w¢_c, Abbott waves propagate inward at arbitrary large speed. Wave propagation along shallow solutions was not treated so far in the literature, as a consequence of Abbott's (as it appears now: erroneous) postulate that ``a line driven wind has no analog to a solar breeze.'' (Abbott 1980)

Abandoning steep solutions. Steep solutions are super-abbottic everywhere, and cannot communicate with the wind base. Hence, their mass loss rate cannot converge to an eigenvalue. Numerical simulations show indeed that steep solutions are stable to almost any perturbations [9]. The perturbations are advected to the outer boundary and leave the mesh. Steep solutions are in the following of not much interest. Still, they may have physical significance, and we encountered them already before: the quasi-stationary rarefaction regions in time-dependent simulations of the de-shadowing instability are steep solution `patches'.

Characteristic analysis. The above dispersion analysis holds for linear waves. We turn now to a characteristic analysis for arbitrary amplitudes. The equations of motion for arbitrary g are written as,

C(r,v) = . r   + vr¢+ rv¢=0, (34) E(r,v) = . v   + vv¢+ g(z) - 2G Ö   v¢/r   =0, (35)

with G º Ö{g_c r_c v_c}. C and E are brought into advection form without further approximation. A first-order system of partial differential equations is called quasi-linear if it is linear in all derivatives of the unknowns (r and v here). Usually, characteristics are defined for quasi-linear systems. Equations (34, 35) are not quasi-linear, due to the presence of Ö[(v¢)]. Courant & Hilbert (1968) show that characteristic directions, a, for general, nonlinear systems are given by

æ ç è -a C[(r)\dot] + Cr¢ -a C[v\dot] + Cv¢ -a E[(r)\dot] + Er¢ -a E[v\dot] + Ev¢ ö ÷ ø =0.

(36)

The symbol a, so far reserved for the sound speed, refers now to any characteristic speed. The distinction should be clear from the context. For quasi-linear systems, the matrix in (36) becomes independent of differentials [(r)\dot], v¢, etc. For the present case, however,

æ ç ç ç è -a+v r 0 -a+v-G/ Ö   rv¢   ö ÷ ÷ ÷ ø =0,

(37)

or, in the observers frame,

a+ = v,             a- = A = v- G rv¢ ,

(38)

with Abbott speed A. For small perturbations,

A = v (1-   ______ Ögc/mw¢   ) = A0,

(39)

with A₀ from the dispersion analysis above. To bring the hydrodynamic equations into characteristic form [10], the system must be made quasi-linear, by applying the following trick (Courant & Hilbert 1968): the Euler equation is differentiated with respect to z, and a new variable f=v¢ is introduced,

. f   + vf¢+ f2 - G Ö fr æ ç è f¢-f r¢ r ö ÷ ø =-g¢.

(40)

This equation is linear in [f\dot] and f¢. The new system consists then of the continuity equation, the Euler equation (40) for f, and the defining relation f=v¢. Re-bracketing (40) and multiplying by r,

r . f   + rAf¢+ rf2 + G Ö fr fr¢ = -rg¢.

(41)

Replacing rf using the continuity equation,

-rg¢ = r . f   + rAf¢- f . r   - æ ç ç ç è v- G Ö fr ö ÷ ÷ ÷ ø fr¢ (42) = r . f   + rAf¢- f . r   - A fr¢ (43) = r2 (¶t + A ¶z ) f r . (44)

The characteristic (or advection) form of the Euler equation is therefore,

(¶t + A ¶z ) v¢ r = - g¢ r .

(45)

We assume that WKB approximation applies, i.e., that the temporal and spatial derivatives on the left hand side are individually much larger than the right hand side, hence, the latter can be neglected. In a frame moving at speed A, the function v¢/r is then constant, and can be interpreted as a wave amplitude propagating with speed A. The Sobolev optical depth is proportional to the inverse of v¢/r, which indicates that Abbott waves are indeed a radiative mode. Introducing f=v¢ in the continuity equation, the latter is already in characteristic form,

(¶t + v ¶z ) r = - fr.

(46)

Since it contains no derivatives of r or f, fr is an inhomogeneous term. The Riemann invariant, r, is no longer constant along the v characteristic: WKB approximation, which would mean r¢/r >> v¢/v, does not apply. Since r scales with gas pressure, this mode can be identified with sound waves.

Green's function for Sobolev line force. Finally, we derive the Green's function for Abbott waves in Sobolev approximation. The Green's function gives the response of the wind to a localized, delta function perturbation in space and time, and is complementary to the harmonic dispersion analysis of Abbott (1980) and Owocki & Rybicki (1984). Since localized perturbations consist of many harmonics, a Green's function describes wave interference. This is clearly seen for water waves, whose Green's function is known from Fresnel diffraction in optics (Lamb 1932, p. 386). For simplicity, we consider a single, optically thick line only, with Sobolev force,

gl(z)=A  ¶v(z) ¶z .

(47)

Density r was absorbed into the constant A. We assume WKB approximation to apply (slowly varying background flow), and consider only velocity perturbations: the Abbott wave amplitude v¢/r from characteristic analysis is not annihilated by this restriction. The linearized Euler equation for small perturbations is

¶ ¶t dv(z,t) = dgl(z,t) = A dv¢(z,t).

(48)

The Green's function problem is posed by specifying as initial conditions,

dv(z,0)=d(z-z0).

(49)

The solution is obtained by Fourier transformation with respect to z. We use the conventions, for an arbitrary function F (bars indicate Fourier transforms),

_ F   (k,t)= ó õ ¥ -¥  dz  e-ikz  F(z,t),       F(z,t)= 1 2p ó õ ¥ -¥  dk eikz  _ F   (k,t).

(50)

Fourier transforming (48) gives

¶ ¶t dv   (k,t) = ikA dv   (k,t).

(51)

The right hand side was obtained by integration by parts, assuming dv(-¥,t) = dv(¥,t)=0. This will be shown a posteriori. The solution of (51) is

dv   (k,t)=beikAt,

(52)

with constant b. Fourier transforming the initial conditions (49),

dv   (k,0)=e-ikz0=b,

(53)

hence,

dv   (k,t)=eik(At-z0).

(54)

Finally, Fourier transforming back to z space,

dv(z,t) = 1 2p ó õ ¥ -¥  dk  eikzeik(At-z0) = d(z-z0+At).

(55)

The delta function perturbation propagates without dispersion towards smaller z, at a speed -A. Also, dv=0 at z=±¥, as assumed above. A is the Abbott speed, as is seen by inserting a harmonic perturbation, dv=[`(dv)]  e^i(wt-kz), in (48), giving for the phase and group speed,

w k = dw dk = -A.

(56)

The result that a delta function peak propagates without dispersion at speed -A could have been foreseen from the phase speed of linear Abbott waves: v_f is independent of l, hence no dispersion occurs and no wave interference between different harmonics. The Green's function, G, is defined by,

F(z,t)= ó õ dz¢ G(z-z¢,t)  F(z¢,0),

(57)

with some arbitrary function F, implying

G(z,t)=d(z+At).

(58)

We have shown that the phase, group, characteristic and Green's function speed for Abbott waves all agree.

The situation seems clear, but isn't.

No Abbott waves for pure absorption. In a landmark paper, Owocki & Rybicki (1986) derived the Green's function for pure absorption line flows subject to the force (1,1). Taking the Sobolev limit in this Green's function, upstream propagating Abbott waves are found.

But this is impossible since, for pure absorption, photons propagate only downstream and an upstream radiative mode cannot exist. Indeed, the bridging law (Owocki & Rybicki 1984) which is used in deriving this Green's function is the sum of a delta function term, expressing de-shadowing instability, and a Heaviside function, expressing downstream shadowing of wind ions caused by velocity perturbations. To derive Abbott's upstream mode from this downstream-only bridging law is a contradiction.

The resolution was given by Owocki & Rybicki (1986), who demonstrated that two limits do not commute in the Green's function analysis: the Sobolev limit l®¥ for long perturbation wavelengths; and the continuum limit å_k®òdk in Fourier transforms of localized functions. Taking the Sobolev limit first, their eqs. (C1, C2) give the same Green's function for inward propagating Abbott waves as we derived above, G(z,t)=d(z+At). If, instead, the continuum limit is taken first (their eqs. B10, C3, C4) and only then the Sobolev limit,

G(z,t)= ¥ å n=0  (At)n n! æ ç è d dz ö ÷ ø n   d(z).

(59)

Formally, this is the Taylor series expansion of a delta function, G(z,t)=d(z+At). The paradox is burried and resolved at this point. The latter equation expresses d(z+At), which is peaked at z=-At, in terms of derivatives of d(z), which vanish exactly at z=-At for t > 0! The series (59) must be taken as it stands. Applying (59) to an arbitrary initial perturbation F(z¢,0) using (57),

F(z,t)= ¥ å n=0  (At)n n!   dn F(z,0) dzn .

(60)

Hence, F(z,t) is constructed from F(z,0), F¢(z,0), F¢¢(z,0), F¢¢¢(z,0), etc. alone, from pure initial conditions at the location z under consideration. No paradox occurs, there are no Abbott waves, only local information is used. The Abbott speed A is a mere constant, without deeper physical significance.

Including scattering. To the present author, this indicates that a bigger frame exists which contains the special case of pure absorption, and gives physical meaning to A. This is achieved when scattering is included. Indeed, numerical Green's function experiments for SSF and EISF forces show an upstream propagating Abbott front (Owocki & Puls 1999). However, due to mathematical difficulties, the Green's function for flows driven by line scattering was not derived analytically so far. And some authors doubt the physical relevance of Abbott waves and the CAK critical point at all, suggesting that they are artefacts of the Sobolev approximation (Lucy 1998). Using an ingenious argument, Lucy (1975) avoided the critical point singularity (0/0) in a numerical quadrature of the Euler equation. He used different numerical representatives for dv/dr in the advection term and in the Sobolev line force. The argument is that dv/dr in the advection term is a true differential, whereas dv/dr in the line force means the velocity difference over the finite Sobolev zone. We adopt here the standpoint: in numerical simulations which apply the Sobolev, SSF, or EISF approximation, Abbott waves must be accounted for in the Courant time step and (outer) boundary conditions. Yet, the physical relevance of Abbott waves remains to be strictly proven by a Green's function analysis including scattering.

Courant time step and outer boundary conditions. In simulations published so far, Abbott waves are not included in the Courant time step, but only sound waves are. We saw above that Abbott waves define the upstream characteristics of line driven winds. Therefore, they determine the numerical time step. If not properly accounted for, Abbott waves cause numerical runaway.

Furthermore, outflow boundary conditions are assumed in the literature on the outermost mesh point. All flow quantities are extrapolated from the interior mesh to the boundary. This is wrong for sub-abbottic shallow solutions, for which Abbott waves enter through the outer boundary. Applying outflow conditions may drive the solution to the critical one, which is super-abbottic at the outer boundary, hence consistent with outflow extrapolation [8].

To maintain shallow solutions, on the other hand, an outer boundary condition must be applied instead of extrapolation. When Abbott waves are included in the Courant time step and non-reflecting boundary conditions are used, we find that shallow solutions are numerically stable [9]. Non-reflecting boundary conditions annihilate any incoming waves (the boundary condition is: `no waves'). Even initial conditions which depart strongly from a shallow wind, e.g., a linear velocity law, converge to a shallow solution.

We now turn the argument around: once numerical stability of shallow solutions is achieved, this allows to explicitly introduce flow perturbations on the interior mesh (away from boundaries) in a controlled way, and to study their evolution, especially stability. We find that a new runaway mechanism exists, which is caused by Abbott waves [8, 10]. In older simulations which applied outflow boundary conditions and did not account for Abbott waves in the Courant time step, this physical runaway was totally outgrown by numerical runaway, and therefore not detected.

Negative velocity gradients. The new, physical runaway occurs in regions where the wind decelerates, v¢ < 0. The Sobolev force is `blind' to the sign of v¢. All what counts is the relative Doppler shift between neighboring gas parcels, which determines the width of the Sobolev resonance zone. The Sobolev line force is then generalized to g_l = 2GÖ{|v¢| /r}, where a = 1/2 is still assumed. However, flow deceleration implies a non-monotonic velocity law, and multiple resonances occur. The incident light is no longer given by the photospheric flux alone, and the constants K in (25) or 2G above become velocity dependent. Rybicki & Hummer (1978) generalized Sobolev theory to account for non-monotonic velocity laws. While introducing interesting, non-local aspects to Abbott waves (action at a distance!), the method of Rybicki & Hummer is not analytically feasible, and we do not consider it further here. Instead, besides the purely local force g_l ~ Ö[(|v¢|)] we also consider g_l ~ Ö{max(v¢,0)}: all incident light is assumed here to be completly absorbed at the first resonance. The true line force should lie between these extremes.

For g_l ~ Ö[(|v¢|)], the Euler equation has characteristic form

(¶t + A ¶z ) v¢ r = - g¢ r ,       where       A = v -± G Ö ±rv¢

(61)

is the Abbott speed in the observers frame. The upper resp. lower sign applies for v¢ > 0 resp. v¢ < 0. Therefore, if the flow switches from acceleration to deceleration, Abbott waves turn from propagating upstream (- characteristic) to propagating downstream (+ characteristic).

For g_l ~ Ö{max(v¢,0)}, the characteristic Euler equation for v¢ < 0 is,

(¶t + v ¶z ) v¢ r = - g¢ r              (v¢ < 0).

(62)

Since the line force vanishes, this is just the ordinary Euler equation at zero sound speed. The upstream Abbott wave turns here into an upstream sound wave. Note that a reversal into a downstream (Abbott) wave does not occur. We conclude that for both types of line force, regions v¢ < 0 cannot communicate with the base at z=0. Abbott (and sound) waves which originate in these regions propagate only outward.

§ 11 Abbott wave runaway

Figure 10, which is taken from [8], shows runaway of a shallow wind velocity law when a sawtooth-like perturbation of amplitude dv and period T is applied at a fixed position z=2 (at one mesh point) at all times. The amplitude dv is sufficiently large that negative v¢ results.

Figure 12: Abbott wave runaway for a shallow wind velocity law (upper panel), and stable Abbott wave propagation along the critical CAK solution (lower panel). The critical point lies at x=1.

The sequence of events which lead to runaway should be as follows [8, 10]: during the half-period of negative velocity perturbations, the slope v¢ < 0 is to the left (at smaller z) of the slope v¢ > 0. Regions with v¢ < 0 propagate outwards, regions with v¢ > 0 inwards. Hence, the two slopes approach and annihilate each other, leaving an essentially unperturbed velocity law. In some more detail: the sawtooth perturbation dv is established via many small, roughly constant dv << dv applied over many time steps dt << T. During the negative half wave -dv, at each time step dt, -dv is applied, and largely annihilates itself during the subsequent `advection' step of the numerical scheme. The total velocity perturbation after T/2 is not the perturbation amplitude dv, but roughly 0! Opposed to this, over the half-period T/2 when dv > 0, the slope v¢ < 0 is to the right of the slope v¢ > 0. The two slopes move apart. The gas in between gets rarefied (is `stretched') at roughly constant speed. (This corresponds to a centered rarefaction wave.) At the next time step, dv is created atop of this region. Over the half-period T/2, the full perturbation amplitude dv builds up.

This is shown in the figure: the negative half wave of the sawtooth annihilates itself, leaving an unperturbed velocity law; the positive half wave is excited at the intended amplitude, and spreads upstream and downstream throughout the wind. Asymmetric evolution, or runaway, of the wind velocity towards larger speeds results over a full period. Note that the runaway is not caused by wave growth, as in a fluid instability, but by missing negative velocity perturbations to compensate for positive perturbations.

The same runaway is also found by introducing a train of sawtooth-like perturbations in the initial conditions, and let them evolve freely, without introducing further perturbations at later times. This is shown in Figure 5 of paper [10]. Neighboring slopes approach and separate in the way described above, which leads to asymmetric evolution of the wind velocity law towards larger speeds.

§ 12 Overloaded winds

Generalized critical points. The runaway stops when Abbott waves can no longer propagate inwards, after a critical point formed in the flow. If the perturbation source is located above the CAK critical point, the wind undergoes runaway until it reaches the CAK solution. The perturbation site comes to lie on the super-abbottic portion of the velocity law. Phases v¢ > 0 and v¢ < 0 of the perturbation cycle combine then to a smooth, outward-propagating Abbott wave, as shown in the figure.

The perturbation source is not required to lie at fixed z. Even when moving outwards with the local wind speed or faster, it can excite Abbott waves which propagate inwards in the observers frame, and therefore alter the whole inner flow.

For coherent perturbations below the CAK critical point, runaway does not terminate at the critical solution [8]. The perturbation site is sub-abbottic on the critical solution, and Abbott waves can further penetrate to the wind base. The wind evolves towards a velocity law which is steeper than the critical one. The wind becomes overloaded, m > 1. The runaway continues until a generalized (non-CAK) critical point, z_a, forms as a barrier for Abbott waves.

What are these generalized critical points, z_a? To see this, we set the (stationary) Abbott speed (39) to zero, A_a = v_a(1-Ö[(g_c/ m w¢_a)]) = 0, giving mw¢_a=g_c. This shows that the square root in the solution for w¢ in (28) vanishes, and a steep and shallow solution merge. Besides at the CAK critical point, this happens only along overloaded winds, when a shallow solution bends back towards smaller z on the steep branch. Below this location, the Euler equation has two real solutions, above it has two imaginary solutions. Physically, this means that in the proximity of the CAK critical point, which is the bottleneck of the flow, the line force can no longer balance gravity, and the wind starts to decelerate. The two real solution branches reappear at some height above the critical point.

Stationary overloading: kinks. Therefore, at z_a the wind jumps to the decelerating branch, v¢ < 0. In Sobolev approximation, the jump is sudden and causes a kink in the velocity law. Since v¢ < 0 on the decelerating branch, Abbott waves propagate outwards beyond z_a, hence z_a is a wave barrier [8].

We have arrived at the somewhat strange conclusion of a stationary solution with kinks. As mentioned before (page pageref and pageref), kinks and other high-order (or weak) discontinuities of hydrodynamic variables lie on characteristics, i.e., move at characteristic speed. This is true for the present kink, which moves at characteristic Abbott speed A_a=0.

Time-dependent overloading: shocks and shells again. Since z_c lies at the gravity maximum, already small super-CAK mass loss rates cause broad deceleration regions. In practice, for mass loss rates only a few percent larger than the CAK value, the deceleration regime is so broad that negative speeds result. Stationary solutions are then no longer possible, as upwards streaming wind gas collides with falling gas. For a periodic perturbation source in the wind, numerical simulations show the occurence of a train of shocks and shells in the wind, which still propagates outwards. This is shown in Fig. 8 of [10].

Does overloading occur in nature? The discussion so far is idealized since it assumes coherent perturbation sources in the wind. Perturbations with a finite life time should still lead to `piecewise' runaway, each perturbation lifting the wind to a slightly faster, shallow solution. The wind waits on the new, stable shallow solution until the next perturbation lifts it further.

The outer wind seems to be the natural seat for runaway perturbations, since v¢ < 0 is easily excited on the outer, flat velocity law. We speculate that in nature, runaway perturbations are most prevalent above the CAK critical point and drive the wind towards the critical, not an overloaded solution [10]. In support of this we add that not even the strong de-shadowing instability leads to significant wind structure below the critical point, which could act as runaway perturbation towards an overloaded solution.

Still, there remains a slight chance that overloading occurs in real winds. Most notably, Lamers (1998, private communication) reports on a broad region in the wind of the luminous blue variable P Cygni, which is indicative of flow deceleration; and therefore, possibly, of overloading.

Instability vs. runaway. We close this chapter with a brief summary of the `driving agents' of wind dynamics in the last two chapters: de-shadowing instability and Abbott wave runaway. The de-shadowing instability is a true fluid instability (an amplification cycle) which acts on infinitesimal velocity perturbations. Even perturbations of short duration will lead to pronounced, non-linear flow structure which is advected outwards with the wind. In Sobolev approximation, the instability occurs from second order on (including velocity curvature terms). Abbott wave runaway, on the other hand, requires finite amplitude perturbations, since v¢ < 0 is required. There is no amplification cycle, but only a kind of perturbation `filtering' as consequence of the asymmetry of the line force with respect to the sign of v¢. Negative velocity perturbations annihilate themselves, and give way to systematic flow acceleration towards the critical solution. The runaway requires persistent perturbations. Yet, it occurs already in first order Sobolev approximation.

We should also mention that the present, hydrodynamic runaway is not related to the plasma runaway caused by frictional decoupling of line driven metal ions and dragged-along protons (Springmann & Pauldrach 1992). Still, a connection may exist between these two runaways: Krticka & Kubát (2000) report that frictional decoupling in thin winds may actually be prevented by the wind jumping to a slow solution with shallow velocity gradient. Issues of Abbott wave propagation and multiple solution branches become again interesting in these two-component fluids.

CHAPTER 4: DISK WINDS

Line driven winds from accretion disks are a relatively new research area, with quantitative modeling starting in the mid 80ies (Shlosman et al. 1985; Weymann et al. 1985). Disk winds are fascinating because of the richness of their environments, including quasars, Seyfert galaxies, cataclysmic variables, and protostars. In these environments, high-energy particle processes, magnetic fields, coronae, jets, accretion-dominated advection, compact objects, hot boundary layers, and dust formation play important roles.

Theoretical modeling of line driven disk winds is still in its infancy. In the present chapter, we treat first a simple, analytic model [6, 7] for line driven winds from cataclysmic variables. Observational evidence for line driving is almost unequivocal in these objects. In later sections, we discuss winds from magnetized accretion disks, mainly aimed at young protostellar objects and cataclysmic variables. Here, interplay of three driving forces (centrifugal, Lorentz, and line force) leads to an intricate wind dynamics.

§ 13 Analytical model

1-d description, flux and gravity. We set up a simplified, 1-D model for stationary, line driven winds from thin accretion disks (no self-gravity) in cataclysmic variables (CVs). These systems consist of a white dwarf and a late-type main-sequence companion, the latter filling its Roche lobe. We consider first non-magnetic systems (DQ Her class), where the accretion disk should (almost) reach the compact object. The clearest indication of line driving comes from the fact that dwarf novae develop P Cygni line profiles during outburst (Krautter et al. 1981; Klare et al. 1982; Córdova & Mason 1982). Hence, the disk radiation field and wind are causally connected.

A number of approximations is made. First, we assume that the disk is geometrically thin, a plane of zero width indeed. We assume that the basic CAK formalism (Sobolev line force; statistical line distribution function) applies for these objects. Our central assumption is that each helical wind trajectory lies in a straight (yet, not vertical) cone. This is sketched in Figure .

Figure 13: Flow geometry in a semi-analytic model of line driven winds from accretion disks in cataclysmic variables. The helical streamlines are assumed to lie in straight cones. Only the velocity component v_l along the cone is solved for. The tilt angle, l, of the cone with the disk is a function of radius, and is derived as an eigenvalue.

A realistic approximation for the temperature run with radius q in the disk is T ~ q^-1/2 (Horne & Stiening 1985; Rutten et al. 1993), slightly shallower than the famos T ~ q^-3/4 of Shakura & Sunyayev (1973). We assume that the disk is optically thick, and that the radiation field at each radius is that of a black body. The radiative flux, F, above a disk with T ~ q^-1/2 can be derived in closed form in cylindrical coordinates [6],

F(r,z)=(Fr,Fz)=pI(r,0) r2  z r2+z2   ×                          (63) æ ç è 3r2-z2-q2 2rÖB - r z2+r2 lnC,   3z2-r2+q2 2zÖB - z z2+r2 lnC ö ÷ ø ê ê ê rd q=rwd  , (64)

where r_wd is the white dwarf radius, r_d the outer disk radius, and

B=(r2+z2-q2)2+4z2 q2, (65) C= (z2-r2)q2+(z2+r2)2+(z2+r2)ÖB q2 . (66)

The flux isocontours are shown in Figures 3 and 4 of [6]. We adopt again the (`radial streaming') approximation that only velocity gradients along the flow direction are important. A qualitatively new feature of disk winds is their run of effective gravity with height, i.e., gravity after subtraction of the centrifugal force. In the disk plane, effective gravity is zero due to Kepler balance. Effective gravity increases linearly with height, reaches a maximum, and drops off with r^-2 far from the disk. The exact formula is given in equation (9) of [6].

At first, it was thought that the existence of a gravity maximum makes a standard CAK critical solution for disk winds impossible, and that ionization gradients have to be included (Vitello & Shlosman 1988; [7]). The argument was that, because of zero effective gravity in the disk, a too large mass is launched. Higher up in the wind, the line force cannot carry the gas over the gravity hill, and the gas falls back to the disk. From the discussion in the last chapter we see why this argument does not apply. While the critical point (at finite height) is indeed the bottleneck of the flow, Abbott waves can adjust the disk (which is the wind base) to the correct, maximum possible mass loss rate.

Wind tilt angle as an eigenvalue. In the terminology of the Lavalle nozzle, the effective `area function' along a wind ray is f²/g, if a = 1/2. Here, f and g are the radiative flux and effective gravity, respectively, normalized to their footpoint values in the disk. The tilt angle, l, of the straight wind cone (with length coordinate l) with the disk midplane is found as a second eigenvalue, besides the mass loss rate. For each disk ring, the 2-D eigenvalue problem is posed as,

d . M   = C max l  min l  f2 g .

(67)

Here, d[M\dot] is the mass loss rate from the ring, and C is a constant. We search for the maximum mass loss rate with respect to l, which can be driven through the bottleneck. Hence, the minimum must be taken with respect to l. Equation (67) defines a new saddle point in the ll plane, besides the standard CAK saddle in the lw¢ plane.

Isocontours of f²/g above an accretion disk are shown in Figures 6 and 9 of [6]. Opposed to the stellar wind case, the solution topology is rather intricate here, with multiple saddles and extrema. This is a familiar situation for complicated area functions, and is also encountered when one allows for energy and momentum input at finite height in the solar wind (Holzer 1977).

We find that the accretion disk wind undergoes a rather sudden transition, at about 4 white dwarf radii in the disk, from steep tilt angles of only 10 degrees with the disk normal, to tilt angles of 30 to 40 degrees with the normal. The overall wind geometry is shown in Fig. 10 of [6]. The wind is strongly bi-conical everywhere, which is one of the main deductions from kinematic model fits to observed P Cygni line profiles (Shlosman & Vitello 1993).

Mass loss rates. The mass loss rate from the disk is obtained by integrating over all rings. For a = 1/2, this can be written as,

. M   = c1(g) c2(f)   QG   L c2 .

(68)

Here, L is the disk luminosity, and L/c² is the mass loss due to a single, optically thick line (the photon mass flux!), Q » 2000 is an effective oscillator strength (Gayley 1995), and G is the disk Eddington factor. The correction factors c₁ and c₂ account for the run of gravity and radiative flux. c₁=c₂=1 gives the mass loss rate from a point star. For a thin accretion disk, we find in [6] that c₁ = 3Ö3/2. This is larger than unity since centrifugal forces assist in launching the wind. c₂ has to be calculated numerically, and is also ³ 1.

We encounter the problem that mass loss rates from the above formula are one to two orders of magnitude smaller than values derived from P Cygni line fits ([7]; Vitello & Shlosman 1993; Knigge, Woods, & Drew 1995). Here, the values for [M\dot] obtained from fitting were already revised downwards, for the following reason. In the original line fit procedure, individual disk rings were assumed to have a blackbody spectrum (as we did above). The strong Lyman continuum, however, leads then to large ionization rates, which must be balanced by large recombination rates, or large [M\dot]. If realistic spectra (Long et al. 1991, 1994) with suppressed (or missing) Lyman continua are used instead, [M\dot] can be reduced by roughly one order of magnitude, while maintaining the ionization parameter.

Still, there remains the discrepancy noted above, between calculated and observationally deduced mass loss rates. One possible resolution of this problem is that distances and luminosities of CVs are systematically underestimated (Drew, Proga, & Oudmaijer 1999). Furthermore, the all-decisive parameter Q was so far only calculated for dense O star winds. The different spectral shapes of accretion disks, and ionization effects like the bistability jump (Pauldrach & Puls 1990) may cause larger Q values in accretion disks than in O stars. The strongest candidate, however, for resolving the discrepancy is the usual suspect: magnetic fields.

§ 14 Numerical model

The analytic model discussed so far is one dimensional, assuming straight wind cones. First time-dependent, numerical 2-D simulations of line driven winds from accretion disks were performed by Proga, Stone, & Drew (1998). Their wind rays are surprisingly straight in the polar plane (rz plane), and mass loss rates agree well with the above values. A new result from numerical simulations is the occurence of wind streamers. The wind above the disk is structured into alternating dense and rarefied regions. The dense streams propagate radially outwards, with a speed relative to the disk which is proportional to the local Kepler speed. During this motion, the angle between dense streamer and disk becomes ever shallower. This is shown in Fig. . A movie can be found at internet URL http://www.astro.physik.uni-potsdam.de/ ~ afeld). The occurence of alternating dense and rarefied regions is most probably related to mass overloading.

Figure 14: Dense streamers in a line driven wind above a cataclysmic variable disk (log density, dense regions in black). The white cone around the disk axis is a nearly gas free region with very fast parcels.

§ 15 Magnetized line driven winds

We turn to accretion disks which are threaded by external magnetic fields. This could be the case in CVs (so-called polars or AM Her systems), young stellar objects (YSOs), and quasars. Our present simulations aim at the first two classes of objects only. The results are still preliminary, and much work remains to be done. The models were calculated using the publically available Zeus 2-D MHD code (Stone & Norman 1992a,b), augmented by an own, Sobolev line force routine (Feldmeier, Drew, & Stone, in preparation).

Boundary conditions. The outcome of numerical hydrodynamic simulations depends strongly on the applied boundary conditions. We must, therefore, discuss boundary conditions in some detail. The reader interested mainly in results on disk winds may skip over the next few paragraphs, to the section entitled `magnetized wind scenarios'.

For the present discussion, we chose again cylindrical coordinates (in the code, however, spherical coordinates are used, to avoid staircasing of the stellar surface). On the polar axis, anti-reflecting boundary conditions are used. Empirically, one finds that above the `dark star' (both in CVs and YSOs, stellar UV fluxes are negligible), an essentially gas-free cone forms near the polar axis. The highly rarefied gas in the cone is accelerated to large speeds (a `steep' CAK solution is adopted), and the Alfvén speed, v_A=B/Ö{4pr}, is very large in this cone. Both effects limit the Courant time step. To avoid effective stopping of the code, we let an artificial, vertical ram pressure jet originate on the star (Krasnopolsky, Li, & Blandford 1999). Along the stellar surface and the outer mesh boundary circle, inflow resp. outflow boundary conditions are chosen. Note that outflow boundary conditions are wrong close to the disk plane, below the Alfvén surface, and exert artificial forces on the wind (Ustyogova et al. 1999). Future work has to account for this.

Remains the disk itself. The inner disk structure is not resolved in our simulations (actually, zero sound speed is assumed), hence the disk is a simple, planar boundary. Given are a total of 7 hydrodynamic fields: r, v_z, v_r, v_f, B_z, B_r, B_f, minus one constraint, div B=0. Hence, 6 wave modes result: 2 poloidal and 2 toroidal Alfvén waves, and 2 magnetosonic waves (up- and downstream mode in each case). Some of these waves are modified by the radiative line force. For simplicity, however, we keep their above names, with the one exception of the fast magnetosonic mode, which is termed a magnetoabbottic mode from now on.

By assumption, the disk shall be supersonic, subabbottic, and subalfvénic. Then 1 poloidal, 1 toroidal Alfvén mode and the slow magnetosonic mode enter the mesh from the disk, and 1 poloidal, 1 toroidal Alfvén mode and the fast magnetoabbottic mode enter the disk from the mesh. Hence, 3 extrapolations can be applied on the disk boundary, and 3 boundary conditions must be specified. For the poloidal components of the gas speed and magnetic field, we follow largely Krasnopolsky et al. (1999). Here, special care is taken to avoid kinks when magnetic field lines enter the disk. These kinks would cause artificial currents and forces. We do not follow Krasnopolsky et al. (1999) in case that gas falls back to the disk. In contrast to their reasoning, this should only affect the slow magnetosonic mode.

Remain the toroidal fields v_f and B_f. With regard to the toroidal Alfvén mode, one extrapolation and one boundary condition have to be applied. We fix v_f to be the Kepler speed. This avoids dramatic events of sub-Keplerian gas falling towards the central object, which are encountered when v_f is left free (Uchida & Shibata 1985; Stone & Norman 1994). B_f is extrapolated. This causes a problem, since Lorentz force terms ~ ¶(r B_f)/¶r occur in the Euler equation for v_f. For B_f ¹ 0, v_f does generally not agree with the Kepler speed. We proceed by postulating B_f ~ 1/r in the disk (Ouyed & Pudritz 1997a,b). Whereas the latter authors leave the toroidal disk field constant at all times, we calculate, at each time step, [`B]<sub>f</sub> as average of B<sub>f</sub> on the first mesh row above the disk. Within the disk (boundary), we set B<sub>f</sub>=[`B]_f R/r, where R is a free parameter.

This introduces an infinite signal speed, since averaging B_f over the whole disk surface, at each time step, implies instantaneous couplings. Still, the above procedure works well empirically, and B_f evolves smoothly in both r and z direction. By contrast, simulations assuming B_f=0 in the disk crash.

Magnetized wind scenarios. As should be clear from this lengthy discussion, B_f is of central importance in our simulations. We do not assume a magnetically dominated, co-rotating, force-free corona a priori. If a strong, line driven flow is launched from the disk, the angular-momentum conserving gas parcels may instead drag the magnetic field along, building up a strong toroidal field.

There exist two complementary, magnetocentrifugal wind scenarios. In the model of Blandford & Payne (1982), rigid poloidal field lines in a co-rotating corona act as lever arms transferring angular momentum to the gas. There act no Lorentz forces in the launching region of the wind. In the model of Contopoulos (1995), on the other hand, vertical gradients of a toroidal field, B_f, cause non-vanishing curl (pointing towards the central object), and a Lorentz force rot B_f ×B_f which points vertically upwards, see Figure . The generation of dominant, large scale toroidal fields in disk-plus-corona simulations is treated in Miller & Stone (2000) and Elstner & Rüdiger (2000).

Figure 15: Two basic scenarios for magnetized disk winds, launched either by poloidal (left) or toroidal (right) magnetic fields. To apply the right-hand rule in the latter case, a constant B_f,0 should be added to the toroidal field, which leaves the curl unaltered.

Model 1: a pure Lorentz wind. First, we consider a pure Contopoulos scenario. The disk has an Eddington factor of 10^-2, which is appropriate for CVs and YSOs. The poloidal field is set to zero. A temporally constant, toroidal field which drops off as r^-1 is chosen in the disk. Note that a pure toroidal field has no Alfvén point. Physically, it is still meaningful to assume that the flow is superalfvénic, in that an initially poloidal field got wound-up, above the Alfvén point, into the toroidal field. Therefore, both B_f and v_f are now fixed in the disk. The toroidal field strength is taken from model 2: the present model 1 serves mainly as a simple reference frame for the more intricate model 2. From the calculation we find that just above the disk, the vertical Lorentz force is a few times larger than the line force. Hence, the wind is essentially magnetically launched, i.e., a Contopoulos flow. Still, the mass loss rate is almost exactly that of a non-magnetic model discussed in the last section. The explanation is simple: since the flow starts superalfvénic, the CAK critical point is the only flow bottleneck, and determines the mass loss rate - whatever force launches the wind.

Model 2: magnetic eddies. Next, we chose a purely line driven models from the foregoing section as initial conditions (regularly spaced, dense streamers alternating with rarefied regions). A poloidal magnetic field is switched on slowly, assuming exponential growth. The field is initially axisymmetric, B_r=0. However, the field lines get tilted in time and get wound-up, creating B_z and B_f components. Typically, we find toroidal fields which, at the rim where the accretion disk touches the stellar surface, are 20 times stronger than the poloidal field. For a not-too big disk, B_f is still the dominant field component at the outer disk rim, only dropping off as r^-1. With B_f dominating, one expects that the mass loss rate is still the one of model 1. Instead, it is 20 times larger!

Figure 16: Eddies in the poloidal magnetic field above an accretion disk with small Eddington factor. The unit vector at the top right corresponds to 1 Gauss.

The reason is apparently the following. In model 2, the B_f bulge or wedge above the disk reaches to significantly larger heights than in model 1. The Lorentz force can assist then in overcoming the critical point, and the mass loss rate increases.

Why is the B_f bulge broader in model 2? This model shows a pronounced vortex sheet in the toroidal magnetic field and gas speed, as is shown in Fig. 11. The eddies are shed-off from the disk-star rim, and carry the B_f component to larger heights above the disk than in model 1. Indeed, the B_f bulge extends now to the top of the spinning eddies.

The plasma gun, and other open questions. What is the origin of the vortex sheet? It is known from plasma physics that toroidal fields are often unsteady and undergo periodic cycles of field build-up and unloading. This is called a `plasma gun'. We speculate that this is the origin of the present field eddies, too. Future work has to clarify this issue.

Other issues which shall be addressed in future work are: 1. Will a Contopoulos wind still form when a Blandford & Payne wind is launched right from model start? Instead of an initially axisymmetric field, tilt angles > 30 degrees with the disk normal shall here be assumed for the initial, poloidal field. 2. In present simulations the Alfvén surface `dives' periodically into the disk before it reaches the central star. This seems to be related to the shedding-off of magnetic eddies. 3. Riemann boundary conditions shall be applied along v_f, B_f characteristics, instead of the above, rather artificial boundary conditions averaging B_f over the disk.

As was already mentioned, the present models are first steps only. They seem, however, to keep the basic premise that magnetic fields can lead to increased mass loss rates in disk winds. These models shall also help to clarify whether magnetic fields can provide confinement of line driven quasar winds, in order to prevent over-ionization by central source radiation (deKool & Begelman 1995); and whether line driving can overcome problems with launching magnetocentrifugal winds from accretion disks. The latter problems were discussed in an important paper by Ogilvie & Livio (1998).

Outlook and Acknowledgments

With the book ``Foundations of Radiation Hydrodynamics'' by Mihalas & Mihalas recently being added to the Dover series of classical science texts, there seems no further need to motivate ``radiating fluids'' (a term coined by Mihalas & Mihalas). In astrophysical hydrodynamics, radiation and magnetic fields are of similar importance.

Radiation hydrodynamics splits into two branches, according to whether the flow is driven by photon absorption in the continuum or in spectral lines. The former class includes flows with Thomson scattering on free electrons and `dusty winds'. The latter class includes winds from hot stars and accretion disks, and is subject of the present writing. A simple, analytic approximation for the line force g_l from a point source of radiation is found by assuming a power law line list and validity of the Sobolev approximation, giving g_l ~ F (r^-1 dv/dr)^a, with radiative flux F and 0 < a < 1. Thus, g_l depends non-linearly on the velocity gradient dv/dr and matter density r, and represents a truly new, hydrodynamic force.

New waves (termed Abbott waves) and a new instability (termed de-shadowing instability) result from this new force, and are discussed in Chapters 3 and 2, respectively. Abbott waves could be responsible for shaping the velocity law of line driven winds, by causing a runaway towards the critical flow solution of Castor et al. (1975). The de-shadowing instability, on the other hand, is responsible for the formation of strong shocks and dense shells in the wind, and is probably the origin of observed X-ray emission from O stars; initial doubts after first Chandra X-ray line observations seem to be overcome. In Chapter 4, we gave a first glance at a rather new area in line driven wind hydrodynamics: accretion disk winds, and we discussed the interplay between radiative, Lorentz, and centrifugal forces.

What could be interesting topics for future work? With regard to the de-shadowing instability, all simulations so far are one-dimensional (with one puzzling exception), and two-dimensional hydrodynamic models are urgently missing. Line synthesis calculations using 2-D hydrodynamic models will allow to test our understanding of X-ray line formation and certain aspects of UV line variability. With regard to Abbott waves, the controversy whether they are an artefact of Sobolev theory is at present only partially overcome by post-Sobolev, numerical simulations. A strict derivation of the Green's function in presence of line scattering is missing. Finally, with regard to disk winds, some central questions are the increase in mass loss rates by combined magnetoradiative driving; magnetic flow confinement to prevent over-ionization by central source radation; and the amount of external disk viscosity caused by the wind.

With all these questions still being unanswered, line driven wind hydrodynamics appears to be a promising field for future research.

And now enough has been said (Aquinas, Sum. Theol., p. II.1, q. 114, a. 10). Almost. It is my pleasure to thank Janet Drew in London, Wolf-Rainer Hamann in Potsdam, Rolf-Peter Kudritzki in Hawaii, Colin Norman in Baltimore, Stan Owocki in Newark, Adi Pauldrach, Joachim Puls, and Christian Reile in Munich, and Isaac Shlosman in Lexington for many stimulating discussions, at blackboards and in restaurants.

Literature

Abbreviations of journals were adopted from the Astrophysics Data System.

Abbott D.C. 1980, ApJ, 242, 1183

Arav N. 1996, ApJ, 465, 617

Baade D. & Lucy L.B. 1987, A&A, 178, 213

Baum E., Hamann W.-R., Koesterke L., & Wessolowski U. 1992, A&A, 266, 402

Berghöfer T.W. & Schmitt J. 1994, Science, 265, 1689

Berghöfer T.W., Baade D., Schmitt J., et al. 1996, A&A, 306, 899

Bjorkman J.E. 1995, ApJ, 453, 369

Bjorkman J.E. & Cassinelli J.P. 1993, ApJ, 409, 429

Blandford R.D. & Payne D.G. 1982, MNRAS, 199, 883

Blondin J.M., Kallman T.R., Fryxell B.A., & Taam R.E. 1990, ApJ, 356, 591

Blondin J.M., Stevens I.R., & Kallman T.R. 1991, ApJ, 371, 684

Cannon C.J. 1973, ApJ, 185, 621

Carlberg R.G. 1980, ApJ, 241, 1131

Cassinelli J.P. 1979, ARA&A, 17, 275

Cassinelli J.P. & Hartmann L. 1977, ApJ, 212, 488

Cassinelli J.P. & Swank J.H. 1983, ApJ, 271, 681

Cassinelli J.P., Cohen D.H., MacFarlane J.J., Sanders W.T., & Welsh B.Y. 1994, ApJ, 421, 705

Cassinelli J.P., Miller N.A., Waldron W.L., MacFarlane J.J., & Cohen D.H. 2001, ApJ, in press

Castor J.I., Abbott D.C., & Klein R.I. 1975, ApJ, 195, 157 (CAK)

Chandrasekhar S. 1950, Radiative Transfer, Dover, New York

Chevalier R.A. & Imamura J.N. 1982, ApJ, 261, 543

Cohen D.H., Cassinelli J.P., & MacFarlane J.J. 1997, ApJ, 487, 867

Colella P. & Woodward R.P. 1984, J. Comp. Phys., 54, 174

Conti P.S. & Ebbets D. 1977, ApJ, 213, 438

Contopoulos J. 1995, ApJ, 450, 616

Cooper R.G. 1994, Ph. D. thesis, Univ. Delaware

Cooper R.G. & Owocki S.P. 1992, ASPC, 22, 281

Córdova F.A. & Mason, K.O. 1982, ApJ, 260, 716

Courant R. & Hilbert D. 1968, Methoden der Mathematischen Physik II, Springer, Berlin

Cox D.P. & Tucker W.H. 1969, ApJ, 157, 1157

Cranmer S.R. & Owocki S.P. 1996, ApJ, 462, 469

deJong J.A., Henrichs H.F., & Kaper L. 2001, A&A, 368, 601

deKool M. & Begelman M. 1995, ApJ, 455, 448

Deubner F.-L. 1973, IAU Symp., 56

Drew J., Proga D., & Stone J. 1998, MNRAS, 296, L6

Drew J., Proga D., & Oudmaijer R.D. 1999, in: Variable and Non-spherical Stellar Winds in Luminous Hot Stars, eds: Wolf B., Stahl O., & Fullerton A.W., Springer, Berlin, p. 140

Eisloffel J., Mundt R., Ray T.P., & Rodriguez L.F. 2000, in: Protostars and Planets IV, eds: Mannings V., Boss A.P., & Russell S.S., Univ. of Arizona Press, Tucson, p. 815

Elstner D. & Rüdiger G. 2000, A&A, 358, 612

Emmering R.T., Blandford R.D., & Shlosman I. 1992, ApJ, 385, 460

Eversberg T., Lépine S., & Moffat A.F. 1998, ApJ, 494, 799

Feldmeier A., Norman C., Pauldrach A., et al. 1997c, ASP Conf. Ser., 128, 258

Field G.B. 1965, ApJ, 142, 531

Friend D.B. & Abbott D.C. 1986, ApJ, 311, 701

Fullerton A.W., Massa D.L., Prinja R.K., Owocki S.P., & Cranmer S.R. 1997, ApJ, 327, 699

Gaetz T.J., Edgar R.J., & Chevalier R.A. 1988, ApJ, 329, 927

Gayley K.G. 1995, ApJ, 454, 410

Gayley K.G. & Owocki S.P. 1994, ApJ, 434, 684

Gayley K.G., Owocki S.P., & Cranmer S.R. 1995, ApJ, 442, 296

Gayley K.G. & Owocki S.P. 1995, ApJ, 446, 801

Gräfener G., Hamann W.-R., & Koesterke L. 2000, ASPC, 204, 215

Hamann W.-R. 1980, A&A, 84, 342

Hamann W.-R. 1981a, A&A, 93, 353

Hamann W.-R. 1981b, A&A, 100, 169

Hamann W.-R. 1985, A&A, 148, 364

Hamann W.-R. & Koesterke L. 1998a, A&A, 333, 251

Hamann W.-R. & Koesterke L. 1998b, A&A, 335, 1003

Hamann W.-R. & Koesterke L. 2000, A&A, 360, 647

Harnden F.R., Branduardi G., Elvis M., et al. 1979, ApJ, 234, L51

Hawley J.F., Smarr L.L., & Wilson J.R. 1984, ApJS, 55, 211

Heap S.R., Boggess A., Holm A. 1978, Nature, 275, 385

Hearn A.G. 1972, A&A, 19, 417

Hearn A.G. 1973, A&A, 23, 97

Hedstrom G.W. 1979, J. Comp. Phys., 30, 222

Henrichs H. 1986, PASP, 98, 48

Henrichs H., Gies D.R., Kaper L., et al. 1990, in: Evolution in Astrophysics, IUE Astronomy in the Era of New Space Missions, p. 401

Heyvaerts J. & Norman C.A. 1989, ApJ, 347, 1055

Hillier D.J. 1984, ApJ, 280, 744

Hillier D.J., Kudritzki R.P., Pauldrach A., et al. 1993, A&A, 276, 117

Hillier D.J. & Miller D.L. 1998, ApJ, 496, 407

Holzer T.E. 1977, J. Geophys. Res., 82, 23

Horne K. & Stiening R.F. 1985, MNRAS, 216, 933

Ignace R., Oskinova L.M., & Fullon C. 2000, MNRAS, 318, 214

Innes D.E., Giddings J.R., & Falle S.A. 1987, MNRAS, 226, 67

Kahn S.M., Leutenegger M.A., Cottam J., et al. 2001, A&A, 365, L312

Kaper L., Hammerschlag-Hensberge G., & Zuiderwijk E. 1994, A&A, 289, 846

Kaper L., Henrichs H., Nichols J., & Telting J. 1999, A&A, 344, 231

Klare G., Wolf B., Stahl O., et al. 1982, A&A, 113, 76

Knigge C., Woods J.A., & Drew J.E. 1995, MNRAS, 273, 225

Koesterke L. & Hamann W.-R. 1995, A&A, 299, 503

Koesterke L. & Hamann W.-R. 1997, A&A, 320, 91

Königl A. 1989, ApJ, 342 208

Königl A. & Pudritz R.E. 2000, in: Protostars and Planets IV, eds: Mannings V., Boss A.P., & Russell S.S., Univ. of Arizona Press, Tucson, p. 759

Koninx J.P. 1992, Ph. D. thesis, Univ. Utrecht

Krasnopolsky R., Li Z.Y., & Blandford R. 1999, ApJ, 526, 631

Krautter J., Vogt N., Klare G., et al. 1981, A&A, 102, 337

Krolik J.H. & Raymond J.C. 1985, ApJ, 298, 660

Krticka J. & Kubát J. 2000, A&A, 359, 983

Kudritzki R.P., Lennon D., & Puls J. 1995, in: Science with the VLT, eds: Walsh J. & Danziger I., Springer, Berlin, p. 246

Kudritzki R.P., Palsa R., Feldmeier A., Puls J., & Pauldrach A. 1996, in: Röntgenstrahlung from the Universe, eds: Zimmermann H., Trümper J., & Yorke H., MPE Report 263, Garching, p. 9

Kudritzki R.P., Puls J., Lennon D.J., et al. 1999, A&A, 350, 970

Kudritzki & Puls J. 2000, ARA&A, 38, 613

Lamb H. 1908, Proc. Lond. Math. Soc. (1), xix, 144

Lamb H. 1932, Hydrodynamics, Dover, New York

Lamers H.J.G.L.M., Cerruti-Sola M., & Perinotto M. 1987, ApJ, 314, 726

Langer N., Garcia-Segura G., & MacLow M.M. 1999, ApJ, 520, L49

Langer S.H., Chanmugam G., & Shaviv G. 1981, ApJ, 245, L23

Langer S.H., Chanmugam G., & Shaviv G. 1982, ApJ, 258, 289

Leitherer C., Schaerer D., Goldader J.D., et al. 1999, ApJS, 123, 3

Lépine S. & Moffat A.F. 1999, ApJ, 514, 909

Long K.S., Blair W.P., Davidsen A.F., et al. 1991, ApJ, 381, L25

Long K.S., Wade R.A., Blair W.P., Davidsen A.F., & Hubeny I. 1994, ApJ, 426, 704

Lucy L.B. 1975, Memoires Societe Royale des Sciences de Liege, 8, 359

Lucy L.B. 1982a, ApJ, 255, 278

Lucy L.B. 1982b, ApJ, 255, 286

Lucy L.B. 1983, ApJ, 274, 372

Lucy L.B. 1984, ApJ, 284, 351

Lucy L.B. 1998, in: Cyclical Variability in Stellar Winds, eds: Kaper L. & Fullerton A.W., Springer, Berlin, p. 16

Lucy L.B. & Solomon P.M. 1970, ApJ, 159, 879

Lucy L.B. & White R.L. 1980, ApJ, 241, 300

Lucy L.B. & Abbott D.C. 1993, ApJ, 405, 738

MacFarlane J.J. & Cassinelli J.P. 1989, ApJ, 347, 1090

MacFarlane J.J, Cohen D.H., & Wang P. 1994, ApJ, 437, 351

MacGregor K.B., Hartmann L., & Raymond J.C. 1979, ApJ, 231, 514

Massa D., Fullerton A.W., Nichols J.S., et al. 1995, ApJ, 452, L53

Meyer F. & Schmidt H. 1967, Zeitschrift für Astrophysik, 65, 274

Mihalas D. 1978, Stellar Atmospheres, Freeman, San Francisco

Mihalas D. & Mihalas B.W. 1984, Foundations of Radiation Hydrodynamics, Dover, New York

Miller K.A. & Stone J.M. 2000, ApJ, 534, 398

Milne E.A. 1926, MNRAS, 86 459

Moffat A.F. 1994, RvMA, 7, 51

Moffat A.F., Owocki S.P., Fullerton A.W., & St-Louis N. (eds) 1994, Instability and Variability of Hot-Star Winds, Ap&SS, 221

Murray N., Chiang J., Grossman S.A., & Voit G.M. 1995, ApJ, 451, 498

Norman M.L., Wilson J.R., & Barton R.T. 1980, ApJ, 239, 968

Norman M.L. & Winkler K.H. 1986, in: Astrophysical Radiation Hydrodynamics, eds: Winkler K.H. & Norman M.L., Reidel, Dordrecht, p. 187

Ogilvie G.I. & Livio M. 1998, ApJ, 499, 329

Oskinova L.M., Ignace R., Brown J.C., & Cassinelli J.P. 2001, A&A, in press

Ouyed R. & Pudritz R.E. 1997a, ApJ, 482, 712

Ouyed R. & Pudritz R.E. 1997b, ApJ, 484, 794

Owocki S.P. 1991, in: Stellar Atmospheres: Beyond Classical Models, eds: Crivellari L., Hubeny I., & Hummer D.G., Reidel, Dordrecht, p. 235

Owocki S.P. 1992, in: Atmospheres of Early-Type Stars, eds: Heber U. & Jeffery S., Springer, Berlin, p. 393

Owocki S.P. 1994, Ap&SS, 221, 3

Owocki S.P. 1999, in: Variable and Non-spherical Stellar Winds in Luminous Hot Stars, eds: Wolf B., Stahl O., & Fullerton A.W., Springer, Berlin, p. 294

Owocki S.P. & Rybicki G.B. 1984, ApJ, 284, 337

Owocki S.P. & Rybicki G.B. 1985, ApJ, 299, 265

Owocki S.P. & Rybicki G.B. 1986, ApJ, 309, 127

Owocki S.P., Castor J.I., & Rybicki G.B. 1988, ApJ, 335, 914 (OCR)

Owocki S.P. & Rybicki G.B. 1991, ApJ, 368, 261

Owocki S.P. & Zank G.P. 1991 ApJ, 368, 491

Owocki S.P., Cranmer S.R., & Blondin J.M. 1994, ApJ, 424, 887

Owocki S.P., Cranmer S.R., & Fullerton A.W. 1995, ApJ, 453, L37

Owocki S.P. & Puls J. 1996, ApJ, 462, 894

Owocki S.P., Cranmer S.R., & Gayley K.G. 1996, ApJ, 472, L115

Owocki S.P. & Puls J. 1999, ApJ, 510, 355

Owocki S.P. & Cohen D.H. 2001, ApJ, in press

Parker E.N. 1958, ApJ, 128, 664

Pauldrach A. & Puls J. 1990, A&A, 237, 409

Pauldrach A., Hoffmann T., & Lennon M. 2001, A&A, in press

Poe C.H., Owocki S.P., & Castor J.I. 1990, ApJ, 358 199

Proga D., Stone J.M., & Drew J. 1998, MNRAS, 295, 595

Proga D., Stone J.M., & Kallman T.R. 2000, ApJ, 543, 686

Pudritz R.E. & Norman C.A. 1983, ApJ, 274, 677

Pudritz R.E. & Norman C.A. 1986, ApJ, 301, 571

Puls J., Owocki S.P., & Fullerton A.W. 1993, A&A, 276, 117

Puls J., Feldmeier A., Springmann U., Owocki S.P., & Fullerton A.W. 1994, Ap&SS, 221, 409

Puls J., Kudritzki R.P., Herrero A., et al. 1996, A&A, 305, 171

Puls J., Springmann U., & Lennon M. 2000, A&AS, 141, 23

Rayleigh J.W. 1890, Scientific Papers, vol. III, Dover, New York 1964, p. 335

Reile C. & Gehren T. 1991, A&A, 242, 142

Rutten R.G., Dhillon V.S., Horne K., Kuulkers E., & Van Paradijs J. 1993, Nature, 362, 518

Rybicki G.B. & Hummer D.G. 1978, ApJ, 219, 654

Rybicki G.B., Owocki S.P., & Castor J.I. 1990, ApJ, 349, 274

Scharmer G.B. 1981, ApJ, 249, 720

Schmidt H. & Zirker J. 1963, ApJ, 138, 1310

Schulz N.S., Canizares C.R., Huenemoerder D., & Lee J.C. 2000, ApJ, 545, L135

Schulz W.D. 1964, J. Math. Phys., 5, 133

Sellmaier F., Yamamoto T., Pauldrach A., & Rubin R.H. 1996, A&A, 305, L37

Seward F.D., Forman W.R., Giacconi R., et al. 1979, ApJ, 234, L55

Shakura N.I. & Sunyaev R.A. 1973, A&A, 24, 337

Shlosman I., Vitello P., & Shaviv G. 1985, ApJ, 294, 96

Shlosman I. & Vitello P. 1993, ApJ, 409, 372

Simon M. & Axford W.I. 1966, Planet. Space Science, 14, 901

Springmann U. 1994, A&A, 289, 505

Springmann U. & Pauldrach A. 1992, A&A, 262, 515

Steidel C.C., Giavalisco M., Pettini M., Dickinson M., & Adelberger K.L. 1996, ApJ, 462, L17

Stein R.F. & Leibacher J. 1974, ARA&A, 12, 407

Stone J.M. & Norman M.L. 1992a, ApJS, 80, 753

Stone J.M. & Norman M.L. 1992b, ApJS, 80, 791

Stone J.M. & Norman M.L. 1994, ApJ, 433, 746

Thompson K.W. 1987, J. Comp. Phys., 68, 1

Thompson K.W. 1990, J. Comp. Phys., 89, 439

Turnshek D.A. 1984, ApJ, 280, 51

Uchida Y. & Shibata K. 1985, Publ. Astron. Soc. Japan, 37, 515

Ulrich R. 1970, ApJ, 162, 993

Ustyugova G.V., Koldoba A.V., Romanova M.M., Chechetkin V.M., & Lovelace R.V. 1999, ApJ, 516, 221

Van Leer B. 1977, J. Comp. Phys. 23, 267

Vitello P. & Shlosman I. 1988, ApJ, 327, 680

Vitello P. & Shlosman I. 1993, ApJ, 410, 815

Waldron W.L. & Cassinelli J.P. 2001, ApJ, 548, L45

Wax N. 1954, Selected Papers on Noise & Stochastic Processes, Dover, New York

Wegner M. 1999, Diploma thesis, Univ. Munich

Wessolowski U. 1996, in: Röntgenstrahlung from the Universe, eds: Zimmermann H., Trümper J., & Yorke H., MPE Report 263, Garching, p. 75

Weymann R.J., Turnshek D.A., & Christiansen W.A. 1985, in: Astrophysics of Active Galaxies and QSOs, ed: Miller J., OUP, Oxford, p. 333

Wolf B., Stahl O., & Fullerton A.W. (eds) 1999, Variable and Non-spherical Stellar Winds in Luminous Hot Stars, Springer, Berlin

Wu K., Chanmugam G., & Shaviv G. 1992, ApJ, 397, 232

Appended Papers

[1] Feldmeier A. 1995, A&A, 299, 523

[2] Feldmeier A., Anzer U., Börner G., Nagase F. 1996, A&A, 311, 793

[3] Feldmeier A., Kudritzki R., Palsa R., Pauldrach A., Puls J. 1997, A&A, 320, 899

[4] Feldmeier A., Puls J., Pauldrach A. 1997, A&A, 322, 878

[5] Feldmeier A. 1998, A&A, 332, 245

[6] Feldmeier A., Shlosman I. 1999, ApJ, 526, 344

[7] Feldmeier A., Shlosman I., Vitello P. 1999, ApJ, 526, 357

[8] Feldmeier A., Shlosman I. 2000, ApJ, 532, L125

[9] Feldmeier A., Shlosman I. 2002, ApJ, 564, 385

[10] Feldmeier A., Shlosman I., Hamann W-R. 2002, ApJ, 566, 392

Footnotes:

¹Only some key references are cited in the introduction.

²The numbers at the start of paragraphs correspond to the numbers in the figure on pages pageref and pageref.