© ESO 2018

1 Introduction

For massive, hot stars of spectral types OBA, scattering and absorption in spectral lines transfer momentum from the star’s intense radiation field to the plasma, and so provide the force necessary to overcome gravity and drive a strong stellar wind outflow (see Puls et al. 2008, for an extensive review). The first quantitative description of such line-driving was given in the seminal paper by Castor et al. (1975; hereafter CAK). Like many wind models to date, CAK used the so-called Sobolev approximation (Sobolev 1960) to compute the radiative acceleration. This assumes that hydrodynamic flow quantities¹ are constant over a few Sobolev lengths ℓ_Sob = v_th∕(dv_n∕dn) (for ion thermal speed v_th and projected velocity gradient dv_n∕dn along a coordinate direction $\widehat{n}$), allowing then a local treatment of the line radiative transfer.

This Sobolev approach ignores the strong line-deshadowing instability (LDI) that occurs on scales near and below the Sobolev length (Owocki & Rybicki 1984); numerical radiation-hydrodynamic modeling of the nonlinear evolution of the LDI shows that the time-dependent wind develops a very inhomogeneous, “clumped” structure (Owocki et al. 1988; Feldmeier et al. 1997; Dessart & Owocki 2003, 2005; Sundqvist & Owocki 2013, 2015). These clumpy LDI models provide a natural explanation for a number of observed phenomena in OB stars, such as the soft X-ray emission and broad X-ray lines observed by orbiting telescopes like Chandra and XMM-Newton (Feldmeier et al. 1997; Berghoefer et al. 1997; Güdel & Nazé 2009; Cohen et al. 2010; Martínez-Núñez et al. 2017), the extended regions of zero residual flux typically seen in saturated UV resonance lines (Lucy 1983; Puls et al. 1993; Sundqvist et al. 2010), and the migrating spectral subpeaks superimposed on broad optical recombination lines (Eversberg et al. 1998; Dessart & Owocki 2005b; Lépine & Moffat 2008).

However, a severe limitation of most of the above-mentioned models is their assumed spherical symmetry. The fact that most LDI simulations in the past have been limited to 1D is mainly a consequence of the computational cost associated with carrying out the nonlocal integrals needed to compute the radiation acceleration at each simulation time step, while simultaneously resolving length scales below ℓ_Sob. Specifically, following the general escape-integral methods developed by Owocki & Puls (1996), approximately n_x ≈ 3v_∞∕v_th ≈ 1000 discrete frequency points are typically needed to properly resolve line profiles and model the expanding flow. In 2D or 3D, a proper treatment of the multi-dimensional wind further requires integrations along a set of oblique rays in order to compute the radiative force in the radial and lateral directions. A major issue then becomes the misalignment of nonradial rays with the discrete numerical grid (i.e., that oblique ray-integrations from any given point in the mesh in general do not intersect at any other point), requiring that all integrations be repeated for each grid node (then also involving complex interpolation schemes to trace the rays).

As an explicit example (see also Dessart & Owocki 2005), for a 2D grid of n_r radial and n_ϕ azimuthal points, n_rn_ϕ integrations on the order of n_rn_x operations are needed for every considered ray; this gives an overall scaling $n_{\mathrm{\text{ray}}}{n}_x{n}_r^2{n}_{\varphi }$, implying for a typical case of n_ray ≈ 5, n_ϕ ≈ 100, and n_x ≈ n_r ≈ 1000 on the order of 10¹¹⁻¹² operations to evaluate the radiative force. Moreover, such a calculation has to be carried out at each time step of the hydrodynamical simulation, which for a typical courant time ~ 5 s in a hot-starwind outflow, and a total simulation time of, e.g., ~ 50 dynamical time scales t_dyn = R_*∕v_∞~ 10 ks, requires approximately ~10⁵ repeated evaluations of the radiative force. This simple example thus illustrates quite vividly the rather daunting task of constructing multi-dimensional LDI wind models.

Nonetheless, a few previous attempts have been performed. Dessart & Owocki (2003) carried out 2D hydro+1D radiation simulations by simply focusing only on the line force from a single radial ray, thus ignoring lateral influences. This led then to an extensive breakup of the spherical shells seen in 1D simulations, resulting in lateral incoherence all the way down to the grid scale. However, these simulations ignore the lateral component of the diffuse radiative force, which linear stability analysis (Rybicki et al. 1990) shows could lead to damping of velocity variations at scales below the lateral Sobolev length ℓ_Sob = rv_th∕v, and as such to more lateral coherence than seen in the single-ray 2D simulations. Dessart & Owocki (2005) made a first attempt to include oblique rays, by using a special, restricted numerical grid in a 2D plane that forced three rays to always intersect at the discrete mesh points (Owocki 1999). But while these simulations did seem to suggest a somewhat larger lateral coherence than comparable one-ray models, the inherent limitations of the method (e.g., in resolving the proper lateral scales) provided uncertain results (Dessart & Owocki 2005).

