© F. A. Driessen et al. 2022

1 Introduction

The outflows and mass loss from hot, massive (OB) stars are predominantly caused by momentum transfer of photons scattering in spectral lines of ionised elements (Castor et al. 1975, hereafter CAK). This ‘line driving’ is very efficient because the Doppler effect allows a spectral line to sweep out (deshadow) and interact with continuum photons over a significant frequency range. Because the continuum photons interact with primarily metal atomic species, there is a direct dependence on the abundance and metallicity for both the wind driving and mass loss, which is confirmed by theoretical (steady-state) line-driven wind models (Vink et al. 2001; Gräfener & Hamann 2008; Krtička & Kubát 2018; Björklund et al. 2021; Vink & Sander 2021), as well as by observations (Mokiem et al. 2007; Tramper et al. 2011; Herrero et al. 2012; Ramírez-Agudelo et al. 2017; Marcolino et al. 2022).

However, line-driven winds are intrinsically unstable, because the deshadowing of wind ions from the Doppler effect leads to a runaway effect, or line-deshadowing instability (LDI; MacGregor et al. 1979; Carlberg 1980; Owocki & Rybicki 1984). In numerical simulations of the LDI, the wind is disrupted due to strong reverse shocks that lead to an inhomogeneous wind consisting of overdense clumps and a sparse interclump medium (Owocki et al. 1988; Feldmeier & Thomas 2017; Sundqvist et al. 2018; Driessen et al. 2019; Lagae et al. 2021).

This wind clumping then also severely affects the radiative transfer modelling needed to derive empirical mass-loss rates from comparisons to observed spectra. Indeed, neglecting clumping may lead to empirical mass-loss rate estimates that differ by more than an order of magnitude for the same star, depending on which observed spectral diagnostics is considered (Fullerton et al. 2006). Very recently, multi-diagnostic studies properly accounting for wind clumping (both optically thin and thick) were performed for a small set of O-supergiant stars in the Galaxy (Hawcroft et al. 2021), enabling first systematically derived constraints on statistical wind clumping properties. Nonetheless, wind clumping in OB-stars outside the Galaxy is still poorly investigated and constrained, both theoretically and observationally (although see Brands et al. 2022). This lack of constraints for clumping as a function of metallicity may then also indirectly affect the above-mentioned mass loss-metallicity relation for line-driven winds. For example, Mokiem et al. (2007) performed a large sample study of O-stars from the Galaxy and the Large and Small Magellanic Cloud (LMC and SMC) to infer empirical mass-loss rates. Using primarily the clumpingsensitive optical Hα line, these authors confirmed that LMC and SMC stars show lower mass-loss rates for a given stellar luminosity. However, Mokiem et al. (2007) did not take into account wind clumping corrections in their mass loss-metallicity relation, meaning that the quantity they effectively really derived was $\dot{M}\sqrt{f_{\text{cl}}}\propto Z^x$, with x a power-law index and f_cl the (quantitatively unknown) clumping factor (as defined, e.g. by Eq. (9) of this paper). With the advent of the Hubble Space Telescope Ultraviolet Legacy Library of Young Stars as Essential Standards (ULLYSES) project (Roman-Duval et al. 2020), large amounts of high-quality spectra of extragalactic massive stars in the ultraviolet (UV) will be added to already existing optical data. Using similar techniques to those in Hawcroft et al. (2021), we expect that these observations will be able to provide for the first time quantitative empirical constraints on the metallicity dependence of wind clumping. As such, it becomes pertinent that this then can be adequately compared to corresponding theoretical predictions stemming from state-of-the-art modelling of the unsteady line-driven wind. To date, all such theoretical investigations of line-driven wind clumping have been performed for Galactic stars, meaning the expected dependence on metallicity remains mostly unexplored.

In this paper, we investigate the dependence of wind clumping on metallicity from a theoretical point of view, using numerical radiation-hydrodynamic simulations of the LDI for O-stars at a fixed luminosity, but at different metallicities. The remainder of this paper is organised as follows: in Sect. 2 we detail our hydrodynamic setup, after which we present our results in Sect. 3. We discuss these in Sect. 4 and conclude in Sect. 5.

2 Method

We adopt the same numerical algorithms and grid setup as presented in Driessen et al. (2021, for specific details) within the open-source astrophysical (magneto-)hydrodynamic code MPI-AMRVAC (Xia et al. 2018; Keppens et al. 2021), except that we do not consider magnetic fields in the present work. This means we solve the equations of hydrodynamics on a two-dimensional (2D) spherical (r, θ) grid for density ρ and velocity v, $$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho v\right)=0,$$(1) $$\frac{\partial \left(\rho v\right)}{\partial t}+\nabla \cdot \left(\rho vv+pI\right)=\rho g_{\text{eff}}+\rho g_{\mathrm{line}},$$(2) $$\begin{array}{cc}p=a^2\rho ,& a=\sqrt{k_{\text{B}}{T}_{\text{wind}}/m}\end{array}.$$(3)

