© ESO 2019

1 Introduction

Determining and understanding mass-loss 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 line-driven 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 line-driven winds, such as the metallicity dependence of mass-loss rate (Vink et al. 2001; Mokiem et al. 2007) and the relation between the stellar luminosity and wind momentum (Kudritzki et al. 1995).

Currently mass-loss rates of OB stars are either calculated from theoretical line-driven wind models (more below) or inferred from observational spectral diagnostics. A key problem regarding theory and observation is the so-called bi-stability 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 bi-stability jump Lamers et al. (1995) observed a drop in terminal wind speed. However, line-driven wind models (Vink et al. 1999) predict an additional increase in mass-loss rate by a factor of several when transiting from the hot to the cool side of the bi-stability 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 mass-loss rate has not been confirmed. This lack of mass-loss increase over the bi-stability 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 bi-stability mass-loss 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 mass-loss rates derived from observations across the bi-stability jump. There is ample evidence that line-driven 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 mass-loss rates derived from observations. In particular, spectral diagnostics of ρ²-dependent atomic processes are affected by wind clumps (e.g. the optical recombination line Hα). To see this, consider how small-scale 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 mass-loss 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 ${\rho }_{\text{clump}}^2=\langle {\rho }^2\rangle /f_{\text{vol}}$. To quantify the wind clumping forthese ρ²-dependent diagnostics, it is commonto define a clumping factor $f_{\text{cl}}\equiv \langle {\rho }^2\rangle /{\langle \rho \rangle }^2=1/f_{\text{vol}}\ge 1$. In the case of recombinations between an electron-proton pair, the opacity χ can be related to the mass-loss rate via $\chi \propto \langle {\rho }^2\rangle ={\langle \rho \rangle }^2{f}_{\text{cl}}\propto {\dot{M}}^2{f}_{\text{cl}}$. For the important Hα mass-loss diagnostic this means then that wind clumping leads to an overestimate of mass-loss rates with a $\sqrt{f_{\text{cl}}}$ dependency,if it is not properly accounted for in the diagnostic modelling (i.e. ${\dot{M}}_{\text{unclumped}}={\dot{M}}_{\text{clumped}}\sqrt{f_{\text{cl}}}$ (Puls et al. 2006; Sundqvist & Puls 2018). This effect could explain the theoretically predicted increase in mass-loss rate at the bi-stability 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 bi-stability jump. Inferring wind clumping properties is, however, a non-trivial 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 time-dependent non-linear hydrodynamic simulations (as presented here).

These hydrodynamic simulations are motivated by linear stability analysis that has shown that line-driving of hot star winds is in fact highly unstable due to a very strong, intrinsic line-deshadowing instability (LDI) occurring on small spatial scales (MacGregor et al. 1979; Carlberg 1980; Owocki & Rybicki 1984). Subsequent time-dependent line-driven 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 small-scale structured clumpy winds also provide basic explanations of key observed signatures in the winds of massive stars, for example soft X-ray emission (Berghoefer et al. 1997; Cohen et al. 2010; Martínez-Núñ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 bi-stability 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 time-dependent line-driven 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 non-linear, numerical simulations of the LDI. We discuss whether wind clumping in OB supergiant winds undergoes a structural change across the bi-stability jump. Specifically, we look at wind clumping in the Hα line-formation region in the inner wind because of the importance of Hα as mass-loss rate diagnostic.

2 Scaling relation for the growth rate of unstable structures using perturbation analysis

Following Owocki & Rybicki (1984) we consider small-amplitude 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 co-moving with the mean flow, and in the limit of vanishing gas pressure, the first-order perturbed momentum equation becomes $$\omega =i\frac{\delta g}{\delta v},$$(1)

with δg the response of the mean line-force 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 line-force, $$\frac{\delta g}{\delta v}\approx \Omega \frac{ik\Lambda }{1+ik\Lambda },$$(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 line-force, 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 steady-state 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 non-Sobolev nature.

It follows that long-wavelength perturbations (kΛ ≪ 1) lead to $\mathcal{R}[\omega ]$ and result in stable radiative-acoustic waves (Abbott 1980), while in the limit of short-wavelength perturbations (kΛ ≫ 1) we retrieve $\mathcal{I}[\omega ]$, which means an instability occurs (MacGregor et al. 1979; Carlberg 1980). This instability, the line-deshadowing instability, is very strong and is inherent to a line-driven wind.

The above bridging law is valid for pure photon absorption (i.e. perturbations in the direct line-force component). Inclusion of line-scattering, captured by the perturbed diffuse line-force, results in a reduction of the growth rate near the stellar base (Owocki & Rybicki 1985). Since line-scattering is part of our hydrodynamical models (detailed in Sect. 3), following Owocki & Puls (1996, hereafter OP96) a scattering-modified bridging law is $$\frac{\delta g}{\delta v}\approx \Omega \frac{kL}{1+{(kL/2x_{\star })}^2}\left[(1-s)\frac{kL}{2x_{\star }}\right],$$(3)

with x_⋆ ≈ 1 the frequency at the blue-absorption edge (see OP96 for definition), and $$s=\frac{1}{1+{\mu }_{\star }},{\mu }_{\star }=\sqrt{1-{\left(\frac{R_{\star }}{r}\right)}^2}.$$

Taking nowthe limit of short-wavelength unstable perturbations (kL ≫ 1), this becomes $$\frac{\delta g}{\delta v}\approx 2\Omega (1-s).$$(4)

A simple scaling relation for the above perturbed line-force then follows from applying the Sobolev theory in the optically thick limit for a single line $$\frac{g_{\text{Sob}}}{v_{\text{th}}}\propto \frac{L_{\star }}{r^2}\frac{(\text{d}v/\text{d}r)/\rho }{v_{\text{th}}}.$$(5)

Recalling that Ṁ ∝ r²ρ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) $$\begin{array}{ll}\frac{\delta g}{\delta v}\approx \Omega \hfill & \propto \frac{v_{\infty }^2}{v_{\text{th}}}\frac{L_{\star }}{R_{\star }}\frac{1}{\dot{M}},\hfill \\ \hfill & \propto v_{\infty }\left(\frac{v_{\infty }}{v_{\text{th}}}\right)\frac{T_{\text{eff}}^4{R}_{\star }}{\dot{M}},\hfill \end{array}$$(6)

where we have used $L_{\star }\propto R_{\star }^2{T}_{\text{eff}}^4$ 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 line-force.

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 small-amplitude ansatz the amplitude of the instability grows with exp Ωt with Ω⁻¹ a characteristic timescale (growth time). It can thus be understood that Ω⁻¹ 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 non-linear, 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 non-linear evolution of the radiatively unstable stellar wind we evolve the equations of 1D spherically symmetric isothermal radiation-hydrodynamics, 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 VH-1¹.

The key challenge in these instability simulations is a proper treatment of the radiation line-force in a time-dependent highly structured stellar wind outflow. As the instability acts on small spatial length scales, expensive non-local line-force integrations have to be carried out at each hydrodynamical time step. It is important that such non-local line-force computations should still closely approximate steady-state 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 non-uniform 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 ρ₀ is fixed andwe set velocity by constant slope extrapolation. The lower boundary density ρ₀ 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 line-driving

For our wind modelling we assume driving by an ensemble of lines described by a power-law 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 $$\frac{\text{d}N}{\text{d}q}=\frac{1}{\Gamma (\alpha )\overline{Q}}{\left(\frac{q}{\overline{Q}}\right)}^{\alpha -2}{e}^{-q/Q_{\text{max}}},$$(7)

with Γ(α) the Gamma function. We further describe the line-driven wind using α, the CAK power-law index denoting the relative contribution of the line-force from optically thick lines to the total line-force, and $\overline{Q}$ as a line-strength normalization factor describing the ratio of the total line-force to the electron scattering force if all lines were optically thin (Gayley 1995)². For the continuum opacities, we take, as is customary, an electron scattering contribution described by an opacity κ_e = 0.34 cm² g⁻¹. Electron scattering effectively lowers gravity and is commonly captured by Eddington’s gamma Γ_e ≡ κ_eL_⋆∕(4πGM_⋆c). The effective escape speed at the stellar surface then becomes $v_{\text{esc},\text{eff}}=\sqrt{2G{M}_{\star }(1-{\Gamma }_e)/R_{\star }}$. This set-up allows the computation of the average mass-loss 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 line-force is numerically computed.

4 Numerical results

4.1 Theoretical mass-loss rate predictions

It is important to point out that the bi-stability jump as predicted from line-driven wind models (Vink et al. 1999) implies an increase in mass-loss 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 non-linear simulations do not fundamentally predict global mass-loss rates. We chose our model wind line-force parameters (e.g. Puls et al. 2000) such that our relaxed CAK simulation mass-loss rates are approximately the same as those predicted from the Vink et al. (2001) recipe. Evolving such relaxed CAK models including the non-linear instability then provides average mass-loss 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 bi-stability jump region of Vink et al. (1999). The reason for this is as outlined above: our steady-state CAK models are required to converge to mass-loss rates close to those predicted by Vink et al. (2001) using a realistic set of wind parameters. However, the Vink et al. (2001) mass-loss rate recipe breaks down in the bi-stability 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, steady-state Sobolev model, and thus any line-driven instability simulation and wind clumping prediction within the bi-stability.

To that end our non-linear simulations do not provide a direct test of the predicted increase in mass-loss rates across the bi-stability 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 mass-loss rates derived from observational diagnostics. The observational lack of evidence for increased mass-loss rates might then be explained in terms of applying incorrect clumping factors for OB supergiants across the bi-stability jump.

4.2 Structure formation in O and B supergiant winds

We illustrate first the basic overall dynamics of an unstable 1D line-driven wind. Figure 2 shows the basic structure resulting from a non-linear 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 high-speed 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 steady-state 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 high-density 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 small-scale 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 Ω⁻¹ shows that the actual development of non-linear 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, re-occurring 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, limb-darkening 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 shallow-slope solutions that are reached in semi-regular 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 bi-stability jump

Having discussed the basic structures of the non-linear 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 $$f_{\text{cl}}=\frac{\langle {\rho }^2\rangle }{{\langle \rho \rangle }^2},$$(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 bi-stability jump. The average over the full wind already shows a significant difference between wind clumping properties on the hot and cool side of the bi-stability 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 mass-loss rates across the bi-stability jump (e.g. Markova & Puls 2008). However, this Hα fitting also depends sensitively on the wind clumping properties, due to the ρ² -dependence (see Introduction). Effectively, when we derive Ṁ from Hα fitting, we are really deriving the product $\dot{M}\sqrt{f_{\text{cl}}}$. This immediately shows that a degeneracy occurs when fitting line profiles because different combinations of Ṁ and f_cl can provide the same invariant $\dot{M}\sqrt{f_{\text{cl}}}$.

On the hot side of the bi-stability jump, our simulations show f_cl ~ 17 in the Hα forming region. On the other hand, the cool side of the bi-stability jump shows much lower factors f_cl ~ 2. We note that our predicted wind clumping factors on the hot side of the bi-stability 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 bi-stability 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 X-ray emission of B supergiant winds. A key observational result is that the X-ray flux drops in B-star winds as compared to their O-star 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 shock-heated regions directly. To facilitate a discussion around shocks, we define a velocity dispersion $$v_{\text{disp}}=\sqrt{\langle v^2\rangle -{\langle v\rangle }^2},$$(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 X-ray modelling, the velocity dispersion is a reasonable proxy for the overall intrinsic X-ray 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 small-scale velocity structures near the stellar base by the LDI. Eventually these small-scale 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⁻¹ 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⁻¹, which is enough toproduce intrinsic X-ray emission, whereas the velocity dispersion in a B supergiant model is several factors lower (v_disp ~ 70 km s⁻¹), and so will not result in much (or any) shock-heated gas hot enough to emit in the X-ray band. Physically speaking, this is a combined consequence of the weaker clumping properties and lower wind speeds of B supergiant line-driven 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 fcl = 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 bi-stability 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 time-dependent line-driven instability simulations that describe the non-linear evolution of the LDI.

Extending these line-driven instability simulations to B supergiants has allowed us to study the behaviour of wind clumping across the predicted bi-stability 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 bi-stability jump. Therefore, this suggests that accounting correctly for wind clumping is important in deriving empirical mass-loss rates across the bi-stability jump, if they are derived by means of the ρ² -sensitive Hα line. Specifically, in these diagnostic Hα studies of mass loss across the bi-stability jump, clumping effects are typically neglected entirely or treated as a constant across the bi-stability 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 mass-loss increase over the bi-stability 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 line-driven 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 luminosity-mass models would also exhibit significantly weaker wind clumping across the bi-stability jump.

Since the modelled velocity dispersion is much lower in BSG than OSG winds, our simulations provide a theoretical rationale to observations showing that X-ray fluxes drop significantly in single B-star winds compared to their O-star 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 high-mass X-ray binaries (HMXBs). These systems often consist of a B supergiant and an orbiting object that accretes material, for example the prototype system Vela X-1 where the donor star is a early B supergiant with a slow and dense wind characteristic of a simulation below the bi-stability jump (Sander et al. 2018). However, in the accretion models by El Mellah et al. (2018) targeting Vela X-1, a multi-D LDI simulation of an O supergiant was used to simulate the clumpy wind accretion and the corresponding time variation in X-ray luminosity. Our preliminary calculations suggest that also in such multi-D 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 follow-up work we plan to study in detail the effect of multi-D LDI models of B supergiants on the accretion rate and the corresponding variations of X-ray 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 bi-stability 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 bi-stability jump and discuss whether this may provide an explanation of the lack of observational evidence for increasing mass-loss rates across the bi-stability jump. A quantitative wind clumping description depends on the treatment of limb-darkening, photospheric base perturbations, the source function, and the amount of lateral fragmentation occurring in much more complex multi-dimensional wind clumping models (Dessart & Owocki 2003; Sundqvist et al. 2018).

Moreover, our simulations only make predictions of clumping factors, but not the overall stellar mass-loss rate. Extending this study to also treat stars within the bi-stability jump region must be left to future work once more reliable predictions for the smooth steady-state 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 bi-stability 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 bi-stability 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 EQUATION-group 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 line-force

As discussed in the Introduction, the Sobolev method is not suitable for computing the small-scale structures arising from the line-deshadowing instability as the instability acts on length scales smaller than the Sobolev length. This means that the line-force no longer depends only on local flow quantities; instead, there is a non-local coupling between the radiation field and outflowing gas. Such full non-local calculations are computationally daunting, and approximate expressions for the line-force 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 $$b_q(\mu ,r)\equiv {\displaystyle {\int }_{-\infty }^{+\infty }\text{d}}x\varphi \left(x-\mu v(r)/v_{\text{th}}\right)e^{-qt(x,p,z,z_b)},$$(A.1)

where ϕ is a profile function (taken to be Gaussian) at the observer’s frame frequency x ≡ (ν∕ν₀ − 1)c∕v_th defined from line-centre ν₀ 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: $$qt(x,p,z_1,z_2)\equiv {\displaystyle {\int }_{z_1}^{z_2}\kappa }\rho (r^{\prime })\varphi \left(x-\mu v(r^{\prime })/v_{\text{th}}\right)\text{d}{z}^{\prime }.$$(A.2)

Here we perform the integrations between two positions z₁ < z₂ along a certain ray p (impact parameter) from the origin where z ≡ μ∕r and r² ≡ p² + z². We note that κ in the above equation is not the usual opacity (expressed in units of [cm² g⁻¹]), but rather is a frequency integrated line-opacity (units [cm² g⁻¹ Hz⁻¹]). Specifically, κ ≡ χ∕ρ where χ is a frequency-integrated line-strength (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 (line-driven instabilities), this integral is inherently non-local 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 line-distribution (see Sect. 3.2, also definition in OP96) we define an escape probability $$\begin{array}{ll}b(\mu ,r)\hfill & \equiv \Gamma {(\alpha )}^{1/(1-\alpha )}{\displaystyle {\int }_0^{+\infty }\text{d}}q\frac{\text{d}N}{\text{d}q}\frac{q}{\overline{Q}}{b}_q(\mu ,r),\hfill \\ \hfill & =\Gamma {(\alpha )}^{(1+\alpha )/(1-\alpha )}{\displaystyle {\int }_0^{+\infty }\text{d}}q\frac{\varphi \left(x-\mu v(r)/v_{\text{th}}\right)}{{\left(\overline{Q}t(x,p,z,z_b)+\overline{Q}/Q_{\text{max}}\right)}^{\alpha }}.\hfill \end{array}$$(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 line-force expressions for the direct and diffuse radiation become $$g_{\text{dir}}=g_{\text{thin}}\frac{\langle \mu D(\mu ,r)b(\mu ,r)\rangle }{\langle \mu D(\mu ,r)\rangle }$$(A.4)

and $$g_{\text{diff}}=-2g_{\text{thin}}\frac{S(r)}{I_{\star }}\langle \mu b(\mu ,r)\rangle ,$$(A.5)

with $g_{\text{thin}}\equiv \overline{Q}{g}_{\text{e}}$ the line-acceleration in the optically thin limit and g_e = κ_eL_⋆∕(4πr²c) the acceleration due to continuum electron scattering. Angle brackets in the line-force expressions denote average angle integrations. The total force due to radiation g_rad is the sum of the above continuum force and line-forces (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 ${\mu }_{\star }=\sqrt{1-{(R_{\star }/r)}^2}$ 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 limb-darkening law (Hubeny & Mihalas 2014, Eq. (17.17)) $$D(\mu ,r)=\frac{1}{2}+\frac{3}{4}\sqrt{\frac{{\mu }^2-{\mu }_{\star }^2}{1-{\mu }_{\star }^2}},$$(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 line-force 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 $$S(r)=\langle I_{\star }D(\mu ,r)\rangle .$$(A.7)

Plugging Eq. (A.6) into the above then gives (Sundqvist & Owocki 2013) $$S^{\text{LD}}(r)=\frac{I_{\star }}{16}\left(7-4{\mu }_{\star }+3{\mu }_{\star }^2\frac{\ln \left[{\mu }_{\star }/(1+\sqrt{1-{\mu }_{\star }^2})\right]}{\sqrt{1-{\mu }_{\star }^2}}\right).$$(A.8)

For convenience it is useful to cast the angle integration of the line-forces in terms of a ray parameter $y\equiv {(p/R_{\star })}^2$ with direction cosine ${\mu }_y=\sqrt{1-y{(R_{\star }/r)}^2}$. Applying this ray formulation and the Eddington limb-darkened source function into the expressions for the line-forces (i.e. Eqs. (A.4) and (A.5)) gives $$g_{\text{dir}}(r)=g_{\text{thin}}\left(\frac{1}{2}{\displaystyle {\int }_0^{1}b}({\mu }_y,r)\text{d}y+\frac{3}{4}{\displaystyle {\int }_0^1\sqrt{1-y}}b({\mu }_y,r)\text{d}y\right)$$(A.9)

and $$g_{\text{diff}}(r)=\frac{S(r)}{S^{\text{LD}}(r)}\frac{g_{\text{thin}}}{2(1+{\mu }_{\star })}{\displaystyle {\int }_0^1\left(b_{-}({\mu }_y,r)-b_{+}({\mu }_y,r)\right)}\text{d}y,$$(A.10)

in which the diffuse line-force follows a “two-stream” approximation of inward (“−”, μ < 0) and outward (“+”, μ > 0) streaming photons. Additionally we correct with a factor ${(r/R_{\star })}^2$ as the ray integration should formally extend to $y={(r/R_{\star })}^2$ (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 frequency-dependent 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) $$t_{-}(-x,y,R_{\text{max}})=t_{+}(2{\mu }_{y}v(R_{\text{max}})/v_{\text{th}}-x,y,R_{\text{max}}).$$(A.11)

With the above equations, the line-forces are computed at each hydrodynamical time step and added together to get the total non-local line-force acting on the outflowing stellar wind.

References

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 fcl = 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 bi-stability 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