This paper introduces a pseudo-planar modeling approach for a multi-dimensional wind subject to the LDI. In this box-in-a-wind method, all sphericity effects of the expanding flow are included in a radial direction r, but some curvature terms are ignored in the lateral direction(s). As we discuss in Sect. 2 (and detail in Appendix A), for a 2D simulation in the r, y plane, this allows us to consider five long-characteristic rays with a computational cost-scaling 3n_xn_rn_y, thus reducing the general scaling above with a factor ~ n_r = 1000 for our standard setup. Using this method, in Sect. 3 we examine the resulting 2D clumpy wind structure in much greater detail than possible before, and in Sect. 4 we discuss the results, we make comparisons with other simulation test runs, and we outline future work.

2 Modeling

2.1 Hydrodynamics

The simulations here use the numerical PPM (Colella & Woodward 1984) hydrodynamics code vh-1² to evolve the conservation equations of mass and momentum for a 2D, isothermal line-driven stellar wind outflow. A key point of this paper is that we keep all sphericity effects of an expanding outflow in the radial direction, but we neglect some curvature terms in the lateral direction(s) (see Appendix A for details). Preserving all properties of a spherical outflow, this pseudo-planar, box-in-a-wind approach allows us to resolve laterally the relevant clump length scales, and to implement nonradial rays for the radiative line-driving in a time efficient way (see further below and Appendix A).

All the presented results adopt the same stellar and wind parameters as in Sundqvist & Owocki (2013, 2015), given here in Table 1, which are typical for an O-star in the Galaxy. The standard setup uses a spatial grid with 1000 discrete radial (r) mesh-points between R_*≤ r ≤ 2R_* and 100 lateral (y) ones that cover in total 0.1R_*. As such, the grid is uniform and has a constant step size Δ = 0.001R_*; a small Δ is required to resolve both the subsonic wind base with effective scale height $H=a^2{R}_{*}^2∕\left(G{M}_{*}\left(1-{\mathrm{\Gamma }}_{\mathrm{\text{e}}}\right)\right)\approx 0.002R_{*}>\mathrm{\Delta }$ and the resulting small-scale 2D clump structures in the supersonic wind (the focus of this paper). Each simulation evolves from a smooth, CAK-like initial condition, computed by relaxing to a steady state a 1D spherically symmetric time-dependent simulation that uses a CAK/Sobolev form for the line force. To prevent artificial structures due to numerical truncation errors, we use an evolution time step that is the minimum of a fixed 2.5 s and a variable 1/3 of the courant time (see discussion in Poe et al. 1990). As in previous work, the lower boundary at the assumed stellar surface fixes the density to a value ~5–10 times that at the sonic point. Moreover, since we are interested in structures that are considerably smaller than the computational box, the lateral boundaries are simply treated as periodic.

2.2 Radiative driving

The central challenge in these simulations is to compute the 2D radiation line force in a highly structured, time-dependent wind with a non-monotonic velocity. This requires nonlocal integrations of the line-transport within each time step of the simulation in order to capture the instability near and below the Sobolev length. To meet this objective, we develop here a multi-dimensional pseudo-planar extension of the smooth source function (SSF; Owocki 1991) method described extensively in Owocki & Puls (1996; see also Sundqvist & Owocki 2013). Appendix A describes in detail this 2D-SSFformulation; below follows a summary of the key features.

Our pseudo-planar 2D-SSF approach allows us to follow the nonlinear evolution of the strong, intrinsic LDI in the radial direction while simultaneously accounting for the potentially stabilizing effect of the scattered, diffuse radiation field, both in the radial and lateral directions (Lucy 1984; Owocki & Rybicki 1985; Rybicki et al. 1990). Moreover, SSF assumes that the line-strength number distribution is given by an exponentially truncated power law. In this formalism, α is the standard CAK power-law index, which can be physically interpreted as the ratio of the line force due to optically thick lines to the total line force; $\bar{Q}$ is a line-strength normalization constant, which can be interpreted as the ratio of the total line force to the electron scattering force when all lines are optically thin; Q_max is the maximum line-strength cutoff³. For typical O-star conditions at solar metallicity, $Q_{\mathrm{\text{max}}}\approx \bar{Q}\approx 2000$ (Gayley 1995; Puls et al. 2000). In practice, keeping the nonlinear amplitude of the instability from exceeding the limitations of the numerical scheme requires a significantly smaller cutoff (Owocki et al. 1988; Sundqvist & Owocki 2013).

As noted in Sect. 1, including oblique rays in a multi-dimensional outflow presents severe computational challenges, largely due to the general misalignment of the rays with the nodes of the numerical grid. While earlier attempts of 2D LDI simulations have either used a 2D hydro+1D radiation approach (Dessart & Owocki 2003) or experimented with a restricted special radial grid setup (Dessart & Owocki 2005), the pseudo-planar method introduced here largely circumvents these issues of grid misalignment. Namely, while radial ray-integrations are here calculated as they were in the original SSF method, for oblique rays both the azimuthal radiation angle ϕ and the ray’s radial directional cosine $\mu \equiv cos\theta =\widehat{r}\cdot \widehat{n}$ become constant throughout the computational domain. To this end, we apply a set of five rays with $\mu ,\varphi =\left(1,1∕\sqrt{3},1∕\sqrt{3},1∕\sqrt{3},1∕\sqrt{3}\right),\left(0,\pi ∕4,-\pi ∕4,3\pi ∕4,-3\pi ∕4\right)$ (see Fig. 1 and Appendix A for a detailed explanation). In addition to the (trivial) radial ray, this thus considers four oblique rays that are also pointing up/down with respect to the 2D equatorial plane in which the hydrodynamical calculations are carried out (in order to avoid certain 2D “flat-land” radiation effects; see Gayley & Owocki 2000). For our assumed grid then, with constant spacings in the radial and lateral directions, information can be used for all grid nodes when the ray integration for a given (μ, ϕ) pair has been performed only once over R_*≤ r ≤ 2R_* for each of the lateral grid points. This means that the solid angle integrations required to compute the line force in the radial and lateral directions then can be performed without the need of any further ray integrations. In addition, because of the symmetry of rays pointing up/down from the equatorial plane, we only have to explicitly carry out the integrations for three of our five angles. With respect to the general situation, this means we have effectively reduced the number of required “long-characteristic” ray integrations at each time step with a factor of ~n_r ( = 10³ for our standard setup here).