For the purpose of this study, we assume an isothermal wind at a temperature¹ T_wind ≈ T_eff with isotropic thermal pressure pI and isothermal sound speed a, which is a good approximation for the dense outflows from O-stars that undergo efficient radiative cooling (but see Feldmeier 1995; Lagae et al. 2021). In the momentum equation, we take into account the effect of stellar gravity reduced by constant electron scatterin Γ_e, g_eff = –GM_★(1 – Γ_e)/r² ê_r, and the radiation line force g_line. The adopted parameters in Table 1 also follow Driessen et al. (2021) and are representative of a high-luminosity Galactic O-star. We vary the line-strength parametrisation $\overline{Q}$ to simulate different metallicity conditions (see Sect. 3 for details). A key difference with Driessen et al. (2021) is that in this work, and for the first time in a 2D LDI wind model, we adopt a realistic line-strength cut-off $Q_{\max }=\overline{Q}$ (see Sect. 3). Choosing a fixed luminosity L_★ for our grid allows us to separate out another, potentially complicated relationship between clumping and L_★, and thus to focus exclusively on the Z dependence of interest here. The model luminosity L_★ has been chosen such that it coincides well with that considered by Mokiem et al. (2007) when they derived their empirical mass-loss metallicity relation for O-stars (neglecting the effects of clumping).

Following previous work on non-linear LDI simulations in a 2D hydrodynamic setting (Dessart & Owocki 2003; Driessen et al. 2021), we formulate the one-dimensional (1D) line force (consisting of direct and diffuse contributions) as ɡ_line = ɡ_line ê_r = (ɡ_dir + ɡ_diff) ê_r in the radial direction for each co-latitudinal cell using the smooth source function approximation (SSF; Owocki & Puls 1996). This avoids complex ray integrations from a bundle of rays with arbitrary directions. The adopted radial ray at each colatitude is then tilted with respect to the radial direction to mimic the finite extent of the stellar disc (Sundqvist & Owocki 2013) and to avoid that each point in the wind becomes a critical point (CAK). We refer the reader to Owocki & Puls (1996, their Sect. 5) and Driessen et al. (2021, their Eqs. (9)–(14)) for an extensive account of the implementation of ɡ_dir, ɡ_diff into the hydrodynamics.

