Issue 
A&A
Volume 611, March 2018



Article Number  A17  
Number of page(s)  10  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/201731718  
Published online  19 March 2018 
2D wind clumping in hot, massive stars from hydrodynamical linedriven instability simulations using a pseudoplanar approach
^{1}
KU Leuven, Instituut voor Sterrenkunde,
Celestijnenlaan 200D,
3001
Leuven, Belgium
email: jon.sundqvist@kuleuven.be
^{2}
Department of Physics and Astronomy, Bartol Research Institute, University of Delaware,
Newark,
DE 19716, USA
^{3}
Universitätssternwarte München,
Scheinerstr. 1,
81679
München, Germany
Received:
4
August
2017
Accepted:
20
October
2017
Context. Clumping in the radiationdriven winds of hot, massive stars arises naturally due to the strong, intrinsic instability of linedriving (the linedeshadowing instability, hereafter LDI). But LDI wind models have so far mostly been limited to 1D, mainly because of the severe computational challenges regarding calculation of the multidimensional radiation force.
Aim. In this paper we simulate and examine the dynamics and multidimensional nature of wind structure resulting from the LDI.
Methods. We introduce a pseudoplanar, boxinawind method that allows us to efficiently compute the line force in the radial and lateral directions, and then use this approach to carry out 2D radiationhydrodynamical simulations of the timedependent wind.
Results. Our 2D simulations show that the LDI first manifests itself by mimicking the typical shell structure seen in 1D models, but that these shells quickly break up into complex 2D density and velocity structures, characterized by smallscale density “clumps” embedded in larger regions of fast and rarefied gas. Key results of the simulations are that density variations in the welldeveloped wind are statistically quite isotropic and that characteristic length scales are small; a typical clump size is ℓ_{cl}∕R_{*}~ 0.01 at 2R_{*}, thus also resulting in rather low typical clump masses m_{cl} ~ 10^{17} g. Overall, our results agree well with the theoretical expectation that the characteristic scale for LDI generated windstructure is on the order of the Sobolev length ℓ_{Sob}. We further confirm some earlier results that lateral “filling in” of radially compressed gas leads to somewhat lower clumping factors in 2D simulations than in comparable 1D models. We conclude by discussing an extension of our method toward rotating LDI wind models that exhibit an intriguing combination of large and smallscale structures extending down to the wind base.
Key words: radiation: dynamics / hydrodynamics / instabilities / stars: earlytype / stars: mass loss / stars: winds, outflows
© 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 linedriving was given in the seminal paper by Castor et al. (1975; hereafter CAK). Like many wind models to date, CAK used the socalled Sobolev approximation (Sobolev 1960) to compute the radiative acceleration. This assumes that hydrodynamic flow quantities^{1} 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 ), allowing then a local treatment of the line radiative transfer.
This Sobolev approach ignores the strong linedeshadowing instability (LDI) that occurs on scales near and below the Sobolev length (Owocki & Rybicki 1984); numerical radiationhydrodynamic modeling of the nonlinear evolution of the LDI shows that the timedependent 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 Xray emission and broad Xray lines observed by orbiting telescopes like Chandra and XMMNewton (Feldmeier et al. 1997; Berghoefer et al. 1997; Güdel & Nazé 2009; Cohen et al. 2010; MartínezNúñ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 abovementioned 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 escapeintegral 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 multidimensional 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 rayintegrations 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_{r}n_{ϕ} integrations on the order of n_{r}n_{x} operations are needed for every considered ray; this gives an overall scaling , implying for a typical case of n_{ray} ≈ 5, n_{ϕ} ≈ 100, and n_{x} ≈ n_{r} ≈ 1000 on the order of 10^{11−12} 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 hotstarwind outflow, and a total simulation time of, e.g., ~ 50 dynamical time scales t_{dyn} = R_{*}∕v_{∞}~ 10 ks, requires approximately ~10^{5} repeated evaluations of the radiative force. This simple example thus illustrates quite vividly the rather daunting task of constructing multidimensional 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 singleray 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 oneray 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 pseudoplanar modeling approach for a multidimensional wind subject to the LDI. In this boxinawind 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 longcharacteristic rays with a computational costscaling 3n_{x}n_{r}n_{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 vh1^{2} to evolve the conservation equations of mass and momentum for a 2D, isothermal linedriven 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 pseudoplanar, boxinawind approach allows us to resolve laterally the relevant clump length scales, and to implement nonradial rays for the radiative linedriving 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 Ostar in the Galaxy. The standard setup uses a spatial grid with 1000 discrete radial (r) meshpoints 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 and the resulting smallscale 2D clump structures in the supersonic wind (the focus of this paper). Each simulation evolves from a smooth, CAKlike initial condition, computed by relaxing to a steady state a 1D spherically symmetric timedependent 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.
Summary of stellar and wind parameters.
2.2 Radiative driving
The central challenge in these simulations is to compute the 2D radiation line force in a highly structured, timedependent wind with a nonmonotonic velocity. This requires nonlocal integrations of the linetransport 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 multidimensional pseudoplanar 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 2DSSFformulation; below follows a summary of the key features.
Our pseudoplanar 2DSSF 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 linestrength number distribution is given by an exponentially truncated power law. In this formalism, α is the standard CAK powerlaw index, which can be physically interpreted as the ratio of the line force due to optically thick lines to the total line force; is a linestrength 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 linestrength cutoff^{3}. For typical Ostar conditions at solar metallicity, (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 multidimensional 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 pseudoplanar method introduced here largely circumvents these issues of grid misalignment. Namely, while radial rayintegrations 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 become constant throughout the computational domain. To this end, we apply a set of five rays with (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 “flatland” 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 “longcharacteristic” ray integrations at each time step with a factor of ~n_{r} ( = 10^{3} for our standard setup here).
Another attractive feature of this pseudoplanar 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 radialray 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 hotstar winds indeed extends to very nearphotospheric 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 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 pseudoplanar, boxinawind approach used in this paper. The upper left illustrates the general situation of nonalignment between oblique rays and the numerical grid points. The lower panel then shows how we create a pseudoplanar 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 n_{y } = 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 highdensityclumps 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 highvelocity streams as in corresponding 1D simulations; however, the velocity now also exhibits extensive lateral variation, reflecting again the breakup of 1D shells into smallscale 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 highspeed 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 v_{rad } , 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 v_{rad } 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 v_{rad } [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 (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 nearphotospheric 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 timedependent massloss rate: (2)
The black line in this plot shows a simple lateral average of the massflux escaping the outermost radial gridpoint at a specifictime. However, since our simulation box only covers 0.1R_{*}, this average very likely overestimates the timedependent 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 As expected, this curve shows a drastically lower temporal variation of Ṁ , despite the highly timedependent flow. This is consistent with, e.g., decadelong observations of spectral lines in Ostars, which typically indicate that time variations in the massloss rate of such stars are small.
To estimate typical clump length scales, the lower left panel of Fig. 6 plots a density autocorrelation length (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 welldeveloped 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 , where the estimated clump density here simply reads off the output of the simulations (see Fig. 5). More generally, such a clump massestimate may be obtained using the Sobolev length and mass conservation (4)
which for the 2D simulation analyzed here indeed gives m_{cl} ≈ 10^{17} 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 (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^{−1} in the outer wind.
Fig. 6 Selected statistical properties of the 2D simulation (see text). The upper left panel plots the clumping factor f_{cl } ; the upper right panel shows the timedependent massloss 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 pseudoplanar, boxinawind approach suitable for carrying out radiationhydrodynamical 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 smallscale 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 smallscale 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 correlationlength analysis also shows that, statistically, density variations in the welldeveloped 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^{17} 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^{−1} at the surface, and an initial condition set by steadystate angular momentum conservation, v_{y } (r) = v_{rot}R_{*}∕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 smallscale 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 smallscale structures needed to explain in parallel various observed phenomena in hotstar winds, such as discrete absorption components (DACs; Kaper et al. 1999) and smallscale 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 v_{y } = 300 km s^{−1} (see text), with colorcoding 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 followup 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 highmass Xray binary system (elMellah 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 pseudoplanar boxinawind.
Uncited reference
Acknowledgments
This work was supported in part by SAO Chandra grant TM314001A awarded to the University of Delaware, and in part by the visiting professor scholarship ZKD133200D01 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 longcharacteristic methods. We finally thank the referee for the useful comments on the paper.
Appendix A 2D pseudoplanar line force
Our development here of a 2D vector form for the line acceleration follows a direct generalization of the 1DSSF 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 profileweighted line optical depth between two wind locations along a ray coordinate z, (A.1)
where κ_{o} is a spatially constant line opacity normalization^{4}, and u_{z} ≡ẑ ⋅ u is the local zprojection 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, , with x = (ν − ν_{o})∕Δν_{D} the observerframe frequency displacement from line center in thermal Doppler units Δν_{D} = ν_{o}v_{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 , and the signtaken to be positive (negative) in the forward (backward) hemisphere. Moreover, since the velocity is purely radial , we have simply u_{z } = μ_{z}u_{r}, with radial projection cosine μ_{z} = z∕r.
In the present 2D pseudoplanar 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 , 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 , with two pairs of rays with , 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 , while μ_{y,0} = 0 and (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}, (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 pseudoplanar 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: (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 1DSSF 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)) (A.4)
where u_{z,k} ≡ μ_{r,k}u_{r}(r, y) + μ_{y,k}u_{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 anglequadrature weight w_{y } is now replaced with the triad w_{k}, with normalized values w_{0} = 0.211325 and w_{±1} = 0.394338. This triad weights the oblique rays such that they cover the full μspace from 0 to , plus half the space from to a radial ray at 1; the radial ray thus has weight .
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 profileaveraged, outward (+) escape probabilities in the prograde/retrograde (±) directions as (A.5)
we can write the direct component of the lateral line acceleration (still scaled by the radial gravity) as (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 (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 (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 radiationhydrodynamics simulations, we apply these total radial and lateral line accelerations in the associated radial and lateral momentum equations, (A.9) (A.10)
where P is the gas pressure, M_{eff} ≡ M(1 − Γ_{e}), and the last term in each equation corrects our pseudoplanar treatment for curvilinear coordinate effects (centrifugal and coriolis forces) in a spherical outflow. The density ρ is evolved according to the mass continuity equation (A.11)
where the source term on the righthand side corrects for the neglect of the spherical divergence within our pseudoplanar treatment of the flow divergence ∇⋅ (ρv). The 1∕r^{2} decline of the radiative flux, which sets the scale of the line accelerations, is accounted for by scaling these accelerations with the inversesquare decline of the stellar gravity.
References
 Berghoefer, T. W., Schmitt, J. H. M. M., Danner, R., & Cassinelli, J. P. 1997, A&A, 322, 167 [NASA ADS] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Cohen, D. H., Leutenegger, M. A., Wollman, E. E., et al. 2010, MNRAS, 405, 2391 [NASA ADS] [Google Scholar]
 Cohen, D. H., Wollman, E. E., Leutenegger, M. A., et al. 2014, MNRAS, 439, 908 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Dessart, L., & Owocki, S. P. 2003, A&A, 406, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dessart, L., & Owocki, S. P. 2005a, A&A, 437, 657 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dessart, L., & Owocki, S. P. 2005b, A&A, 432, 281 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 elMellah, I., Sundqvist, J. O., & Keppens, R. 2018, MNRAS, 475, 3240 [NASA ADS] [CrossRef] [Google Scholar]
 Eversberg, T., Lepine, S., & Moffat, A. F. J. 1998, ApJ, 494, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Feldmeier, A., Puls, J., & Pauldrach, A. W. A. 1997, A&A, 322, 878 [NASA ADS] [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Gayley, K. G., & Owocki, S. P. 2000, ApJ, 537, 461 [NASA ADS] [CrossRef] [Google Scholar]
 Güdel, M., & Nazé, Y. 2009, A&ARv, 17, 309 [NASA ADS] [CrossRef] [Google Scholar]
 Kaper, L., Henrichs, H. F., Nichols, J. S., & Telting, J. H. 1999, A&A, 344, 231 [NASA ADS] [Google Scholar]
 Lépine, S., & Moffat, A. F. J. 2008, AJ, 136, 548 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1983, ApJ, 274, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1984, ApJ, 284, 351 [NASA ADS] [CrossRef] [Google Scholar]
 MartínezNúñez, S., Kretschmar, P., Bozzo, E., et al. 2017, Space Sci. Rev., 212, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P. 1991, in Stellar Atmospheres – Beyond Classical Models, eds. L. Crivellari, I. Hubeny, & D. G. Hummer, NATO ASIC Proc., 341, 235 [Google Scholar]
 Owocki, S. P. 1999, in Variable and Nonspherical Stellar Winds in Luminous Hot Stars, Lect. Notes Phys., 523, eds. B. Wolf, O. Stahl, & A. W. Fullerton (Berlin: Springer Verlag), IAU Colloq., 169, 294 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & Puls, J. 1996, ApJ, 462, 894 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & Rybicki, G. B. 1984, ApJ, 284, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & Rybicki, G. B. 1985, ApJ, 299, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & Sundqvist, J. O. 2018, MNRAS, 475, 814 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Poe, C. H., Owocki, S. P., & Castor, J. I. 1990, ApJ, 358, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Puls, J., Owocki, S. P., & Fullerton, A. W. 1993, A&A, 279, 457 [NASA ADS] [Google Scholar]
 Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Vink, J. S., & Najarro, F. 2008, A&ARv, 16, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rybicki, G. B., Owocki, S. P., & Castor, J. I. 1990, ApJ, 349, 274 [NASA ADS] [CrossRef] [Google Scholar]
 Sobolev, V. V. 1960, Moving Envelopes of Stars (Cambridge: Harvard University Press) [CrossRef] [Google Scholar]
 Sundqvist, J. O., & Owocki, S. P. 2013, MNRAS, 428, 1837 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., & Owocki, S. P. 2015, MNRAS, 453, 3428 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., Puls, J., & Feldmeier, A. 2010, A&A, 510, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
The vh1 hydrodynamics computercode package has been developed by J. Blondin and collaborators, and is available for download at: http://wonka.physics.ncsu.edu/pub/VH1/
We have recast the line force using the notation of Gayley & Owocki (2000) rather than the κ_{0} notation of OP96. has the advantage of being a dimensionless measure of line strength that is independent of the thermal speed. The relation between the two parameter formulations is given in Appendix A.
In the notation of Gayley (1995), the line normalization here is given by , where Γ(α) is the complete Gamma function, and the numerical values used here are given in Table 1.
All Tables
All Figures
Fig. 1 Sketch illustrating the basic idea of the pseudoplanar, boxinawind approach used in this paper. The upper left illustrates the general situation of nonalignment between oblique rays and the numerical grid points. The lower panel then shows how we create a pseudoplanar 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 n_{y } = 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 v_{rad } , 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 v_{rad } 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 v_{rad } [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 f_{cl } ; the upper right panel shows the timedependent massloss 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 v_{y } = 300 km s^{−1} (see text), with colorcoding 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.