Another attractive feature of this pseudo-planar model is that it preserves all properties for a 1D purely radial outflow. As detailed in Appendix A, this is achieved by preserving the general scaling of the flux with radius for a spherical outflow by including a sink term for the density to mimic spherical divergence and by including in the force equations terms to account for stellar rotation along the lateral axis y. As such, our approach allows for easy testing and benchmarking, and we have verified that a simulation run with n_y = 1 and integration weights for all oblique rays set to zero indeed gives the same results as a “normal” 1D spherical radial-ray SSF simulation. However, since such radial models are also subject to the global wind instability associated with nodal topology (Poe et al. 1990; Sundqvist & Owocki 2015), they exhibit clumpy structure all the way down to the lower boundary (Sundqvist & Owocki 2013, 2015). While there are strong indications that clumping in hot-star winds indeed extends to very near-photospheric layers (e.g., Cohen et al. 2014), in these first 2D simulations we nonetheless opt to stabilize the wind base by introducing a small radial increase in $\bar{Q}$ between R_* < r < 1.5R_*. This allows us to study the emerging clump formation and structure in a somewhat more controlled environment with respect to simulations that lie on the nodal topology branch (see Sect. 4).

Fig. 1 Sketch illustrating the basic idea of the pseudo-planar, box-in-a-wind approach used in this paper. The upper left illustrates the general situation of non-alignment between oblique rays and the numerical grid points. The lower panel then shows how we create a pseudo-planar box in the wind by cutting out a small but representative fraction of the wind volume. For illustration purposes, it shows projections onto the equatorial plane of rays in the prograde (blue), retrograde (red), and radial (black) directions for two lateral periods of a simple case with just ny = 2 zones in lateral direction y. The right panel then illustrates how extension out of the equatorial plane involves a total of five rays: one radial plus two oblique pairs that extend up/down from the plane. See Appendix A for a detailed explanation and for further illustrations of the assumed ray geometry.

3 Simulation results

3.1 General properties

Figure 2 illustrates a key result of our simulations, namely the spatial and temporal variation in log ρ relative to the initial, smooth CAK steady state. The figure showsclearly how a radial shell structure first develops, but then quickly breaks up into laterally complex density variations. The upper panel displays snapshots during the first 100 ks of the simulation, illustrating how already after a few dynamical flow times t_dyn ≈ R_*∕⟨v_max⟩≈ 11 ks the characteristic shells, seen in all 1D LDI simulations, break up in what initially seem to resemble Rayleigh–Taylor structures. The lower panel then shows how over time the structures eventually develop into a complex but statistically quite steadyflow, characterized now by localized density enhancements (clumps) of very small spatial scales embedded in larger regions of much lower density.

Figure 3 zooms in on the same log density in a small 0.1R_* square box over a short time sequence long after the initial condition. This illustrates in greater detail the quite complex 2D density structure showing a range of scales, as well as high-densityclumps with different shapes. The figure also demonstrates that although the structures are small, they are clearly resolved by our numerical grid.

Figure 4 displays temporal and spatial variations in radial velocity, illustrating essentially the same kind of outer wind shock structure and high-velocity streams as in corresponding 1D simulations; however, the velocity now also exhibits extensive lateral variation, reflecting again the breakup of 1D shells into small-scale 2D clumps.

Figure 5 emphasizes some similarities between these 2D simulations and the corresponding 1D simulations by showing a radial cut through the simulation box at a time snapshot (again taken long after the simulation has developed into a statistically steady flow). The figure demonstrates how such radial cuts indeed still show the characteristic structure of the nonlinear growth of the LDI, namely high-speed rarefactions that steepen into strong shocks and wind plasma compressed into spatially narrow clumps separated by rather large regions of rarefied gas. There are some differences though; in addition to the lateral breakup of shells discussed above, another key distinction between 1D and 2D simulations is that the radial density variations are a bit lower in the latter, which occurs because of the lateral “filling in” of radial rarefactions (see also Dessart & Owocki 2003) and is discussed further in the following section.

