Issue 
A&A
Volume 631, November 2019



Article Number  A172  
Number of page(s)  9  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/201936331  
Published online  19 November 2019 
Theoretical wind clumping predictions of OB supergiants from linedriven instability simulations across the bistability jump
Institute of Astronomy, KU Leuven,
Celestijnenlaan 200D box 2401,
3001
Leuven,
Belgium
email: florian.driessen@kuleuven.be
Received:
17
July
2019
Accepted:
23
October
2019
Context. The behaviour of mass loss across the socalled bistability jump, where iron recombines from Fe IV to Fe III, is a key uncertainty in models of massive stars. Specifically, while an increase in mass loss is theoretically predicted, this has not yet been observationally confirmed. However, radiationdriven winds of hot massive stars are known to exhibit clumpy structures triggered by the linedeshadowing instability (LDI). This wind clumping severely affects empirical massloss rates inferred from ρ^{2}dependent spectral diagnostics. Thus, if clumping properties differ significantly for O and B supergiants across the bistability jump, this may help alleviate current discrepancies between theory and observations.
Aims. We investigated with analyt ical and numerical tools how the onset of clumpy structures behave in the winds of O supergiants (OSG) and B supergiants (BSG) across the bistability jump.
Methods. We derived a scaling relation for the linear growth rate of the LDI for a single optically thick line and applied it in the OSG and BSG regime. We ran 1D timedependent linedriven instability simulations to study the nonlinear evolution of the LDI in clumpy OSG and BSG winds.
Results. Linear perturbation analysis for a single line shows that the LDI linear growth rate Ω scales strongly with stellar effective temperature and terminal wind speed: Ω∝v_{∞}^{2}T_{eff}^{4}. This implies significantly lower growth rates for (the cooler and slower) BSG winds than for OSG winds. This is confirmed by the nonlinear simulations, which show significant differences in OSG and BSG wind structure formation, with the latter characterized by significantly weaker clumping factors and lower velocity dispersions. This suggests that lower correction factors due to clumping should be employed when deriving empirical massloss rates for BSGs on the cool side of the bistability jump. Moreover, the nonlinear simulations provide a theoretical background towards explaining the general lack of observed intrinsic Xray emission in single Bstar winds.
Key words: radiation: dynamics / hydrodynamics / instabilities / stars: earlytype / stars: massloss / stars: winds, outflows
© ESO 2019
1 Introduction
Determining and understanding massloss rates of massive stars is of great importance because the lives of these stars are heavily dependent on the amount of mass they lose throughout their lifetime (Chiosi & Maeder 1986). The mass expelled from the surface of hot massive OB stars is known to be driven by the scattered radiation off spectral lines (Castor et al. 1975, hereafter CAK) via a linedriven wind. Subsequent refinements to the CAK theory (Friend & Abbott 1986; Pauldrach et al. 1986) provided a basic quantitative understanding of the global properties of these linedriven winds, such as the metallicity dependence of massloss rate (Vink et al. 2001; Mokiem et al. 2007) and the relation between the stellar luminosity and wind momentum (Kudritzki et al. 1995).
Currently massloss rates of OB stars are either calculated from theoretical linedriven wind models (more below) or inferred from observational spectral diagnostics. A key problem regarding theory and observation is the socalled bistability jump (Pauldrach & Puls 1990) in the effective temperature range of 22.5–27.5 kK. At these temperatures iron recombines to a lower ionization state (Fe IV → Fe III) and has then more driving lines leading to a stronger radiation force. Across the bistability jump Lamers et al. (1995) observed a drop in terminal wind speed. However, linedriven wind models (Vink et al. 1999) predict an additional increase in massloss rate by a factor of several when transiting from the hot to the cool side of the bistability jump. Whereas a decrease in terminal wind speed is observationally established (Benaglia et al. 2007; Markova & Puls 2008; Haucke et al. 2018), the predicted increase in massloss rate has not been confirmed. This lack of massloss increase over the bistability jump might also have rather severe consequences for massive star evolution studies (Vink et al. 2010). Recently, Keszthelyi et al. (2017) concluded that without a bistability massloss increase, massive stars do not spin down fast enough to explain the rather large observed population of slowly rotating B supergiants.
However, there might be a bias regarding empirical massloss rates derived from observations across the bistability jump. There is ample evidence that linedriven winds are both spatially and temporally structured and variable (i.e. they are inhomogeneous clumpy winds (see e.g. Puls et al. 2008; Sundqvist et al. 2011, for reviews). The wind clumping can then seriously affect empirical massloss rates derived from observations. In particular, spectral diagnostics of ρ^{2}dependent atomic processes are affected by wind clumps (e.g. the optical recombination line Hα). To see this, consider how smallscale density inhomogeneities redistribute matter into overdense clumps and an effective void interclump medium. Using conservation of mass, the overall average density in the wind scales with the massloss rate as ⟨ρ⟩ ∝ Ṁ. The density of the clumps can then be written as ρ_{clump} = ⟨ρ⟩∕f_{vol} with f_{vol} accounting for the clump volume. Similarly, it follows that . To quantify the wind clumping forthese ρ^{2}dependent diagnostics, it is commonto define a clumping factor . In the case of recombinations between an electronproton pair, the opacity χ can be related to the massloss rate via . For the important Hα massloss diagnostic this means then that wind clumping leads to an overestimate of massloss rates with a dependency,if it is not properly accounted for in the diagnostic modelling (i.e. (Puls et al. 2006; Sundqvist & Puls 2018). This effect could explain the theoretically predicted increase in massloss rate at the bistability jump, if clumping factors change significantly across the temperature region where this occurs. A first clue that this might be true has been provided by Petrov et al. (2014), who found that Hα line profiles switch from strong to weak emission across the bistability jump. Inferring wind clumping properties is, however, a nontrivial task, and can be either done by a multitude of empirical spectral wind diagnostics (e.g. Puls et al. 2006)or by performing theoretical complex timedependent nonlinear hydrodynamic simulations (as presented here).
These hydrodynamic simulations are motivated by linear stability analysis that has shown that linedriving of hot star winds is in fact highly unstable due to a very strong, intrinsic linedeshadowing instability (LDI) occurring on small spatial scales (MacGregor et al. 1979; Carlberg 1980; Owocki & Rybicki 1984). Subsequent timedependent linedriven instability models (Owocki et al. 1988; Feldmeier 1995; Dessart & Owocki 2003; Sundqvist et al. 2018) have confirmed the picture of a highly inhomogeneous clumpy wind. Such smallscale structured clumpy winds also provide basic explanations of key observed signatures in the winds of massive stars, for example soft Xray emission (Berghoefer et al. 1997; Cohen et al. 2010; MartínezNúñez et al. 2017), extended regions of zero residual flux in some UV resonance lines (Lucy 1983; Sundqvist et al. 2010), or migrating spectral subpeaks in optical recombination lines (Eversberg et al. 1998; Dessart & Owocki 2005).
A key point, however, is that all these previous LDI simulations were carried out for stellar parameters typical of O supergiant winds. To date, no LDI model has been calculated for a B supergiant. The bistability jump happens where spectroscopic classification approximately separates the O and B star regimes, requiring new LDI models for the latter. These models provide a first step toward a better characterization of clumping properties and mass loss for B supergiants. In turn, this characterization of clumping and mass loss might then also provide better constraints regarding these objects’ rather poorly understood evolutionary status (Langer 2012).
Inspired by this, we perform timedependent linedriven instability simulations of both O and B supergiant winds. We carefully examine and contrast their clumping behaviour and discuss the implications. By means of analytical perturbation analysis we illustrate the difference in linear growth rate of the LDI between O and B supergiant winds. These analytical calculations are verified with nonlinear, numerical simulations of the LDI. We discuss whether wind clumping in OB supergiant winds undergoes a structural change across the bistability jump. Specifically, we look at wind clumping in the Hα lineformation region in the inner wind because of the importance of Hα as massloss rate diagnostic.
2 Scaling relation for the growth rate of unstable structures using perturbation analysis
Following Owocki & Rybicki (1984) we consider smallamplitude perturbations in velocity δv ∝ e^{i(kx−ωt)} applied to the unperturbed 1D momentum equation. Here k denotes the wavenumber and ω the angular frequency that can be either real, imaginary, or complex. For a spherically symmetric flow locally comoving with the mean flow, and in the limit of vanishing gas pressure, the firstorder perturbed momentum equation becomes (1)
with δg the response of the mean lineforce to the perturbed velocity. By considering an optically thick line in the mean flow and perturbations in wavelength and optical depth, Owocki & Rybicki (1984) derived a general bridging law for the perturbed lineforce, (2)
with Λ the bridging length, which is of the order of a radial Sobolev length L ≡ v_{th}∕(dv∕dr) ~ Λ, while Ω ≈ g_{Sob}∕v_{th} the growth rate of the instability, g_{Sob} the Sobolev lineforce, and v_{th} the thermal speed in the wind. We note that the perturbed quantities δg and δv are taken with respect to the background unperturbed steadystate mean flow described by the Sobolev theory for a single optically thick line. Mathematically this gives rise to the Sobolev dependency of the bridging length and growth rate of the instability, even though the instability itself is of nonSobolev nature.
It follows that longwavelength perturbations (kΛ ≪ 1) lead to and result in stable radiativeacoustic waves (Abbott 1980), while in the limit of shortwavelength perturbations (kΛ ≫ 1) we retrieve , which means an instability occurs (MacGregor et al. 1979; Carlberg 1980). This instability, the linedeshadowing instability, is very strong and is inherent to a linedriven wind.
The above bridging law is valid for pure photon absorption (i.e. perturbations in the direct lineforce component). Inclusion of linescattering, captured by the perturbed diffuse lineforce, results in a reduction of the growth rate near the stellar base (Owocki & Rybicki 1985). Since linescattering is part of our hydrodynamical models (detailed in Sect. 3), following Owocki & Puls (1996, hereafter OP96) a scatteringmodified bridging law is (3)
with x_{⋆} ≈ 1 the frequency at the blueabsorption edge (see OP96 for definition), and
Taking nowthe limit of shortwavelength unstable perturbations (kL ≫ 1), this becomes (4)
A simple scaling relation for the above perturbed lineforce then follows from applying the Sobolev theory in the optically thick limit for a single line (5)
Recalling that Ṁ ∝ r^{2}ρv, and considering the maximum growth rate in the limit far away from the star (s → 0.5; i.e. the supersonic regime in which the Sobolev approximation holds (Sobolev 1960): v ≫ v_{th}) (6)
where we have used in the last expression. A basic dependence of v_{∞}∕v_{th} appears in this relation that has also been found in the analysis of Owocki & Rybicki (1984), who considered a scaling of the growth rate compared to the typical wind expansion rate by approximating v(d v∕dr) ~ g_{CAK}, with g_{CAK} the total mean CAK lineforce.
If we plot this maximum Ω dependency for our models it can be clearly seen that the relative growth rate of the LDI differs significantly between O and B supergiants (Fig. 1). Between the hottest and coolest model there is a decrease in growth rate of about 90%. From our smallamplitude ansatz the amplitude of the instability grows with exp Ωt with Ω^{−1} a characteristic timescale (growth time). It can thus be understood that Ω^{−1} denotes the characteristic time at which the wave becomes amplified by a factor e, providing a relevant timescale of the problem. From this analysis it follows then that the onset of structure formation in O supergiants should occur much faster than in B supergiants, and that the resulting structure should also be more vigorous in O supergiants. As shown further below, this is exactly what we observe in our nonlinear, numerical simulations.
Fig. 1 Comparing maximum growth rates for OB supergiants as a function of the effective temperature using the scaling relation from Eq. (6). All models are normalized with respect to the hottest model to show the relative difference. The parameters used are listed in Table 1. 
3 Numerical simulations of clumpy winds
3.1 Hydrodynamics
To perform simulations of the nonlinear evolution of the radiatively unstable stellar wind we evolve the equations of 1D spherically symmetric isothermal radiationhydrodynamics, assuming the wind is at the same temperature as the stellar effective temperature. Table 1 shows the adopted stellar and wind parameters over a range of representative effective temperatures for OB supergiants in the Galaxy (see Sect. 4.1 for a further discussion of parameter choices). The simulations are performed using the piecewise parabolic method (PPM; Colella & Woodward 1984) Lagrangian remap hydrodynamics code VH1^{1}.
The key challenge in these instability simulations is a proper treatment of the radiation lineforce in a timedependent highly structured stellar wind outflow. As the instability acts on small spatial length scales, expensive nonlocal lineforce integrations have to be carried out at each hydrodynamical time step. It is important that such nonlocal lineforce computations should still closely approximate steadystate Sobolev models in order to be consistent. To this end we rely on the smooth source function (SSF) formulation (Owocki 1991) that has been extensively described by OP96.
All our model simulations start from a smooth, relaxed CAK wind model covering a nonuniform grid (1− 10 R_{⋆}) of 2000 cells. We apply a grid stretching factor s between subsequent radial cells i and i+ 1 such that s ≡ Δr_{i+1}∕Δr_{i} = 1.002. This allows the grid to start from a subsonic base that resolves the effective photospheric scale height. Inner boundary conditions are chosen such that the density ρ_{0} is fixed andwe set velocity by constant slope extrapolation. The lower boundary density ρ_{0} is about 5–7 times the density at the sonic point (Sundqvist & Owocki 2013). At the (supersonic) outer boundary all hydrodynamical variables are set by simple extrapolation assuming constant gradients.
3.2 Radiation linedriving
For our wind modelling we assume driving by an ensemble of lines described by a powerlaw line distribution, as originally proposed by CAK, in line strength q according to the formalism of Gayley (1995). Additionally, following OP96, we assume an exponential truncation of this line distribution limited to a maximum strength Q_{max} (7)
with Γ(α) the Gamma function. We further describe the linedriven wind using α, the CAK powerlaw index denoting the relative contribution of the lineforce from optically thick lines to the total lineforce, and as a linestrength normalization factor describing the ratio of the total lineforce to the electron scattering force if all lines were optically thin (Gayley 1995)^{2}. For the continuum opacities, we take, as is customary, an electron scattering contribution described by an opacity κ_{e} = 0.34 cm^{2} g^{−1}. Electron scattering effectively lowers gravity and is commonly captured by Eddington’s gamma Γ_{e} ≡ κ_{e}L_{⋆}∕(4πGM_{⋆}c). The effective escape speed at the stellar surface then becomes . This setup allows the computation of the average massloss rate ⟨Ṁ ⟩ and average terminal wind speed ⟨v_{∞}⟩ in our simulations.
We include a photospheric linear Eddington limb darkening description as motivated both by perturbation analysis and numerical simulations (Owocki & Rybicki 1985; Sundqvist & Owocki 2013). As studied extensively in the latter paper, this also leads to a net instability at the stellar surface and so to structure formation somewhat closer to the photosphere than in simulations assuming a uniformly bright stellar disk.
We refer the reader to Appendix A for a full description of how the perturbed lineforce is numerically computed.
Input stellar and wind parameters of our models.
4 Numerical results
4.1 Theoretical massloss rate predictions
It is important to point out that the bistability jump as predicted from linedriven wind models (Vink et al. 1999) implies an increase in massloss rate most directly when considering models with fixed mass and luminosity for which effective temperature is varied. Given that the internal parameter space of B supergiants is very scattered in fundamental stellar parameters like L_{⋆}, M_{⋆}, and R_{⋆} (Benaglia et al. 2007; Markova & Puls 2008; Haucke et al. 2018), we do not consider such fixed luminosity–mass models. Instead, we consider models where we select an effective temperature and choose a stellar radius that is typical in the observed range of OB supergiants. From this radius and effective temperature a stellar luminosity is determined. Similarly, we pick a reasonable stellar mass for each model based on observations of OB supergiants.
We should also point out that our nonlinear simulations do not fundamentally predict global massloss rates. We chose our model wind lineforce parameters (e.g. Puls et al. 2000) such that our relaxed CAK simulation massloss rates are approximately the same as those predicted from the Vink et al. (2001) recipe. Evolving such relaxed CAK models including the nonlinear instability then provides average massloss rates ⟨Ṁ ⟩ that are comparable to those of Vink et al. (2001). Similarly, we extract average terminal wind velocities ⟨v_{∞} ⟩ from our simulation and compute the ratio to the effective escape speed. The retrieved simulation quantities are shown in Table 2 for comparison.
Inspection of Table 1 shows that we have not computed any models that are part of the predicted bistability jump region of Vink et al. (1999). The reason for this is as outlined above: our steadystate CAK models are required to converge to massloss rates close to those predicted by Vink et al. (2001) using a realistic set of wind parameters. However, the Vink et al. (2001) massloss rate recipe breaks down in the bistability region such that there is no clear prediction for what should actually happen here. Therefore, it is hard to assess and justify the correctness of any initial, steadystate Sobolev model, and thus any linedriven instability simulation and wind clumping prediction within the bistability.
To that end our nonlinear simulations do not provide a direct test of the predicted increase in massloss rates across the bistability jump as found by Vink et al. (1999). However, our simulations can provide evidence of whether wind clumping properties change between O and B supergiants, which affects massloss rates derived from observational diagnostics. The observational lack of evidence for increased massloss rates might then be explained in terms of applying incorrect clumping factors for OB supergiants across the bistability jump.
4.2 Structure formation in O and B supergiant winds
We illustrate first the basic overall dynamics of an unstable 1D linedriven wind. Figure 2 shows the basic structure resulting from a nonlinear LDI simulation, taken after a long enough simulation time which we set here by evolving over sufficient dynamical timescales τ_{dyn} = R_{⋆}∕v_{∞}. This 1D snapshot illustrates the typical formation of (strong) shocks that compress matter into dense, spatially narrow clumps (really shells in 1D) separated by a highspeed rarefied medium. The models shown start from a relaxed CAK model and then evolve to a point where the wind reaches a statistically steady flow.
In Fig. 2 we compare our hottest O supergiant model (OSG1) with that of the coolest B supergiant (BSG2). The figure directly shows the structural changes that occur in shock formation and resulting clump formation, as predicted by linear perturbation analysis in the previous section (see also Fig. 1). Not only is the strength of structure formation different, but also its onset (i.e. the dynamical timescale between O and B supergiants) is quite different. A direct computation shows that for the O supergiant τ_{dyn} ≈ 14 ks and for the B supergiant τ_{dyn} ≈ 65 ks, meaning that B supergiants take longer to relax into a statistically steadystate flow.
Figure 3 shows a 2D contour of wind density ρ and wind velocity v over time for all our models. As motivated in the next section we concentrate here on the inner wind region where the Hα line is formed (r ≤ 2 R_{⋆}). All simulations extend to 300 ks for a direct comparison, but this only captures the start of structure formation in the B supergiant wind models, due to the longer dynamical timescale. In density space (top panels) there is a significant contrast between the O and B supergiant winds. Specifically, the B supergiant winds display much lower density and velocity contrasts between highdensity clumps (shells) and the rarefied regions in between them, as can be directly seen in the lower panel of Fig. 2. An O supergiant has regions of very rarefied gas compared to a B supergiant. The gaseousmaterial is much more accelerated in an O supergiant wind than in a B supergiant wind (bottom panels).
In accordance with our derived scaling relation above, the instability sets in much faster for the O supergiants than for the B supergiants. The O supergiants start to form smallscale structures at around 10 ks while the B supergiants only exhibit unstable behaviour at around 50 ks. This result comes quite close to the prediction made from Eq. (6): computing Ω^{−1} shows that the actual development of nonlinear structure starts about several ks to several 10 ks for O and B supergiants, respectively.If we let these simulations run on dynamical timescales that are long enough, reoccurring periodic structures appear in density and velocity space. This is already clearly visible in the O supergiant winds, while periodic structures still have to develop in the B supergiant winds. Essentially, limbdarkening lowers scattering contributions in the transonic wind region such that the wind solution becomes of unstable nodal type instead of stable X type. The nodal solution topology branch has a degenerate set of shallowslope solutions that are reached in semiregular intervals of the simulation (see Sundqvist & Owocki 2015, for more quantitative details).
Fig. 2 Velocity and density evolution of LDI in the wind of OB supergiants. The comparison is between our hottest and coolest model. The dashed line represents the average wind velocity and density profile. 
4.3 Wind clumping across the bistability jump
Having discussed the basic structures of the nonlinear evolution of the LDI we here present statistically computed properties for our OB supergiant models. The different dynamical properties between O and B supergiants also alter the wind clumping properties. To describe the clumping properties we use a clumping factor (8)
with brackets denoting time averages after a statistically steady flow has developed. Essentially, the O and B supergiant models are all averaged between 10–25τ_{dyn}, which translates to about 100–300 ks for the O supergiant models and about 700 ks–2 Ms for the B supergiant models. Following this approach we take into account the difference in dynamical timescale and make all the simulations insensitive to initial conditions.
The left panel of Fig. 4 shows the relative change of the radially varying wind clumping factor f_{cl} between ourmodels. The middle and right panels show average wind clumping factors over the full wind and inner wind (1− 2 R_{⋆}), respectively.We compute these averages for each of our models and show them as a function of effective temperature to illustrate the behaviour across the bistability jump. The average over the full wind already shows a significant difference between wind clumping properties on the hot and cool side of the bistability jump. As the inner wind is the region where the Hα line is being formed, it is important to consider a wind clumping average here. The Hα line is a prime spectral diagnostic used to derive empirical massloss rates across the bistability jump (e.g. Markova & Puls 2008). However, this Hα fitting also depends sensitively on the wind clumping properties, due to the ρ^{2} dependence (see Introduction). Effectively, when we derive Ṁ from Hα fitting, we are really deriving the product . This immediately shows that a degeneracy occurs when fitting line profiles because different combinations of Ṁ and f_{cl} can provide the same invariant .
On the hot side of the bistability jump, our simulations show f_{cl} ~ 17 in the Hα forming region. On the other hand, the cool side of the bistability jump shows much lower factors f_{cl} ~ 2. We note that our predicted wind clumping factors on the hot side of the bistability jump agree reasonably well with typicalvalues inferred from spectral diagnostical studies (Najarro et al. 2011; Sundqvist & Owocki 2013). This makes our study quite robust in a relative sense, even though there are still many uncertainties related to a quantitative treatment of the LDI (e.g. Sundqvist et al. 2018). On the cool side we find very low f_{cl} ≈ 2, suggesting that it is important to account correctly for the different wind clumping properties across the bistability jump. From a differential viewpoint, a key result here is thus that our models suggest that the winds of O and B supergiants have different structures in terms of clumping.
The lack of structure formation and the corresponding shocks is also important for intrinsic Xray emission of B supergiant winds. A key observational result is that the Xray flux drops in Bstar winds as compared to their Ostar counterparts (Cassinelli et al. 1994; Berghoefer et al. 1997; Güdel & Nazé 2009). This is in agreement with our simulations, although here we are not computing the shockheated regions directly. To facilitate a discussion around shocks, we define a velocity dispersion (9)
which is a measure of the standard deviation of wind velocity and describes the typical velocity around the mean wind flow speed (brackets denote again a time average). In the absence of direct Xray modelling, the velocity dispersion is a reasonable proxy for the overall intrinsic Xray emission expected from LDI shocks.
Figure 5 shows the velocity dispersion computed for our OB supergiant models. Clearly there is a large difference between the velocity dispersion in the models. For the O supergiant wind a steep increase in velocity dispersion results from the initial strong amplification of smallscale velocity structures near the stellar base by the LDI. Eventually these smallscale velocity structures situate themselves in the outer wind where almost no further amplification occurs. It is therefore not surprising to see that the B supergiant wind does not have a steep initial increase nor does it achieve large values in velocity dispersion.
Typical shock jump velocities in O supergiant winds are of the order of v ~ 300 km s^{−1} giving rise to shock temperatures of the order of T ~ 1 MK. A simple visual inspection of Fig. 5 shows our O supergiant models have v_{disp} ~ 300 km s^{−1}, which is enough toproduce intrinsic Xray emission, whereas the velocity dispersion in a B supergiant model is several factors lower (v_{disp} ~ 70 km s^{−1}), and so will not result in much (or any) shockheated gas hot enough to emit in the Xray band. Physically speaking, this is a combined consequence of the weaker clumping properties and lower wind speeds of B supergiant linedriven winds (see also previous figures).
Finally, we also note that the lower levels of velocity dispersion in B supergiant winds than in O supergiant winds are generally consistent with quantitative spectroscopic modelling of saturated UV P Cygni lines. Specifically, in this kind of diagnostic modelling the velocity dispersion is conventionally modelled using “wind microturbulence”; in order to match observationsthis microturbulence is typically assumed to scale with the local wind speed and reaches maxima of the order of 0.1v_{∞} (e.g. Puls et al. 1993). Since B supergiant winds have much lower v_{∞} than O supergiant winds, this agrees well with our findings here.
Fig. 3 Top row: logarithmic wind density structure. From left to right: OSG1, OSG2, BSG1, BSG2. Bottom row: wind velocity structure of the corresponding models. The rarefied medium in the density plot is characterized by high velocities, while the dense clumps move at much lower velocities. All models extend to 300 ks, but the B supergiant models actually run for a much longer time (see text). Also illustrated is the difference in dynamical timescale that governs the onset of structure formation. 
Fig. 4 Left panel: relative change in wind clumping between our models. The dashed lines indicate the inner wind where Hα line formation occurs and the typical value f_{cl} = 10 often used in diagnostic modelling. Middle panel: averaged clumping factor over the whole wind. Right panel: averaged clumping factor in the Hα region. The grey area represents the bistability region as predicted by Vink et al. (1999). 
Fig. 5 Velocity dispersion inside the wind computed from our simulations of OB supergiant models. 
5 Conclusions and future work
The central idea of this work is that theoretical predictions of clumping in O and B supergiant winds differ significantly from each other. With the help of analytical perturbation analysis we computed a scaling relation for the linear growth rate of the LDI for a single optically thick line. This calculation shows that there is a significant difference in instability growth rates between O and B supergiants. Likewise, we find that the onset of structure formation occurs at different times. We confirm these analytical findings using timedependent linedriven instability simulations that describe the nonlinear evolution of the LDI.
Extending these linedriven instability simulations to B supergiants has allowed us to study the behaviour of wind clumping across the predicted bistability jump in massive star winds (Pauldrach & Puls 1990; Lamers et al. 1995; Vink et al. 1999). A key result of our study is that simulations show a significant difference in wind clumping properties between the hot and cool sides of the bistability jump. Therefore, this suggests that accounting correctly for wind clumping is important in deriving empirical massloss rates across the bistability jump, if they are derived by means of the ρ^{2} sensitive Hα line. Specifically, in these diagnostic Hα studies of mass loss across the bistability jump, clumping effects are typically neglected entirely or treated as a constant across the bistability jump (e.g. Markova & Puls 2008). In turn, this may then also help alleviate current problems with massive star evolution models that do not assume a massloss increase over the bistability jump (Vink et al. 2010; Keszthelyi et al. 2017).
Building on the results of this paper, it would next be interesting to study in more detail linedriven instabilities across the internal parameter space of B supergiants. The stellar parameter space of B supergiants is very scattered (Benaglia et al. 2007; Markova & Puls 2008; Haucke et al. 2018) making it natural to ask to what extent such different sets of stellar parameters also affect clumping properties in the wind. Along these lines it would then also be interesting to consider models at fixed luminosity and mass like those of Vink et al. (1999). As directly evident from the linear scaling relation (see Eq. (6)), we indeed expect that such fixed luminositymass models would also exhibit significantly weaker wind clumping across the bistability jump.
Since the modelled velocity dispersion is much lower in BSG than OSG winds, our simulations provide a theoretical rationale to observations showing that Xray fluxes drop significantly in single Bstar winds compared to their Ostar counterparts (Cassinelli et al. 1994; Berghoefer et al. 1997; Güdel & Nazé 2009). This then also naturally explains the diagnostic finding that the wind microturbulence needed to fit observations of strong UV wind lines typically scales with the terminal speed. Moreover, the lack of strong shocks and structure formation might also affect wind accretion in highmass Xray binaries (HMXBs). These systems often consist of a B supergiant and an orbiting object that accretes material, for example the prototype system Vela X1 where the donor star is a early B supergiant with a slow and dense wind characteristic of a simulation below the bistability jump (Sander et al. 2018). However, in the accretion models by El Mellah et al. (2018) targeting Vela X1, a multiD LDI simulation of an O supergiant was used to simulate the clumpy wind accretion and the corresponding time variation in Xray luminosity. Our preliminary calculations suggest that also in such multiD models the structure formation in B supergiants is still much weaker than in O supergiants. Consequently, the amount of LDI generated accretion rate variations may be somewhat overpredicted in the models of El Mellah et al. (2018). In a followup work we plan to study in detail the effect of multiD LDI models of B supergiants on the accretion rate and the corresponding variations of Xray luminosities in such HMXBs.
Finally, we note that this work is a first theoretical study of the behaviour of wind clumping on both sides of the bistability jump. Our primary aim was to point out the difference in wind clumping properties in the O and B supergiant regimes. We show how the wind clumping properties change across the bistability jump and discuss whether this may provide an explanation of the lack of observational evidence for increasing massloss rates across the bistability jump. A quantitative wind clumping description depends on the treatment of limbdarkening, photospheric base perturbations, the source function, and the amount of lateral fragmentation occurring in much more complex multidimensional wind clumping models (Dessart & Owocki 2003; Sundqvist et al. 2018).
Moreover, our simulations only make predictions of clumping factors, but not the overall stellar massloss rate. Extending this study to also treat stars within the bistability jump region must be left to future work once more reliable predictions for the smooth steadystate winds in this domain are in place. However, we calibrated our 1D models so that they approximately agree with observed clumping factors on the hot side of the bistability jump. This means that, even thoughthere are still some significant uncertainties regarding the input physics used to model the LDI, the simulations presented here should be quite robust in a relative sense. Thus this paper provides a good first study of the differential effect of theoretical predictions of wind clumping of O and B supergiants across the bistability jump.
Acknowledgements
Many thanks to the anonymous referee for raising helpful comments that improved the manuscript. The authors wish to thank Stan Owocki for insightful discussions. We also thank the KU Leuven EQUATIONgroup for their valuable input as well as weekly supply of cake. F.A.D. and J.O.S. acknowledge support by the Belgian Research Foundation Flanders (FWO) Odysseus program under grant number G0H9218N. N.D.K. acknowledges support from the KU Leuven C1 grant MAESTRO C16/17/007.
Appendix A Computing the perturbed lineforce
As discussed in the Introduction, the Sobolev method is not suitable for computing the smallscale structures arising from the linedeshadowing instability as the instability acts on length scales smaller than the Sobolev length. This means that the lineforce no longer depends only on local flow quantities; instead, there is a nonlocal coupling between the radiation field and outflowing gas. Such full nonlocal calculations are computationally daunting, and approximate expressions for the lineforce have to be found. One general approximation is based on the escape integral probability method (OP96).
We define the escape probability for a single line, following the Gayley (1995) parametrization, as (A.1)
where ϕ is a profile function (taken to be Gaussian) at the observer’s frame frequency x ≡ (ν∕ν_{0} − 1)c∕v_{th} defined from linecentre ν_{0} in terms of the frequency broadening from ion thermal motions, μ is a local direction cosine at radius r, and qt is a frequency dependent optical depth for a line at strength q: (A.2)
Here we perform the integrations between two positions z_{1} < z_{2} along a certain ray p (impact parameter) from the origin where z ≡ μ∕r and r^{2} ≡ p^{2} + z^{2}. We note that κ in the above equation is not the usual opacity (expressed in units of [cm^{2} g^{−1}]), but rather is a frequency integrated lineopacity (units [cm^{2} g^{−1} Hz^{−1}]). Specifically, κ ≡ χ∕ρ where χ is a frequencyintegrated linestrength (e.g. as defined in Puls et al. 1993, see their Eq. (2)).
In the Sobolev approximation the velocity variation dominates the integral (v ≫ v_{th}) such that κρ becomes approximately constant over a Sobolev length L ≡ v_{th}∕(dv∕dr) (i.e. the integral is analytic). For physics occurring on spatial scales smaller than L (linedriven instabilities), this integral is inherently nonlocal and requires numerical evaluation.
The stellar wind off a hot, luminous star is not driven by a single line, but by a large collection of lines. For an ensemble of lines modelled by a linedistribution (see Sect. 3.2, also definition in OP96) we define an escape probability (A.3)
Starting from the definition of the radiative line acceleration it is possible to recast it in terms of the above ensemble escape probability (detailed derivation in OP96). The lineforce expressions for the direct and diffuse radiation become (A.4)
with the lineacceleration in the optically thin limit and g_{e} = κ_{e}L_{⋆}∕(4πr^{2}c) the acceleration due to continuum electron scattering. Angle brackets in the lineforce expressions denote average angle integrations. The total force due to radiation g_{rad} is the sum of the above continuum force and lineforces (i.e. g_{rad} ≡ g_{e} + g_{dir} + g_{diff}).
The function D(μ, r) contains any reductions of stellar intensity in case the star is not assumed to be a uniformly bright disk: D(μ, r) = 1 (μ_{⋆} < μ < 1) and D(μ, r) = 0 (μ < μ_{⋆}), where gives the angular size of the star. To account for the fact that the star is not a uniformly bright disk, we include such a D(μ, r) that describes a photospheric linear Eddington limbdarkening law (Hubeny & Mihalas 2014, Eq. (17.17)) (A.6)
which has the additional effect of reducing the stabilizing effect of diffuse radiation near the stellar base.
Additionally, in order to compute the diffuse lineforce in Eq. (A.5), we require an expression for the source function S(r). The optically thin limit turns out to be a good approximation and the source function can be closely approximated (OP96) by (A.7)
Plugging Eq. (A.6) into the above then gives (Sundqvist & Owocki 2013) (A.8)
For convenience it is useful to cast the angle integration of the lineforces in terms of a ray parameter with direction cosine . Applying this ray formulation and the Eddington limbdarkened source function into the expressions for the lineforces (i.e. Eqs. (A.4) and (A.5)) gives (A.9)
in which the diffuse lineforce follows a “twostream” approximation of inward (“−”, μ < 0) and outward (“+”, μ > 0) streaming photons. Additionally we correct with a factor as the ray integration should formally extend to (Owocki 1991). It has been shown that y = 0.5 is a good ray parameter choice and approximates the full stellar disk integration within a few percent for most of the wind (Sundqvist & Owocki 2013).
The outward escape probability b_{+}(μ_{y}, r) is computed according to Eq. (A.3) using an outward frequencydependent optical depth t_{+} ≡ qt(+x, y, r). The inward escape probability b_{−} (μ_{y}, r) requires a corresponding inward frequency dependent optical depth t_{−}, which can be computed from t_{+} by using (OP96) (A.11)
With the above equations, the lineforces are computed at each hydrodynamical time step and added together to get the total nonlocal lineforce acting on the outflowing stellar wind.
References
 Abbott, D. C. 1980, ApJ, 242, 1183 [NASA ADS] [CrossRef] [Google Scholar]
 Benaglia, P., Vink, J. S., Martí, J., et al. 2007, A&A, 467, 1265 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berghoefer, T.W., Schmitt, J. H. M. M., Danner, R., & Cassinelli, J. P. 1997, A&A, 322, 167 [NASA ADS] [Google Scholar]
 Carlberg, R. G. 1980, ApJ, 241, 1131 [NASA ADS] [CrossRef] [Google Scholar]
 Cassinelli, J. P., Cohen, D. H., Macfarlane, J. J., Sanders, W. T., & Welsh, B. Y. 1994, ApJ, 421, 705 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Chiosi, C., & Maeder, A. 1986, ARA&A, 24, 329 [NASA ADS] [CrossRef] [Google Scholar]
 Cohen, D. H., Leutenegger, M. A., Wollman, E. E., et al. 2010, MNRAS, 405, 2391 [NASA ADS] [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. 2005, A&A, 432, 281 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 El Mellah, I., Sundqvist, J. O., & Keppens, R. 2018, MNRAS, 475, 3240 [NASA ADS] [CrossRef] [Google Scholar]
 Eversberg, T., Lépine, S., & Moffat, A. F. J. 1998, ApJ, 494, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Feldmeier, A. 1995, A&A, 299, 523 [NASA ADS] [Google Scholar]
 Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701 [NASA ADS] [CrossRef] [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Güdel, M., & Nazé, Y. 2009, A&A Rev., 17, 309 [CrossRef] [Google Scholar]
 Haucke, M., Cidale, L. S., Venero, R. O. J., et al. 2018, A&A, 614, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton, NJ: Princeton University Press) [Google Scholar]
 Keszthelyi, Z., Puls, J., & Wade, G. A. 2017, A&A, 598, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kudritzki, R. P., Lennon, D. J., & Puls, J. 1995, in Science with the VLT, eds. J. R. Walsh, & I. J. Danziger (Berlin: Springer), 246 [Google Scholar]
 Lamers, H. J. G. L. M., Snow, T. P., & Lindholm, D. M. 1995, ApJ, 455, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Langer, N. 2012, ARA&A, 50, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1983, ApJ, 274, 372 [NASA ADS] [CrossRef] [Google Scholar]
 MacGregor, K. B., Hartmann, L., & Raymond, J. C. 1979, ApJ, 231, 514 [NASA ADS] [CrossRef] [Google Scholar]
 Markova, N., & Puls, J. 2008, A&A, 478, 823 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MartínezNúñez, S., Kretschmar, P., Bozzo, E., et al. 2017, Space Sci. Rev., 212, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Mokiem, M. R., de Koter, A., Vink, J. S., et al. 2007, A&A, 473, 603 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Najarro, F., Hanson, M. M., & Puls, J. 2011, A&A, 535, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owocki, S. P. 1991, in NATO Advanced Science Institutes (ASI) Series C, eds. L. Crivellari, I. Hubeny, & D. G. Hummer (Berlin: Springer), 341, 235 [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., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Pauldrach, A. W. A., & Puls, J. 1990, A&A, 237, 409 [NASA ADS] [Google Scholar]
 Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] [Google Scholar]
 Petrov, B., Vink, J. S., & Gräfener, G. 2014, A&A, 565, A62 [NASA ADS] [CrossRef] [EDP Sciences] [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., Markova, N., Scuderi, S., et al. 2006, A&A, 454, 625 [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]
 Sander, A. A. C., Fürst, F., Kretschmar, P., et al. 2018, A&A, 610, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sobolev, V. V. 1960, Moving Envelopes of Stars (Cambridge, US: Harvard University Press) [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. 2018, A&A, 619, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Puls, J., & Feldmeier, A. 2010, A&A, 510, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Puls, J., Feldmeier, A., & Owocki, S. P. 2011, A&A, 528, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Owocki, S. P., & Puls, J. 2018, A&A, 611, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 1999, A&A, 350, 181 [NASA ADS] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vink, J. S., Brott, I., Gräfener, G., et al. 2010, A&A, 512, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Here we deviate from the OP96 (κ, κ_{0}, κ_{max}) parametrization of the line distribution, i.e. our Eq. (7). Gayley’s parametrization (q, , Q_{max}) has the advantage of being a dimensionless measure of line strength that is independent of a fiducial thermal speed dependence. The relation between our parametrization and OP96’s follows from .
All Tables
All Figures
Fig. 1 Comparing maximum growth rates for OB supergiants as a function of the effective temperature using the scaling relation from Eq. (6). All models are normalized with respect to the hottest model to show the relative difference. The parameters used are listed in Table 1. 

In the text 
Fig. 2 Velocity and density evolution of LDI in the wind of OB supergiants. The comparison is between our hottest and coolest model. The dashed line represents the average wind velocity and density profile. 

In the text 
Fig. 3 Top row: logarithmic wind density structure. From left to right: OSG1, OSG2, BSG1, BSG2. Bottom row: wind velocity structure of the corresponding models. The rarefied medium in the density plot is characterized by high velocities, while the dense clumps move at much lower velocities. All models extend to 300 ks, but the B supergiant models actually run for a much longer time (see text). Also illustrated is the difference in dynamical timescale that governs the onset of structure formation. 

In the text 
Fig. 4 Left panel: relative change in wind clumping between our models. The dashed lines indicate the inner wind where Hα line formation occurs and the typical value f_{cl} = 10 often used in diagnostic modelling. Middle panel: averaged clumping factor over the whole wind. Right panel: averaged clumping factor in the Hα region. The grey area represents the bistability region as predicted by Vink et al. (1999). 

In the text 
Fig. 5 Velocity dispersion inside the wind computed from our simulations of OB supergiant models. 

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.