2 General description of models

2.1 Conservation equations

Consider an inviscid flow along the radial direction r from a central star. In Eulerian form, the one-dimensional (1D) time-dependent equations for conservation of mass, momentum, and energy are:

$$% \frac{\partial\rho}{\partial t} +% \frac{1}{r^2}\frac{\partial(r^2 \rho v)}{\partial r}=0$$ (1)

$$\frac{\partial(\rho v)}{\partial t}+% \frac{1}{r^2}\frac{\par... ...% -\frac{\partial p}{\partial r} - \rho g_* + \rho g_{\rm rad}$$ (2)

$$\frac{\partial e}{\partial t} +% \frac{1}{r^2}\frac{\partial ... ...rac{p}{r^2}\frac{\partial (r^2 v)}{\partial r} - Q_{\rm rad} .$$ (3)

Here e is the internal energy density (in erg/cm³), g_* the effective stellar gravity (reduced by electron scattering), $g_{\rm rad}$ the acceleration due to line driving and $Q_{\rm rad}$ the power emitted by radiation per unit volume. The other symbols have their usual meaning. The set of equations is supplemented by the perfect gas law.

Since both the radiative driving and its associated instability are strongest in the radial direction (Rybicki et al. 1990), the assumption here of an effectively 1D, radial flow represents a reasonable first approximation. But more generally, it seems quite likely that the lateral coherence of structure generated by the radial, line-driven instabilities would be disrupted by other, lateral instabilities (e.g., Rayleigh-Taylor) that require a more complete multidimensional treatment. Unfortunately, the non-local integrations required for computation of the line-force become quite formidable to carry out in more than 1D, and as such, there have so far been only some quite limited explorations of the possible multidimensional nature of the structure arising from the line-driven instability (Owocki 1999). In this sense, the 1D simulations here can be thought as modelling the flow dynamics in some typical solid-angle patch of the overall spherical outflow, with the lateral coherence size of this patch yet to be determined (Dessart & Owocki 2001).

The conservation Eqs. (1)-(3) are solved using VH-1, a numerical code developed by J. M. Blondin and colleagues (Blondin, personal communication), by explicitly time-stepping the Lagrangian form of these equations and then remapping the evolved hydrodynamic variables back onto a fixed Eulerian grid after every time-step. In both the Lagrangian update and the remap, the Piecewise Parabolic Method (PPM; Colella & Woodward 1984) is used.

  
2.2 Line driving

The forces in the momentum equation are due to the pressure gradient, gravity and line driving. This last force is the most important in two ways: it is the force that drives the wind and its evaluation is what dominates the computing time of the simulation. In this paper, it is evaluated using the smooth source function (SSF) method (Owocki 1991). A more sophisticated approach (EISF, for "escape-integral source function") exists, but is computationally much more expensive and leads to structure that is qualitatively similar to the results of the cheaper SSF method (OP99).

The force from all lines together is obtained as the integral of the single-line force over a spectrum-averaged line-number distribution

$$\frac{{\rm d}N (\kappa)}{\rm d\kappa} = \frac{1}{\kappa_0} \... ...appa_0}\right)^{\alpha-2} {\rm e}^{-\kappa/\kappa_{\rm max}}.$$ (4)

The normalisation constant $\kappa_0$ is related to the strength of the strongest line. The exponent $\alpha$ can be interpreted as the ratio of the force from optically thick lines to the total line force. The cut-off parameter $\kappa _{\rm max}$limits the maximum line strength (see Sect. 7). Although the CAK line-force assumes $\kappa_{\rm max} \rightarrow \infty$, its implicit assumption of a continuous power-law distribution breaks down for the discrete number of strong lines, i.e. when $\kappa \,{\rm d} N / {\rm d}\kappa \approx 1$. More realistically, one should thus generally impose an opacity cut-off with $\kappa_{\rm max} \approx \kappa_{0}$.

For sufficiently thick wids, the above parametrisation leads to a force that is proportional to $\tau_0^{-\alpha}$, where $\tau_0$ is the optical thickness for a single line of opacity $\kappa_0$. For winds that are optically thin in all lines the force no longer increases with decreasing optical depth (Abbott 1982).