Fig. 2 Spatial and temporal variations of log density relative tothe initial, smooth CAK steady state at t = 0, with color ranging from densities a decade below the t = 0 value (blue) to a decade above (red). The vertical variation extends from the subsonic wind base at the stellar surface R* to a height of one R* above. For clarity, the lateral variation is displayed over twice the horizontal box length 0.1R* . The upper row shows time evolution over the initial 100 ks after the CAK initial condition in steps of 10 ks; the bottom row uses the same step size of 10 ks to show the evolution between 300 and 400 ks, long after the initial condition has developed into a statistically steady turbulent flow.

Fig. 3 As in Fig. 2, spatial and temporal variations of log density relative to the initial, smooth CAK steady state at t = 0 are shown, with color ranging from densities a decade below the t = 0 value (blue) to a decade above (red). Here the vertical variation only extends between 1.9R* and 2.0R* and the lateral variation is displayed over one horizontal box of 0.1R* ; there are thus 100 × 100 discrete mesh points in each of the displayed squares. From left to right are shown a 2 ks time evolution long after the initial condition, in steps of 0.5 ks.

Fig. 4 Spatial and temporal variations of radial velocity vrad , with color ranging from 0 (blue) to 2000 kms−1 (red). As in Fig. 2, the vertical variation extends from the subsonic wind base at the stellar surface R* to a height of one R* above, and the lateral variation is displayed over twice the horizontal box length 0.1R* . The frames from left to right show the time evolution of vrad over 400 ksafter the CAK initial condition, in steps of 50 ks.

Fig. 5 Radial cuts through the 2D simulation box of density ρ [g cm−3 ] (left) and radial velocity vrad [cm s−1] (right). Thered curves are taken at a snapshot long after the simulation has developed into a statistically quite steady flow; the black curves compare this to average values.

3.2 Statistical properties

Figure 6 summarizes some statistical results of the simulations. All averaging here starts at t = 250 ks in order to separate out any dependence on the initial conditions and the adjustment to a new radiative force balance. The upper left panel in Fig. 6 shows the clumping factor $$f_{\mathrm{\text{cl}}}=\frac{⟨{\rho }^2⟩}{{⟨\rho ⟩}^2},$$(1)

where anglebrackets denote averaging both laterally and over time in order to separate out the f_cl primary dependence on radius. The plot illustrates how the lateral filling in of radially compressed gas (see above) decreases the quantitative clumping factor significantly in a 2D simulation compared to earlier 1D models where f_cl ≳10 (e.g., Sundqvist & Owocki 2013); this is also consistent with the previous 2D results by Dessart & Owocki (2003). We note, however, that the actual values of f_cl in our 2D simulation are likely somewhat underestimated, due to our choice of stabilizing the wind base against instability caused by nodal topology (see previous section). As discussed extensively by Sundqvist & Owocki (2013), in these near-photospheric layers the quantitative clumping factor is very sensitive to the choices made for the calculation of the radiative acceleration, and to any variability that may be assumed for the photospheric lower boundary. Regardless of these caveats, the basic qualitative result here that 2D simulations yield relatively lower values of f_cl than comparable 1D simulations is quite robust.

The upper right panel of Fig. 6 then shows the time-dependent mass-loss rate: $$Ṁ\equiv 4\pi r^2\rho v_{\mathrm{\text{rad}}}.$$(2)

The black line in this plot shows a simple lateral average of the mass-flux escaping the outermost radial grid-point at a specifictime. However, since our simulation box only covers 0.1R_*, this average very likely overestimates the time-dependent mass loss significantly. To compensate for this, the red curve in the plot instead uses an average over all grid points at r ≥ 1.5R_* at a specific time, which approximates averaging over a full stellar surface $4\pi {\left(2R_{*}\right)}^2\approx 50R_{*}^2.$ As expected, this curve shows a drastically lower temporal variation of Ṁ , despite the highly time-dependent flow. This is consistent with, e.g., decade-long observations of spectral lines in O-stars, which typically indicate that time variations in the mass-loss rate of such stars are small.

To estimate typical clump length scales, the lower left panel of Fig. 6 plots a density autocorrelation length $$f_{\mathrm{\text{c}}}\left(\mathrm{\Delta }\right)=\sum _{\mathrm{\text{time}}}\sum _i\left({\rho }_i-⟨\rho ⟩\right)\left({\rho }_{i-\mathrm{\Delta }}-⟨\rho ⟩\right),$$(3)

where ⟨ρ⟩ averages laterally and over time. A lateral correlation length is calculated at each of the Δ = 0 − 99 lateral mesh points and is normalized to its Δ = 0 value. The figure then plots an average of this lateral correlation length between r∕R_* = 1.9 − 2.0 (black curve), as well as a radial correlation length (red curve) defined analogously. The lateral and radial density correlation lengths are very similar, and as such illustrate how a statistical ensemble of clumps is quite isotropic in these simulations. This does not imply that any given clump is isotropic (see Fig. 3), but rather that, on average, the well-developed density variations in the simulations do not have a strong preferred direction.