Inner and outer radial boundary conditions on the hydrodynamic variables are set in the same fashion as in Driessen et al. (2021), except that for the radiation field emerging from the star we do not assume a uniformly bright disc, but apply an Eddington limb-darkening law, $$D\left(\mu ,r\right)=\frac{1}{2}+\frac{3}{4}\mu \text{'}=\frac{1}{2}+\frac{3}{4}\sqrt{\frac{{\mu }^2-{\mu }_{\star }^2}{1-{\mu }_{\star }^2}},$$(4)

with μ' = cos θ' the angle between a star-centred radius vector and a ray emanating from the stellar surface μ, the impact parameter, and μ_★ the dilution factor. This modifies at the lower radial boundary the direct radiation field ɡ_dir contribution (Sundqvist & Owocki 2013; Driessen et al. 2019, for more details). To avoid numerical complications at the poles, all hydrodynamic variables at the (singular) polar axes are set using π-boundary conditions to connect both poles and allow flow from one side to the other without an artificial boundary.

As the radiation transport is only solved for in the radial direction (with a fiducial tilt, see above) this means that nonradial radiation effects are inhibited within our simulations. In particular, we are ignoring the line-drag effect that linear analyses suggest may damp lateral velocity fluctuations on scales below the lateral Sobolev length (Rybicki et al. 1990; Driessen et al. 2020). Recent 2D radiation transport simulations of the LDI by Sundqvist et al. (2018) (albeit on a spatially restricted grid) indeed show that typical lateral structures seem to be on the order of the Sobolev length (~0.01R_★ at r ~ 2 R_★). This then motivates our choice here of a lateral grid, ensuring that a typical lateral Sobolev length (or arc length angle), L_θ/r ≡ v_th/v_r ≈ 0-5°, is resolved within our simulations.

The simulations start from a 1D, spherically symmetric, steady-state CAK wind that applies the Sobolev approximation in the computation of ɡ_line. The corresponding initial smooth density and radial velocity profiles are distributed over the meridional plane while setting lateral velocities to zero. From this initial state, we evolve the radiation-hydrodynamic equations using the SSF formalism for a total physical time of approximately 300 ks (corresponds to about 30 wind dynamical timescales: τ_dyn = R_★/v_∞ ≈ 10 ks).

3 Simulation results

In this work, we do not rely on the original CAK line-force parametrisation but use instead the conceptually advantageous but equivalent $\overline{Q}$ parametrisation proposed by Gayley (1995). For our purposes here, a key advantage of the line strength $\overline{Q}$ is its direct proportionality to metallicity Z (Gayley 1995, his Eq. (60); Puls et al. 2000, their Eq. (81)), $$\overline{Q}\propto Z.$$(5)

A reference value is chosen such that the line-strength parameter takes the standard value $\overline{Q}=2000$ found for O-stars at solar metallicity Z_⊙ (Gayley 1995; Puls et al. 2000). In order to then investigate how LDI-generated structure depends on Z, we simply vary $\overline{Q}$ accordingly; that is, for example a simulation at $\overline{Q}=1000$ approximately corresponds to a metallicity Z/Z_⊙ = 0.5 (LMC) and one at $\overline{Q}=400$ to Z/Z_⊙ = 0.2 (SMC).

A priori, the direct scaling $\overline{Q}$ of purely with Z might appear surprising. However, the ensemble line strength $\overline{Q}$ is only dependent on the properties of the radiation field, w_Vi = F_viv_i/F, for flux F and frequency of transition Vi, and the strength of a given single spectral line q_i (Gayley 1995, his Eq. (14); Puls et al. 2000, their Eq. (37)), $$\overline{Q}\equiv {\displaystyle \sum _{i=1}^{\text{all\hspace{0.17em}lines}}{w}_{v_i}{q}_i}.$$(6)

The latter here depends on the atomic level occupation number density ni (neglecting stimulated emission), $$q_i\equiv \frac{{\sigma }_{\text{cl}}f{n}_l}{{\sigma }_{\text{Th}}{v}_i{n}_e},$$(7)

with f the oscillator strength, σ_cl = πe²/(mec) the classical frequency-integrated line cross-section for an electron with charge e and mass me, σ_Th the Thomson cross-section, and ne the electron number density. As the vast majority of driving elements are metals in O-stars, it becomes clear from this definition that q_i ∝ n_i ∝ A ∝ Z, that is q_i (and thus $\overline{Q}$) directly relates to the global metallicity Z via the elemental abundance A. On the other hand, for cooler A-stars (T_eff ~ 10 000 K), contrary to the hotter O-stars studied in this work, line-driving from hydrogen lines might also be important, meaning that this simple linear relation between $\overline{Q}$ and Z may not be completely valid in that regime (Puls et al. 2000, see discussion below their Eq. (81)). Therefore, determining the exact dependence of the $\overline{Q}\left(Z\right)$ relation for any given effective temperature is a topic for future research.

Fig. 1 Spatial and temporal averaged mass-loss rate <Ṁ> as a function of line strength found from our simulations (dots). The dashed line is a power-law fit to the wind models.

3.1 Model calibration and dynamics of the structured wind

As we keep all stellar parameters except the metallicity fixed, our simulations should recover some standard scaling relations for the mass-loss rate as a function of metallicity (or more precisely as a function of our metallicity proxy $\overline{Q}$; see above). To probe the consistency of our wind models in this respect, we closely monitor the mass-loss rate in the different simulations. It is important to stress here that we do not predict mass-loss rates from our LDI wind models. On the contrary, we choose our lineforce parameters $\left(\alpha ,\overline{Q}\right)$ such that the corresponding 1D relaxed CAK simulation provides a reasonable mass flux for a typical O-star at the model luminosity. Performing non-linear LDI simulations then results in a time-averaged mass-loss rate 〈Ṁ〉 that should approximate this mass flux.

Figure 1 shows the variation of 〈Ṁ〉 with line strength. As expected, the mass-loss rate decreases with decreasing line strength because of the weaker wind. The decrease in mass loss between different line strengths can be compared to the theoretically expected reduction following analytic CAK theory. We take as reference the mass-loss rate ${\dot{M}}_{{\overline{Q}}_{2000}}$ for a Galactic star at solar metallicity, whereby:

$${\dot{M}}_{\overline{Q}}={\left(\frac{\overline{Q}}{{\overline{Q}}_{2000}}\right)}^{\left(1-\alpha \right)/\alpha }{\dot{M}}_{{\overline{Q}}_{2000}}.$$(8)

As all parameters are fixed within our model, except for $\overline{Q}$, the mass-loss rate at a given metallicity and fixed α = 0.65 should scale as $\dot{M}\propto {\overline{Q}}^{0.53}$. When performing a power-law fit of the dependency of mass loss on line strength in Fig. 1, we indeed recover a power law with an index of 0.58. Given the non-linear nature of our simulations, the variation from the analytical CAK value is well within the acceptable range. Moreover, our $\overline{Q}=2000$ reference model yields log₁₀ 〈Ṁ〉 = −5.61, which also compares reasonably well with the detailed steady-state predictions by Björklund et al. (2021) for a star of the same luminosity and Z = Z_⊙ (log₁₀ Ṁ = −5.75). We do note that Björklund et al. (2021) formally include millions of individual spectral lines to compute the radiative acceleration in their 1D stationary models. Within such detailed time-independent models, it takes typically ≈30 minutes to converge the non-local thermodynamic equilibrium radiative transfer and compute the line force for one hydrodynamical iteration (a converged wind model then typically needs ≈50 such iterations). Following a similar approach within a time-dependent code is computationally unfeasible because it requires solving the formal integral of radiative transfer for every relevant spectral line at every radial grid point and at each hydrodynamic time-step. For example, a typical simulation presented here needs ≈1.5 × 10⁵ iterations to march towards t_end ≈ 300 ks, which shows the impractical nature of such detailed line-transfer calculation for the radiative acceleration². Instead, the line-ensemble parametrisation adapted here avoids such complications and statistically characterises the cumulative effect of all lines with an exponentially truncated power-law distribution and parameters $\left(\alpha ,\overline{Q},Q_{\max }\right)$. As the obtained reference mass-loss rate agrees well with more detailed computations, we consider our wind models to be well-calibrated with respect to the chosen Galactic standard model of $\overline{Q}=2000$, such that our further predictions with respect to the wind clumping should be robust in a relative sense.

Figure 2 shows the resulting wind density at a snapshot taken well after (at 300 ks) the LDI simulation has developed its characteristic structure consisting of high-density clumps and rarefied interclump gas. Specifically, we compare the Galactic wind model $\overline{Q}=2000\left(Z\approx Z_{\odot }\right)$ with a model mimicking the metallicity of the LMC, $\overline{Q}=1000\left(Z=0.5Z_{\odot }\right)$ and with one mimicking the SMC, $\overline{Q}=400\left(Z=0.2Z_{\odot }\right)$. All simulations display the characteristic LDI-generated features of a rather disrupted flow with quite slowly accelerated wind clumps that are separated by fast rarefied interclump gas. The lateral break-up of the LDI-generated structure is within these SSF models due to shearing motions together with Rayleigh-Taylor and thin-shell instabilities (Dessart & Owocki 2003) because the LDI itself only operates in the radial direction. Moreover, across this metallicity range, the wind clumps remain of similar size. However, this does not imply that the statistically computed overdensities of the clumps are the same (see Sect. 3.2). Within the 1D SSF approximation, on what spatial scales LDI-generated clumps are expected to form is still unclear. Fragmentation might result in relatively large scales due to the lateral line-drag (not present in these models, Rybicki et al. 1990; Driessen et al. 2020) or photospheric turbulent structures interacting with the LDI may also influence the clump spatial scales in the wind (see discussion in Sect. 4). Nonetheless, within the present model assumptions, the wind has a high degree of lateral incoherence whereby some clumpy structures only span a spatial scale close to the lateral cell size (i.e. they are on the order of the lateral Sobolev length).

The effect of lateral break-up is likely enhanced here with our adaptation of a realistic line-strength cut-off: $Q_{\max }=\overline{Q}$ (Gayley 1995). This cut-off was originally introduced by Owocki et al. (1988) to modify the CAK line-ensemble distribution and avoid an infinite radiative acceleration in the limit of optical depths approaching zero. However, in LDI simulations hitherto performed the cut-off has still been lowered to an artificial value of $Q_{\max }\approx {10}^{-3}\overline{Q}$ to avoid severe numerical issues due to the overly steep acceleration of gas (Feldmeier 1995; Sundqvist & Owocki 2013; Sundqvist et al. 2018; Driessen et al. 2019, 2021). In our new 2D hydrodynamic simulations presented here, we are finally able to increase this cut-off to a physically realistic value, $Q_{\max }=\overline{Q}$. We speculate that the previous artificially low cut-off was inherent to the 1D LDI simulations that form the majority of the models present in the literature. Indeed, in these 1D models, gas parcels are forced to run into each other when the relative accelerations become too strong. The resulting strong shocks lead to strong density compressions, but also extremely rarefied gas that would get accelerated even more by virtue of its low optical depth, such that an artificially low line-strength cut-off is a necessity. These problems are alleviated in 2D models such as those presented here, where gas parcels are also able to flow past each other. As a result, some strongly accelerated rarefied gas is able to overtake an overdense clump instead of running into it and causing an even stronger shock. This then leads to an overall smoothing effect of the structures, which in turn means that we no longer have to cut off the line strength in order to prevent further growth. Hence, these new 2D simulations with $Q_{\max }=\overline{Q}$ should also represent a significant step forward regarding quantitative evaluation of clumping structures caused by the LDI (albeit still with the caveats discussed above concerning the lateral coherence scales).

Fig. 2 Comparison of logarithmic density variations with respect to the mean wind density between a Galactic metallicity $(\overline{Q}=2000)$, LMC metallicity $(\overline{Q}=1000)$, and SMC metallicity $(\overline{Q}=400)$ stellar wind.

3.2 Spatial dependence of wind clumping

During our simulations we compute several statistical variables (i.e. temporally averaged quantities) to characterise the wind. These statistical quantities are computed upon every iteration once the wind has sufficiently developed that it now displays its characteristic structure. Typically, this is reached within our simulation after a time of t_stat = 200ks, or about 20τ_dyn, such that about 100 ks in our simulations is devoted to computing the temporal averages.

In this paper, we focus on the wind clumping factor, $$f_{\text{cl}}=\frac{\langle {\rho }^2\rangle }{{\langle \rho \rangle }^2},$$(9)

to quantify the amount of overdensity with respect to the mean wind density. As our models are 2D, we have a wind clumping factor in each point of the meridional plane. To that end, we compute from each model a corresponding lateral-averaged (1D) wind clumping profile in order to separate out the primary dependence on radius from the more stochastic lateral variation; the resulting f_cl(r) is displayed in Fig. 3 for all simulations. The models confirm that, at lower metallicity, not only does the wind driving at a fixed luminosity become weaker, but so do the reverse shocks and density compressions that are caused by the LDI.

We further find that, although structure formation generally starts already at the wind base, the growth in the inner wind is somewhat slower than in previous 1D models. To rule out that differences in wind clumping are due to the use of different codes (mpi-amrvac here; VH-1 in previous studies³), we compare in Fig. 4 the wind clumping stratification stemming from a standard 1D model with an artificial cut-off of $Q_{\max }={10}^{-3}\overline{Q}$, a similar 2D model with $Q_{\max }={10}^{-3}\overline{Q}$, and the 2D model with $Q_{\max }=\overline{Q}$ presented in this work, all for the reference line strength $\overline{Q}=2000$ and all computed with mpi-amrvac. We find a more modest degree of wind clumping for the 2D models in the region r ≤ 2R_★ as compared to 1D models (Sundqvist & Owocki 2013; Driessen et al. 2019). Nonetheless, our 1D model computed with mpi-amrvac agrees well with previous 1D models such that we cannot attribute the different wind clumping stratifications to the use of different codes. Specifically, in those 1D simulations of high-luminosity O-stars, the degree of wind clumping is found to strongly peak at ≈ 2R_★ with fcl ~ 20–30 (Sundqvist & Owocki 2013; Driessen et al. 2019) after which it starts to decrease significantly towards the outer wind. By contrast, from our new 2D simulations we find somewhat lower absolute values for the clumping factor as well as a less pronounced peak (see Figs. 3, 4). Again we suspect that the same physics that governs the effect of the line-strength cut-off Q_max (see Sect. 3.1) manifests itself also in setting the absolute degree of wind clumping. indeed, the 2D nature of our models does not necessarily result in local gas-clump collisions that would lead to stronger shocks and enhanced density contrasts (thus higher wind clumping), and with an artificial cut-off to almost no degree of clumping. Ina1D model, this is a natural outcome such that 1D LDi wind models generally predict higher fcl than corresponding 2D models (see also Dessart & Owocki 2003; Sundqvist et al. 2018).

Fig. 3 Radial stratification (away from the star) of the wind clumping factor fcl for each metallicity model.

Fig. 4 Radial stratification (away from the star) of the wind clumping factor fcl for the reference model $\overline{Q}=2000$ in different dimensions. The solid curve is a 1D simulation with $Q_{\max }={10}^{-3}\overline{Q}$ as found in the literature so far (but now computed with MPI-AMRVAC), while the dashed curve shows a comparable 2D simulation with $Q_{\max }={10}^{-3}\overline{Q}$. The dotted curve displays our novel 2D simulation with $Q_{\max }=\overline{Q}$. The stratification curves from the 2D models appear smoother because they are laterally averaged to retrieve a single radial stratification for comparison with a purely 1D model.

Fig. 5 Spatially- and temporally averaged wind clumping factor 〈fCl〉 as a function of line strength found from our simulations (dots). The dashed line is a power-law fit to the wind models while the gray-shaded area gives the standard deviation σ = 0.01 on the slope.

3.3 Power-law dependence of wind clumping and metallicity

The results shown above (Fig. 3) already suggest that the amount of wind clumping decreases somewhat with decreasing metallicity. Such an effect might indeed be expected, because a lower metallicity wind will undergo weaker wind driving that also results in weaker amplification of velocity fluctuations and subsequent shock formation. We here quantify to what degree wind clumping is dependent on metallicity by assuming a simple power-law relation between fcl and $\overline{Q}\left(\propto Z\right)$.

Thus, we assume, $$f_{\text{cl}}=a{\overline{Q}}^b\iff {\log }_{10}{f}_{\text{cl}}={\log }_{10}a+b{\log }_{10}\overline{Q},$$(10)

which represents a straight line in logarithmic space. The constant a gives the intercept with the ordinate axis and the powerlaw index b gives the slope of the straight line. To obtain one wind clumping estimate for the total wind, we compute from our 2D simulations the resulting lateral-averaged wind clumping (e.g. Fig. 3) after which we further average this over the radial direction to get a characteristic 〈fcl〉 for the model. The results of this procedure are displayed in Fig. 5 for each metallicity model, showing a clear increasing trend for h 〈fcl〉 with increasing $\overline{Q}$.

We apply a standard linear regression routine available from the Python scipy library (Virtanen et al. 2020) to fit the simulation results to a power law with the above-mentioned form. This yields: $$f_{\text{cl}}\propto {\overline{Q}}^{0.15\pm 0.01}\propto Z^{0.15\pm 0.01}.$$(11)

The resulting fit therefore suggests a clear but still rather weak dependence of wind clumping on metallicity.

4 Discussion

As mentioned in Sect. 2, a key reason for our choice of characteristic model luminosity was the study by Mokiem et al. (2007), who derived their empirical mass-loss-metallicity relation for stars in the same luminosity range. More specifically, Mokiem et al. (2007) performed an extensive survey of O-type stars in the Galaxy, LMC, and SMC using the optical Hα line, and derived an empirically calibrated mass-loss rate-metallicity relation: $$\dot{M}\propto Z^{0.83\pm 0.16}.$$(12)

However, this relation was obtained by means of spectral fitting assuming smooth winds. As Hα is a recombination line in O-stars, its emission will be very sensitive to the amount of clumping present in the wind. Namely, in an inhomogeneous wind, the mean opacity (per unit length) for Hα scales as 〈 x〉 ~ 〈ρ²〉 ~ f_clṀ² for an assumed mass-loss rate Ṁ ~ 〈p〉. That is, the scaling invariant is no longer proportional to Ṁ but rather to $\dot{M}\sqrt{f_{\text{cl}}}$. As such, the empirical relation above should be modified according to: $$\dot{M}\sqrt{f_{\text{cl}}}\propto Z^{0.83\pm 0.16}.$$(13)

Inserting our theoretically derived wind clumping dependence on metallicity into this relation then finds a marginal correction to the mass loss-metallicity dependence: $$\dot{M}\propto Z^{0.83}{Z}^{-0.15/2}=Z^{0.76}.$$(14)

Within uncertainties, this agrees fairly well with predictions made from steady-state line-driven wind models that do not rely on any underlying assumptions about the line-ensemble parameters when deriving the mass-loss rates (Vink et al. 2001; Krtička & Kubát 2018; Björklund et al. 2021). As mentioned in the previous section, we further find somewhat lower degrees of clumping than corresponding 1D models in the inner wind, r ≳ 2R_★, and also do not observe any prominent peak for the clumping factor within our simulated domain. Using multi-wavelength diagnostic studies Puls et al. (2006), Najarro et al. (2011), and Rubio-Díez et al. (2022) find that the inner wind of O-stars in the luminosity range considered here seem to be rather more clumped than the outermost radio-emitting wind. However, we note that the radio photosphere of such O-stars typically lies well beyond the outer boundary radius (r/R_★ = 6) of our 2D simulations. In Galactic 1D LDI models with much larger outer boundary radii, f_cl ≈ 4 has been found for the radio emission region at r ≳ 50 R_★ (Runacres & Owocki 2002); this is indeed lower than the 〈f_cl〉 ≈ 12 derived here from our Galactic 2D simulation in the range r/R★ = 1–6, which is in qualitative accordance with the observational findings mentioned above.

Nevertheless, some recent studies focusing on multidiagnostic spectral line fitting in the UV and optical found average clumping factors f_cl ≈ 25 for O-supergiants (Hawcroft et al. 2021; Brands et al. 2022). Quantitatively, this is higher than the average values derived here, and may tentatively point towards some still missing physics in our simulations. In particular, considering that the 1D → 2D transition has led to overall lower degrees of clumping, it seems unlikely that, for example, a full three-dimensional simulation would reverse this trend. Alternatively, we might expect that adding a third dimension to the simulations would further decrease the quantitative levels of clumping, because then the high-density gas parcels will have yet another dimension to escape through that may lead to even more smoothing effects.

In this respect, relaxing our assumption of a hydrostatic boundary at the stellar surface that connects to the wind might aid in inducing such higher levels of clumping also in the innermost wind. Specifically, for high-luminosity Galactic O-stars, the radiative acceleration might exceed gravity even in deep subsurface layers (at the so-called “iron-opacity bump” at T ≈ 200 kK), and so it is possible that this could trigger a very turbulent lower atmosphere. Recently, Jiang et al. (2015, 2018) performed three-dimensional radiation-hydrodynamic simulations of massive-star envelopes and found that velocity and density fluctuations originating from the deep stellar layers are so vigorous that they can protrude into the stellar photosphere (Schultz et al. 2020). Although these “blobs” can gain significant upward momentum from the subsurface opacity peaks, they typically stagnate and lead to ‘failed winds’ falling back down below the surface, simply because these first simulations of subsurface layers relied on Rosseland mean opacities calibrated for static media and thus not including any line-driving effects. Recently, our group implemented suitable modules into MPI-AMRVAC for computing the relevant radiative energy and momentum terms (Moens et al. 2022), and further developed a hybrid opacity scheme able to simultaneously model the deep subsurface layers and a line-driven supersonic outflow in one unified approach (Poniatowski et al. 2021, and in prep.). With this, we will be able to investigate the effects of opacity-driven convective blobs on the line-driven wind structure and clumping properties in OB-type stars in future work.

However, in this paper our primary interest has not been to obtain any absolute wind clumping factors, but rather to investigate their dependence on metallicity in a relative sense. As such, we consider our main simulation result to be robust, namely that the overall wind clumping decreases with metallicity for an O-star with fixed luminosity. In particular, effects stemming from subsurface layers will decrease with metallicity (primary opacity peak due to iron; see above), and so this should not change this overall qualitative behaviour.

Nonetheless, extensions to 3D and inclusion of a turbulent photosphere may slightly modify the quantitative behaviour of the clumping power law displayed in Fig. 3. However, in this respect, we note again that our 2D models are calibrated with respect to a standard Galactic reference model and so should be relatively robust. Proper consideration of the above-discussed effects may therefore change the wind clumping stratification in Fig. 3, but because we compute average wind clumping factors, any such difference in stratification will presumably be reflected in each model average. In other words, we expect at least the qualitative behaviour of wind clumping on metallicity to be relatively insensitive to further improvements of the underlying physics in radiation-hydrodynamic simulations of line-driven O-star winds.

Finally, the models considered in this paper are all taken at a fixed reference luminosity in order to establish a first the oretical f_cl(Z) dependence. Following this approach allows a direct comparison with the empirical mass-loss study of Mokiem et al. (2007), who considered a similar reference luminosity, as discussed above. Nonetheless, it is likely that varying the stellar luminosity would also affect the amount of wind clumping. In the linear analysis of Driessen et al. (2019), it was already found that the growth rate of the LDI (these authors did not consider metallicity effects in their simplified analysis) shows a complex dependence on stellar and wind parameters, and thus the amount of clumping in the non-linear regime may also depend on this indirectly. Yet, a potential luminosity dependence would add another free parameter to the current analysis such that it becomes difficult to disentangle the basic metallicity dependency we aim to investigate here. Whether the f_cl(Z) dependence derived above would also still hold for (significantly) different L_★, or would show a more complex dependence on stellar parameters therefore remains an open question.

5 Conclusions and outlook

We present 2D radiation-hydrodynamic simulations of the LDI using the SSF approximation to evaluate the 1D line force (Owocki & Puls 1996). An important improvement of these 2D models is that we are able to, for the first time, assume a physically realistic line-strength cut-off when evaluating the radiation line force.

In particular, we used our new simulations to investigate whether or not there is a dependence of wind clumping on metallicity for O-stars with a fixed luminosity. Our simulations suggest that such a dependence indeed exists, albeit a rather weak one; fitting a power law to our models yields f_cl ≈ Z^0-15. Therefore, empirically inferred mass-loss rate-metallicity relations derived from diagnostics that probe $\dot{M}\sqrt{f_{\text{cl}}}$ will only be modestly affected by our predicted f_cl(Z) relation; for example, the empirical Ṁ ∝ Z^0.83 derived by Mokiem et al. (2007) would be shifted to Ṁ ∝ Z^0.76.

The fcl(Z) relation derived in this work can be tested with upcoming spectra acquired as part of the ULLYSES project (Roman-Duval et al. 2020). In particular, following Hawcroft et al. (2021) and Brands et al. (2022), novel model atmosphere calculations combined with multi-diagnostic genetic algorithm fitting for large samples of stars should be able to put firm constraints on the empirical fcl(Z) dependence. This would then allow us to assess whether our current 2D LDI simulations require further extensions and/or more input physics (see discussion in previous section), or if our predicted behaviour is robust against such further model improvements.

Additionally, such multi-diagnostic spectral fitting will also be able to establish better constraints on absolute wind clumping factors. This is of particular importance because our current 2D LDI wind models, contrary to previous 1D models, seem to quantitatively predict somewhat lower levels of fcl than inferred from recent studies of Galactic O-supergiants (Hawcroft et al. 2021).

We further plan to re-examine our models by relaxing the assumption of a hydrostatic lower atmosphere in order to study the effect of opacity-driven convective and turbulent motions originating in the massive-star outer envelope protruding into the lower atmosphere (Jiang et al. 2015, 2018). The combined effect of such a time-dependent lower boundary and the interplay of the resulting ‘blobs’ with the LDI can in principle be studied with our group’s recently developed hybrid opacity-method (Poniatowski et al. 2021, and in prep.) and new radiationhydrodynamic module in MPI-AMRVAC (Moens et al. 2022). This will then allow us to assess how such subsurface motions may also affect quantitative predictions of wind clumping, in particular in the inner wind. Furthermore, in published time dependent line-driven wind models, the parameters $\left(\alpha ,\overline{Q},Q_{\max }\right)$ (or different, but equivalent parametrisations) describing the driving from an ensemble of lines is taken as fixed both in space and time. To that end, here we focused our attention on a typical high-luminosity, early O-star for which the line-ensemble parameters are quite well constrained (Gayley 1995; Puls et al. 2008). However, this fixing of the line-ensemble parameters for the full wind is a stringent condition that may have a significant impact on the wind dynamics. In future work, we plan to update and improve our LDI wind models (and clumping predictions) using tabulated line-driven wind parameters from detailed line-list computations that self-consistently compute $\left(\alpha ,\overline{Q},Q_{\max }\right)$ based on the local hydrodynamic conditions of the wind (Poniatowski et al., in prep.). This would then also allow a larger systematic and self-consistent study of LDI-generated structures as a function of stellar parameters (L_★, M_★, R_★).

In summary, we consider the 2D models presented here to be the first step towards a more quantitative clumping analysis using numerical line-driven wind simulations at different metallicities. As discussed above, in that respect our study opens up many new avenues for the theoretical and observational study of structure formation in the winds of massive stars.

Acknowledgements

FAD and JOS acknowledge support from the Odysseus program of the Belgian Research Foundation Flanders (FWO) under grant G0H9218N. We thank the anonymous referee for giving valuable comments that helped to clarify this work in several aspects.

References

All Tables

All Figures

	Fig. 1 Spatial and temporal averaged mass-loss rate <Ṁ> as a function of line strength found from our simulations (dots). The dashed line is a power-law fit to the wind models.
In the text

	Fig. 2 Comparison of logarithmic density variations with respect to the mean wind density between a Galactic metallicity $(\overline{Q}=2000)$, LMC metallicity $(\overline{Q}=1000)$, and SMC metallicity $(\overline{Q}=400)$ stellar wind.
In the text

	Fig. 3 Radial stratification (away from the star) of the wind clumping factor fcl for each metallicity model.
In the text

Fig. 4 Radial stratification (away from the star) of the wind clumping factor fcl for the reference model $\overline{Q}=2000$ in different dimensions. The solid curve is a 1D simulation with $Q_{\max }={10}^{-3}\overline{Q}$ as found in the literature so far (but now computed with MPI-AMRVAC), while the dashed curve shows a comparable 2D simulation with $Q_{\max }={10}^{-3}\overline{Q}$. The dotted curve displays our novel 2D simulation with $Q_{\max }=\overline{Q}$. The stratification curves from the 2D models appear smoother because they are laterally averaged to retrieve a single radial stratification for comparison with a purely 1D model.

In the text

	Fig. 5 Spatially- and temporally averaged wind clumping factor 〈fCl〉 as a function of line strength found from our simulations (dots). The dashed line is a power-law fit to the wind models while the gray-shaded area gives the standard deviation σ = 0.01 on the slope.
In the text