For the above line distribution, OP96 derive from first principles the expressions for the line force in the SSF approximation, giving a full discussion of the underlying assumptions.

2.3 Cooling and heating

In the energy equation, the two terms on the right-hand-side represent adiabatic and radiative cooling. The power radiated per unit volume can be written as

$$Q_{\rm rad} = n_{\rm e} N_{\rm H} \Lambda(T),$$ (5)

where $n_{\rm e}$ and $N_{\rm H}$ respectively represent the number density of electrons and hydrogen (atoms and ions), and the cooling coefficient $\Lambda$ describes the temperature dependence of the cooling (Raymond et al. 1976).

The combined effect of this cooling is balanced by photoionisation heating from the star's UV radiation, which in a relatively smooth, steady outflow has the general effect of thus preventing the wind from cooling much below the stellar effective temperature (e.g., Drew 1989). In this paper, we mimic this effect of photoionisation heating by simply requiring that the wind temperature never fall below a floor value, typically given by the stellar effective temperature. A more complete treatment of such heating in the kind of highly structured, time-dependent flows modelled here is beyond the current scope, and so will be deferred to future study.

Indeed in the structured models here, adiabatic compression and shock conversion of flow kinetic energy often lead to regions of strong compressive heating, which are then subsequently cooled by the radiative emission, as modelled by the cooling function $\Lambda$. Feldmeier (1995) describes how numerical hydrodynamical flows with such cooling are subject to various effects that cause the cooling zone to effectively collapse, thus leading to an underresolution that causes the amount of hot gas to be substantially underestimated. In an effort to identify sufficient hot material to explain observed levels of hot-star X-ray emission, Feldmeier (1995) artificially limited the slope of the cooling function in a way that partially mitigates this underresolution. In the present context, resolving the hot gas is less vital, and we have opted to retain a tabulated cooling function without such modifications. Moreover, our use here of a constant, relatively fine spatial grid means that such underresolution is less severe (see Fig. 1) than in previous models with grid spacing that increases with radius.

2.4 Boundary and initial conditions

The simulation volume extends from $r= 1\;R_*$, the photospheric radius where the optical depth for electron scattering is approximately one, to an outer radius $r=R_{\rm max}=101\;R_*$. At the inner boundary we fix the density $\rho$ and the temperature T but allow the velocity v to float, with the requirement that its second derivative be zero. At the outer boundary, both $\rho$ and v may vary, but must have zero second derivatives, while T is taken constant. Unlike some previous simulations, we do not explicitly perturb the base of the stellar wind by a sound wave or turbulence. Nonetheless, we find the flow in the outer wind develops a persistent, "intrinsic'' variability. As discussed by OP99, such self-excitation of intrinsic variability in the absence of any explicit perturbations has its physical origin in the backscattering of radiation from outer wind structure, which thus seeds small variations at the wind base that subsequently become amplified to maintain the strong structure of the outer wind. The exact reason why self-excited structure is stochastic in nature, and not e.g. quasi-periodic, is not fully understood and requires further study.

Each simulation starts from a smooth wind initial condition computed from a CAK/Sobolev approximation for the line-force. It is then run using the SSF line-force for 2 Msec to allow the response to the initial condition to die away (the time needed to cross the simulation volume at 2000 kms^-1 is ${\sim}0.7$ Msec). The model is then evolved for detailed study over a further 500 ksec. Time-averaged quantities needed for the statistical description of structure (Sect. 3) are calculated by summing the relevant variables at every time step over this 500 ksec interval.

A further important detail is that the time evolution is computed with a fixed time-step, rather than using a time-step that varies to satisfy a fixed Courant number (see, e.g., Laney 1998). The latter approach can introduce an additional, artificial feedback between outer and inner wind structure. (See, e.g., discussion surrounding Eq. (31) of Poe et al. 1990.) Our choice of a fixed time-step of 5 s satisfies the Courant condition with a Courant number never exceeding 0.5.