The Gaussian fit plotted in the blue dashed curve provides an estimate of the autocorrelation length in terms of the Gaussian FWHM ≈ 0.01R_*. Such small characteristic scales agree well with the theoretical expectation (see Sect. 1) that the critical length scale for these clumpy wind simulations is on the order of the Sobolev length ℓ_Sob, which for the lateral direction at 2R_* is ℓ_Sob∕R_* = 2v_th∕v ≈ 0.01. Identifying this as a typical clump length scale ℓ_cl, we can further make a simple estimate of the typical clump mass ${\ell }_{\mathrm{\text{cl}}}^3{\rho }_{\mathrm{\text{cl}}}\approx 10^{-6}{R}_{*}^{3}7\times 10^{-14}gc{m}^{-3}\approx 10^{17}g$, where the estimated clump density here simply reads off the output of the simulations (see Fig. 5). More generally, such a clump mass-estimate may be obtained using the Sobolev length and mass conservation $$m_{\mathrm{\text{cl}}}\approx {\ell }_{\mathrm{\text{Sob}}}^3{\rho }_{\mathrm{\text{cl}}}\approx \frac{v_{\mathrm{\text{th}}}^3Ṁf_{\mathrm{\text{cl}}}r}{v^{4}4\pi },$$(4)

which for the 2D simulation analyzed here indeed gives m_cl ≈ 10¹⁷ g for typical values at 2R_*. Quite generally, Eq. (4) shows explicitly how rather low clump masses are expected to emerge from the LDI.

Finally, the lower right panel in Fig. 6 plots the radial and lateral velocity dispersions $$v_{\mathrm{\text{disp}}}=\sqrt{⟨v^2⟩-{⟨v⟩}^2},$$(5)

where averages are constructed as they were for the clumping factor above. These plots show how, as expected (see also Dessart & Owocki 2003), the lateral velocity dispersion is on the order of the isothermal sound speed, whereas the radial dispersion is much higher and is expected to rise above several hundreds km s⁻¹ in the outer wind.

Fig. 6 Selected statistical properties of the 2D simulation (see text). The upper left panel plots the clumping factor fcl ; the upper right panel shows the time-dependent mass-loss rate, Ṁ [M⊙ yr−1] vs. s, computed in two different ways for the red and black curves(see text); the lower left panel displays lateral (black) andradial (red) density correlation lengths, and the corresponding Gaussian fit (blue, dashed); the lower right panel then finally plots radial (left) and lateral (right) velocity dispersions.

4 Summary and future work

4.1 Summary

We have introduced a pseudo-planar, box-in-a-wind approach suitable for carrying out radiation-hydrodynamical simulations in situations where the computation of the radiative acceleration is challenging and time consuming. The method is used here to simulate the2D nonlinear evolution of the strong LDI that causes clumping in the stellar winds from hot, massive stars. Accounting fully for both the direct and diffuse radiation components in the calculations of both the radial and lateral radiative accelerations, we examine in detail the small-scale clumpy wind structure resulting from our simulations.

Overall,the 2D simulations show that the LDI first manifests itself by mimicking the typical shell structure seen in 1D simulations, but these shells then quickly break up because of basic hydrodynamic instabilities (e.g., Rayleigh–Taylor) and the influence of the oblique radiation rays. This results in a quite complex 2D density and velocity structure, characterized by small-scale density clumps embedded in larger regions of fast and rarefied gas.

While inspection of radial cuts through the 2D simulation box confirms that the typical radial structure of the LDI is intact, quantitatively the lateral filling in of gas leads to lower values of the clumping factor than for corresponding 1D models. A correlation-length analysis also shows that, statistically, density variations in the well-developed wind are quite isotropic; identifying then the computed autocorrelation length with a typical clump size gives ℓ_cl∕R_*~ 0.01 at 2R_*, and thus also quite low typical clump masses m_cl ~ 10¹⁷ g. This agrees well with the theoretical expectation that the important length scale for LDI generated wind structure is on the order of the Sobolev length ℓ_Sob.

4.2 Influence of rotation and topology

As noted in Sects. 2 and 3.1, the level of structure in near photospheric layers is likely underestimated in the simulation analyzed above because we chose to stabilize the wind base. To demonstrate this further, Fig. 7 shows a test run with the same 2D setup as before, but now introducing stellar rotation with a fixed v_rot = 300 km s⁻¹ at the surface, and an initial condition set by steady-state angular momentum conservation, v_y (r) = v_rotR_*∕r. The figure shows that once the simulation has adjusted to its new force conditions, radial streaks of high density now appear at the surface; in other test runs, we have found that such structures are typical for simulations with an unstable base and nodal topology. The radial streaks in this rotating model migrate along with the surface rotation; embedded in the larger scale structures are the typical small-scale clumps discussed previously. As speculated previously in Sundqvist & Owocki (2015), these tentative first results thus suggest that rotating LDI models may quite naturally lead to the type of combined large- and small-scale structures needed to explain in parallel various observed phenomena in hot-star winds, such as discrete absorption components (DACs; Kaper et al. 1999) and small-scale wind clumping (Eversberg et al. 1998). Future work will examine in detail the connections between these rotating LDI models and the presence of various types of wind substructure.

Fig. 7 Spatial and temporal variations of log density, radial velocity, and lateral velocity for a modelwith stellar rotation at the surface vy = 300 km s−1 (see text), with color-coding as in Figs. 2 and 4. The vertical variation in this simulation extends only from 1.0 to 1.5 R* , but the lateral variation is displayed as before over twice the horizontal box length 0.1R* . From left to right are shown the time evolution over 350 ks after the CAK initial condition, in steps of 50 ks.

4.3 Future work

The simulations presented in this paper also lead naturally to a number of follow-up investigations; already in the pipeline are the development of a formalism for characterizing porosity effects in turbulent media (Owocki & Sundqvist 2018) and a study on the influence of the clumpy wind on the accretion properties of an orbiting neutron star in a high-mass X-ray binary system (el-Mellah et al. 2018). More directly related to this paper, in addition to further analyzing the effects of rotation and topology, we also plan to extend the current simulations to 3D and to higher wind radii, and also to develop a more general radiative transfer scheme (allowing for an arbitrary number of rays) for the computation of the line acceleration within a pseudo-planar box-in-a-wind.

Uncited reference

Dessart and Owocki (2005a).

Acknowledgments

This work was supported in part by SAO Chandra grant TM3-14001A awarded to the University of Delaware, and in part by the visiting professor scholarship ZKD1332-00-D01 for S.P.O. from KU Leuven. S.P.O. acknowledges sabbatical leave support from the University of Delaware, and we also thank John Castor for helpful discussions on long-characteristic methods. We finally thank the referee for the useful comments on the paper.

Appendix A 2D pseudo-planar line force

Our development here of a 2D vector form for the line acceleration follows a direct generalization of the 1D-SSF method detailed in Owocki & Puls (1996; hereafter OP96), as further developed in Sundqvist and Owocki (2015; hereafter SO15). As discussed in OP96 (see their Eq. (3)), a key step is to compute efficiently the profile-weighted line optical depth between two wind locations along a ray coordinate z, $$\mathrm{\Delta t}\left(x,z_1,z_2\right)={\int }_{z_1}^{z_2}{\kappa }_{\mathrm{\text{o}}}\rho \left(z\right)\varphi \left[x-u_z\left(z\right)\right]\mathrm{\text{d}}z,$$(A.1)

where κ_o is a spatially constant line opacity normalization⁴, and u_z ≡ẑ ⋅ u is the local z-projection of the vector flow velocity normalized to the ion thermal speed, u ≡v∕v_th. The line profile function is taken here to have a normalized Gaussian form, $\varphi \left(x\right)={\mathrm{\text{e}}}^{-x^2}∕\sqrt{\pi }$, with x = (ν − ν_o)∕Δν_D the observer-frame frequency displacement from line center in thermal Doppler units Δν_D = ν_ov_th∕c.

In 1D spherically symmetric models where variables only depend on the radius r, the ray direction is defined in terms of the local r and a stellar impact parameter p, with $\left|z\right|=\sqrt{r^2-p^2}$, and the signtaken to be positive (negative) in the forward (backward) hemisphere. Moreover, since the velocity is purely radial $\mathrm{u}={\mathit{u}}_{\mathit{r}}\widehat{\mathrm{r}}$, we have simply u_z = μ_zu_r, with radial projection cosine μ_z = z∕r.

In the present 2D pseudo-planar formulation,variations can occur in both radius r and a lateral orthogonal direction y, taken to lie in the equatorial plane of symmetry. The ray directions z now have local projection cosines μ_r and μ_y relative to the r and y axes, with thus $u_z={\mu }_r{u}_r+{\mu }_y{u}_y=\mu u_r+sin\varphi \sqrt{1-{\mu }^2}{u}_y$, with θ ≡ arccos μ and ϕ the customary radiation angles in Sect. 2. Our computations include one purely radial ray, with μ_r = 1 and μ_y = 0, so that u_z = u_r(r, y); as noted in Sect. 2, we also formally account for four additional rays that all have ${\mu }_r=1∕\sqrt{3}$, with two pairs of rays with ${\mu }_y=\pm 1∕\sqrt{3}$, but each pair forming mirror projections above or below the r − y plane. In practice, because of the mirror symmetry about this plane, explicit computation is only needed for one pair; the other pair is simply accounted for by doubling the quadrature weights (see Fig. A.1). For notational convenience, let us denote this triad with an index k = −1, 0, 1, such that μ_r,0 = 1 and ${\mu }_{r,\pm 1}=1∕\sqrt{3}$, while μ_y,0 = 0 and ${\mu }_{y,\pm 1}=\pm 1∕\sqrt{3}$ (see Fig. A.1).

Fig. A.1 Ray trajectories in the prograde (k = +1; blue), radial (k = 0; black), and retrograde (k = −1; red) directions, crossing grid nodes (black dots) that neighbor a central node with spatial indices {i, j}. The upper panel shows the full 3D geometry of the radiation rays. Since conditions are assumed constant in x (and thus symmetric about x = 0), ray integrations computed along dashed and solid lines of the same color are identical, and so can be accounted for by simply doing one prograde (blue) and one retrograde (red) integration, and then giving these double weight in the angle quadrature. The lower panel shows thesefinal three distinct rays projected on the 2D r − y calculation plane.

For our uniform spatial grid with fixed spacings Δ r = Δy = Δ = 0.001R_*, we have coordinates r_i = R_* + iΔ and y_j = jΔ, for grid indices i =1 to n_r = 1000 and j =1 to n_y = 100. At each grid node {i, j}, the outward (+) increment in optical depth Δ t_+,ijk(x) along each of the directional triads k is computed from (A.1), assuming a piecewise linear variation of density ρ and velocities u_r and u_y to the next outer grid node, with indices {i + 1, j + k}. Summation from the lower boundary at the stellar surface then gives the associated outwardly integrated optical depths along each direction k to a node with coordinates {r, y}, $$t_{+,k}\left(x,r,y\right)=\sum _{i,j}\mathrm{\Delta }{t}_{+,\mathrm{ijk}}\left(x\right)+t_{+,k}\left(x,R_{*},y_{*}\right),$$(A.2)

where the summation is understood to be over all i below the index for r, and over the associated j variation for each particular ray k; the assumed periodic variation in y means that j indices are simply mapped onto the range 0 < j < 100 by taking mod(j, 100). This means that all rays considered here in the pseudo-planar model hit the stellar surface at the lower boundary of the grid. The surface boundary value allows one to account for a photospheric line profile and the effect of a cutoff at a maximum opacity κ_max in the line distribution, as given by Eq. (66) in OP96: $$t_{+,k}\left(x,R_{*},y_{*}\right)=\frac{{\kappa }_{\mathrm{\text{o}}}}{{\kappa }_{\mathrm{\text{max}}}}+\frac{{\kappa }_{\mathrm{\text{o}}}\varphi \left(x\right)}{{\kappa }_{\mathrm{\text{e}}}}.$$(A.3)

With the outward optical depths t_+,k(x, r, y) in hand, the computation of the resulting radial component of the line acceleration follows much the same approach as for the 1D-SSF formalism given in Sect. 5.3 of OP96, as further elaborated in Sect. 2.1 of SO15. The direct absorption component of gravitationally scaled line acceleration thus takes the form (SO15, Eq. (13)) $${\mathrm{\Gamma }}_{\mathrm{\text{dir}},r}\left(r,y\right)={\mathrm{\Gamma }}_{\mathrm{\text{thin}}}\sum _{x,k}{w}_x{w}_k\varphi \left(x-u_{z,k}\right)t_{+,k}^{-\alpha }\left(x,r,y\right),$$(A.4)

where u_z,k ≡ μ_r,ku_r(r, y) + μ_y,ku_y(r, y), and the optically thin normalization Γ _thin is given by SO15, Eq. (4). As in 1D models, the quadrature in frequency x is uniform with equal weights w_x = Γ(α)∕n_x, and a resolution of three points per thermal doppler width, Δ x = 1∕3; but the 1D angle-quadrature weight w_y is now replaced with the triad w_k, with normalized values w₀ = 0.211325 and w_±1 = 0.394338. This triad weights the oblique rays such that they cover the full μ-space from 0 to $1∕\sqrt{3}$, plus half the space from $1∕\sqrt{3}$ to a radial ray at 1; the radial ray thus has weight $w_0=\left(1-1∕\sqrt{3}\right)∕2=0.211235$.

Within this SSF formalism, the diffuse (scattering) component of the line force is again (see OP96) formed by using t₊ to build an associated inward optical depth t₋ , and using this to form a difference between the outward vs. inward escape probability, as given in Eq. (6) in SO15. To avoid the variability of a nodal topology at the wind base (see Sect. 3 of SO15), we assume a simple optically thin source function computed from a uniformly bright surface without limb darkening, as given in Eq. (8) in SO15.

Our computation of the lateral (y) component of the line force warrants some further elaboration. While the overall formulation is similar, this now depends the difference between the escape probabilities in prograde (k = +1) and retrograde (k = −1) directions, applied to both the direct and diffuse components (radial ray k = 0 plays no role). Defining the profile-averaged, outward (+) escape probabilities in the prograde/retrograde (±) directions as $$(\mathrm{A}.5)$$(A.5)

we can write the direct component of the lateral line acceleration (still scaled by the radial gravity) as $${\mathrm{\Gamma }}_{\mathrm{\text{dir}},y}\left(r,y\right)={\mathrm{\Gamma }}_{\mathrm{\text{thin}}}f\left(r\right)\left[b_{+,+}\left(r,y\right)-b_{+,-}\left(r,y\right)\right].$$(A.6)

To account for the additional radial dropoff associated with angular shrinking of the stellar core (e.g., Gayley & Owocki 2000), we include here a correction factor $$f\left(r\right)=w_1{\mu }_{y,+1}\frac{{R_{*}}^2}{r^2}.$$(A.7)

Defining inward (–) escape probabilities in a way analogous to (A.5), we can write the associated diffuse component of the lateral line acceleration as $${\mathrm{\Gamma }}_{\mathrm{\text{diff}},y}\left(r,y\right)={\mathrm{\Gamma }}_{\mathrm{\text{thin}}}s\left(r\right)\left[b_{-,+}+b_{+,-}-b_{+,+}-b_{-,-}\right],$$(A.8)

where the optically thin source function factor s(r) is given by Eq. (8) in SO15. For both the radial (r) and lateral (y) components, the associated total acceleration is given by the sum of the direct and diffuse contributions, Γ _tot = Γ_dir + Γ_diff.

In our numerical radiation-hydrodynamics simulations, we apply these total radial and lateral line accelerations in the associated radial and lateral momentum equations, $$(\mathrm{A}.9)$$(A.9) $$(\mathrm{A}.10)$$(A.10)

where P is the gas pressure, M_eff ≡ M(1 − Γ_e), and the last term in each equation corrects our pseudo-planar treatment for curvilinear coordinate effects (centrifugal and coriolis forces) in a spherical outflow. The density ρ is evolved according to the mass continuity equation $$\frac{\partial \rho }{\mathrm{\partial t}}+\frac{\partial \left(\rho v_r\right)}{\mathrm{\partial r}}+\frac{\partial \left(\rho v_y\right)}{\mathrm{\partial y}}=-\frac{2\rho v_r}{r},$$(A.11)

where the source term on the right-hand side corrects for the neglect of the spherical divergence within our pseudo-planar treatment of the flow divergence ∇⋅ (ρv). The 1∕r² decline of the radiative flux, which sets the scale of the line accelerations, is accounted for by scaling these accelerations with the inverse-square decline of the stellar gravity.

References

All Tables

All Figures

Fig. 1 Sketch illustrating the basic idea of the pseudo-planar, box-in-a-wind approach used in this paper. The upper left illustrates the general situation of non-alignment between oblique rays and the numerical grid points. The lower panel then shows how we create a pseudo-planar box in the wind by cutting out a small but representative fraction of the wind volume. For illustration purposes, it shows projections onto the equatorial plane of rays in the prograde (blue), retrograde (red), and radial (black) directions for two lateral periods of a simple case with just ny = 2 zones in lateral direction y. The right panel then illustrates how extension out of the equatorial plane involves a total of five rays: one radial plus two oblique pairs that extend up/down from the plane. See Appendix A for a detailed explanation and for further illustrations of the assumed ray geometry.

In the text

Fig. 2 Spatial and temporal variations of log density relative tothe initial, smooth CAK steady state at t = 0, with color ranging from densities a decade below the t = 0 value (blue) to a decade above (red). The vertical variation extends from the subsonic wind base at the stellar surface R* to a height of one R* above. For clarity, the lateral variation is displayed over twice the horizontal box length 0.1R* . The upper row shows time evolution over the initial 100 ks after the CAK initial condition in steps of 10 ks; the bottom row uses the same step size of 10 ks to show the evolution between 300 and 400 ks, long after the initial condition has developed into a statistically steady turbulent flow.

In the text

Fig. 3 As in Fig. 2, spatial and temporal variations of log density relative to the initial, smooth CAK steady state at t = 0 are shown, with color ranging from densities a decade below the t = 0 value (blue) to a decade above (red). Here the vertical variation only extends between 1.9R* and 2.0R* and the lateral variation is displayed over one horizontal box of 0.1R* ; there are thus 100 × 100 discrete mesh points in each of the displayed squares. From left to right are shown a 2 ks time evolution long after the initial condition, in steps of 0.5 ks.

In the text

Fig. 4 Spatial and temporal variations of radial velocity vrad , with color ranging from 0 (blue) to 2000 kms−1 (red). As in Fig. 2, the vertical variation extends from the subsonic wind base at the stellar surface R* to a height of one R* above, and the lateral variation is displayed over twice the horizontal box length 0.1R* . The frames from left to right show the time evolution of vrad over 400 ksafter the CAK initial condition, in steps of 50 ks.

In the text

	Fig. 5 Radial cuts through the 2D simulation box of density ρ [g cm−3 ] (left) and radial velocity vrad [cm s−1] (right). Thered curves are taken at a snapshot long after the simulation has developed into a statistically quite steady flow; the black curves compare this to average values.
In the text

Fig. 6 Selected statistical properties of the 2D simulation (see text). The upper left panel plots the clumping factor fcl ; the upper right panel shows the time-dependent mass-loss rate, Ṁ [M⊙ yr−1] vs. s, computed in two different ways for the red and black curves(see text); the lower left panel displays lateral (black) andradial (red) density correlation lengths, and the corresponding Gaussian fit (blue, dashed); the lower right panel then finally plots radial (left) and lateral (right) velocity dispersions.

In the text

Fig. 7 Spatial and temporal variations of log density, radial velocity, and lateral velocity for a modelwith stellar rotation at the surface vy = 300 km s−1 (see text), with color-coding as in Figs. 2 and 4. The vertical variation in this simulation extends only from 1.0 to 1.5 R* , but the lateral variation is displayed as before over twice the horizontal box length 0.1R* . From left to right are shown the time evolution over 350 ks after the CAK initial condition, in steps of 50 ks.

In the text

Fig. A.1 Ray trajectories in the prograde (k = +1; blue), radial (k = 0; black), and retrograde (k = −1; red) directions, crossing grid nodes (black dots) that neighbor a central node with spatial indices {i, j}. The upper panel shows the full 3D geometry of the radiation rays. Since conditions are assumed constant in x (and thus symmetric about x = 0), ray integrations computed along dashed and solid lines of the same color are identical, and so can be accounted for by simply doing one prograde (blue) and one retrograde (red) integration, and then giving these double weight in the angle quadrature. The lower panel shows thesefinal three distinct rays projected on the 2D r − y calculation plane.

In the text