© L. G. Poniatowski et al. 2022

1 Introduction

Radiation plays a dominant role in the dynamics of various astrophysical environments. In addition to exchanging energy with matter via heating and cooling processes, radiation also exchanges momentum with the medium through absorption, emission, and scattering processes. Such momentum provided by radiation is dynamically important within, for example, accretion discs around stars, in active galactic nuclei, in cataclysmic variables, and in the atmospheres and winds of stars (Puls et al. 2008; Higginbottom et al. 2014; Kee et al. 2016; Jiang et al. 2019).

In massive-star atmospheres, the outward push from radiation causes the star to lose mass through powerful radiation-driven stellar winds. Prominent spectral features such as P-Cygni lines provide observational evidence of radiation-driven wind outflows (for a review, see Puls et al. 2008; Vink 2022), which Lucy & Solomon (1970) theorised are mainly powered by absorption and scattering within spectral lines. For a given atom, there can be a vast number of spectral lines for photons to interact with, and in accelerating and transonic media, these effects are further enhanced due to Doppler shifts, leading to strong radiation line forces (see e.g. Puls et al. 2008; Owocki 2015; Sander et al. 2020). However, accurately accounting for such ‘line driving’ in simulations of supersonic flows is a very challenging and computationally expensive problem.

Indeed, state-of-the-art radiation-hydrodynamic (RHD) applications always rely on various approximations when computing the acceleration due to spectral lines. On the one hand, applications where the line force is computed more accurately (i.e. by solving the multi-line radiative transfer problem in a reference frame), have exclusively focused on 1D steady-state cases (e.g. Sander et al. 2017; Sundqvist et al. 2019). On the other hand, multi-dimensional, time-dependent studies either entirely neglect the dynamic effects of line driving in supersonic media, that is, they only consider opacities calibrated for static media such as Rosseland mean opacities (e.g. Jiang et al. 2018, 2019; Goldberg et al. 2021), or invoke further approximations (e.g. Proga et al. 1998; Petrenz & Puls 2000; Sundqvist et al. 2018).

In Poniatowski et al. (2021), we proposed a hybrid model wherein Rosseland mean opacities were combined with opacities that account for line driving in supersonic regions (as also suggested for the ‘velocity-stretch’ opacity methods discussed by Castor 2004, chapter 6). For the line driving, the hybrid opacity model by Poniatowski et al. (2021) relied on a parametrisation first proposed by Castor et al. (1975, hereafter CAK). They used the so-called Sobolev approximation (Sobolev 1960), assuming that hydrodynamic quantities¹ are constant over the photon-line interaction region and approximated the accumulative acceleration from an ensemble of lines as a power law of the local projected velocity gradient. CAK computed fitting parameters to their assumed power law for a set of C III lines assuming local thermodynamic equilibrium (LTE), and several later studies then updated these fits using more complete lists of spectral lines (e.g. Abbott 1982; Puls et al. 2000). However, these approaches typically only provide the CAK fit parameters as a function of stellar properties, thus limiting their range of potential applications to (most often 1D and steady-state) simulations of hot-star winds (see also Gormaz-Matamala et al. 2019).

In time-dependent numerical RHD simulations based on a CAK parameterisation, line-force parameters are typically assumed to take some pre-described values; for example, in Poniatowski et al. (2021) an ad hoc radial stratification for the CAK power-law index parameter α was assumed. However, in a more realistic treatment, the parameters describing the line force should rather be self-consistently updated within the simulation itself according to the local gas and radiation conditions. A key advantage of our proposed method, compared for example to the similar tabulations by Lattimer & Cranmer (2021, hereafter LC), is that it also offers direct insight into the underlying linedistribution function. Since most general RHD codes require an opacity as input (see e.g. Davis et al. 2012; Moens et al. 2022b), we further chose to work directly with (flux-weighted) opacities here so that our resulting tables can be used directly in such RHD codes in a way quite analogous to how Rosseland mean opacities (e.g. ‘OPAL’, Iglesias & Rogers 1996) are used in the static limit.

For example, two applications already using our new method regard multi-dimensional extensions of the 1D Wolf-Rayet (WR) models presented in Moens et al. (2022b) and an extension of the current method to corresponding non-Sobolev, observer’s frame line-force calculations (essentially following Owocki & Puls 1996). This latter will allow us, for the first time, to investigate effects of self-consistent line-force parameters on the strong ‘line deshadowing instability’ (LDI) (Owocki & Rybicki 1984); thus far, all such LDI simulations have simply used pre-described, prototypical line-force parameters that are assumed to be fixed in time and space (e.g. Owocki et al. 1988; Sundqvist et al. 2018; Driessen et al. 2021).

In this paper, we use the Munich atomic database (Pauldrach et al. 1998, 2001; Puls et al. 2005), which contains over 4 million spectral lines, to calculate the flux-weighted line opacity (Sect. 2). We then fit our results from the complete line list to an expression derived from the assumed line-strength distribution function, tabulating them by means of our (three) fitting parameters (Sect. 2). To benchmark our method, we apply our tables in a 1D time-dependent RHD showcase study of the calibrated grid of line-driven O-star models from Björklund et al. (2021), comparing the resulting wind outflows to predictions for this regime from the literature (Sect. 3). In the final section, we discuss our results and give conclusions (Sect. 4).

2 Flux-Weighted Mean Opacity from Spectral Lines in an Accelerating Supersonic Medium

As mentioned above, using a statistical approach, one may approximate the flux-weighted mean opacity in both the Sobolev and non-Sobolev limits (Owocki & Puls 1996). In this paper, we use the Sobolev approximation to evaluate the flux-weighted mean opacity from all contributing lines and then compute the corresponding values of the statistical parameters. Applying this procedure to a range of assumed radiation and gas states, we can then tabulate the resulting parameters as a function of the input state.

2.1 Flux-Weighted Mean Opacity from a Spectral Line in the Sobolev Approximation

Assuming isotropic emission, the radiation acceleration on a medium illuminated by a frequency-dependent specific intensity l_v is

$$g_{\text{r}}=\frac{1}{\rho c}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \oint \text{d}v\text{d}\Omega k_v\widehat{n}{I}_v,}}$$(1)

where $\widehat{n}$ is the unit vector in the direction of radiation propagation, ρ is the mass density in g cm⁻³, and c is the speed of the light in km s⁻¹. The integration is performed for all frequencies dv and over solid angle dΩ = −dμ dϕ, where μ = cos θ. Furthermore, k_v is the extinction coefficient (units of inverse length), related to the mass absorption coefficient κ_v through κ_v = ρκ_v For various numerical applications, the radiation acceleration is often formulated using a flux-weighted mean opacity ${\overline{\kappa }}_{\text{F}}$ and radiation flux F with magnitude F = |F| Then g_r simply reads:

$$g_{\text{r}}=\frac{{\overline{\kappa }}_{\text{F}}·F}{c},$$(2)

with:

$$F\equiv {\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \oint \text{d}v\text{d}\Omega \widehat{n}{I}_v}.}$$(3)

A simple comparison of the equations above reveals that the flux-weighted mean opacity is a rank two diagonal tensor:

$${\overline{\kappa }}_{\text{F}}=\frac{1}{\rho }\frac{F}{F^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \oint \text{d}v\text{d}\Omega k_v\widehat{n}{I}_v}}.$$(4)

For a single spectral line with line centre frequency v₀ and gf values corresponding to a specific l − u transition between two atomic levels l (‘lower’) and u (‘upper’), the extinction coefficient reads:

$$k_v={\sigma }_{\text{cl}}\left(gf\right)\left(\frac{n_1}{g_1}-\frac{n_{\text{u}}}{g_{\text{u}}}\right)\varphi \left(v-v_0\right)=k_{\text{L}}{\varphi }_v,$$(5)

where n and g are the atomic occupation number density and statistical weight of the respective levels. Here, σ_cl = πe²/(m_ec) is the classical frequency-integrated line cross-section, with e and m_e being the charge and mass of the electron, and ϕ_v = ϕ(v − v₀) is the normalised line profile function, with v₀ the line-centre frequency; k_L then represents the frequency-integrated line-extinction coefficient. For such a single line, the components of Eq. (4) are

$${\kappa }_{\text{F}\left(i,i\right)}=\frac{k_{\text{L}}}{\rho }\frac{F_i}{F^2}{\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \oint \text{d}v\text{d}\Omega {\varphi }_v{\widehat{n}}_i{I}_v}}.$$(6)

In general, the specific intensity along some direction s can be computed from a formal solution to the radiative transfer equation:

$$I_v\left({\tau }_v\right)=I_v^{\text{core}}{\text{e}}^{-{\tau }_v}+{\displaystyle \underset{0}{\overset{{\tau }_v}{\int }}S\left({\tilde{\tau }}_v\right){\text{e}}^{-\left|{\tau }_v-{\tilde{\tau }}_v\right|}\text{d}{\tilde{\tau }}_v=I_v^{\text{dir}}+I_v^{\text{diff}}}$$(7)

with the line optical depth given by

$${\tau }_v={\displaystyle \underset{0}{\overset{\infty }{\int }}\text{d}s{k}_{\text{L}}{\varphi }_v}.$$(8)

We have here divided the specific intensity into a direct component $I_v^{\text{dir}}$, associated with the attenuated specific intensity $I_v^{\text{core}}$, and a diffuse component $I_v^{\text{diff}}$, related to the source function $S\left({\tilde{\tau }}_v\right)$ and describing reprocessed radiation. Using the Doppler formula, we transform the optical-depth equation above from physical to frequency space. Then, under the Sobolev approximation (Sobolev 1960), where the line-of-sight velocity gradient dv_s/ds and the line extinction k_L are assumed to be constant over the line resonance region, one finds

$${\tau }_v=\frac{k_{\text{L}}c}{v_0}{\left|\frac{\text{d}{v}_{\text{s}}}{\text{d}s}\right|}^{-1}{\displaystyle \underset{0}{\overset{\infty }{\int }}\text{d}v{\varphi }_v}={\tau }_{\text{S}}{\displaystyle \underset{0}{\overset{\infty }{\int }}\text{d}v{\varphi }_v},$$(9)

where τ_S is the directional-dependent Sobolev optical depth

$${\tau }_{\text{S}}\left(\widehat{n}\right)\frac{k_{\text{L}}c}{v_0}{\left|\frac{\text{d}{v}_{\text{s}}}{\text{d}s}\right|}^{-1}=\frac{k_{\text{L}}c}{v_0\left|\widehat{n}·\nabla \left(\widehat{n}·v\right)\right|}.$$(10)

For notational simplicity, here we introduce dimensionless parameters: line strength q (Gayley 1995), an illumination function w_v, and a characteristic scale τ_t for the Sobolev optical depth:

$$q=\frac{k_{\text{L}}}{v_0\rho {\kappa }_0},$$(11)

$$w_v=\frac{\pi v_0{I}_v^{\text{core}}}{F},$$(12)

$${\tau }_{\text{t}}\left(\widehat{n}\right)\frac{c\rho {\kappa }_0}{v_0\left|\widehat{n}·\nabla \left(\widehat{n}·v\right)\right|},$$(13)

where κ₀ = 0.34 cm² g⁻¹ is a nominal mass absorption coefficient set to a constant for convenience of the normalisation², and using this notation, τ_S = qτ_t. Assuming radially streaming photons from $I_v^{\text{core}}$ we directly take the integrals in Eq. (6), where the diffuse component vanishes due to isotropicity. The radial component of the flux-weighted mean opacity then gives the so-called force multiplier for the single spectral line:

$$M_{\text{line}}\left({\tau }_{\text{t}}\right)=\frac{{\kappa }_{\text{F}\left(r,r\right)}}{{\kappa }_0}=w_{v}q\frac{1-{\text{e}}^{-q{\tau }_{\text{t}}}}{q{\tau }_{\text{t}}}.$$(14)

We note that here the force multiplier is defined for the flux-weighted opacity, rather than the standard definition using the radiation force (see e.g. CAK, LC).

2.2 Force Multiplier for a List of Discrete Lines

In reality, there can be many thousands of spectral lines contributing towards the full radiation force under various conditions. Lists of spectral lines can be found in appropriate databases. For example, the Munich database we use in this work contains about 4.2 × 10⁶ lines (for more information, see Pauldrach et al. 1998, 2001). For a discrete list of intrinsically non-overlapping lines, we can evaluate the total force multiplier as a sum over all contributing lines (see Gayley 1995, LC):

$$M\left({\tau }_{\text{t}}\right)={\displaystyle \sum _{i=1}^{\text{all\hspace{0.17em}lines}}{w}_{v_i}{q}_i\frac{1-{\text{e}}^{-q_i{\tau }_{\text{t}}}}{q_i{\tau }_{\text{t}}}.}$$(15)

Here the subscript i indexes each line in our line list; to evaluate this sum for a given τ_t, we thus need to compute w_vi, and q_i for each line in our dataset.

We calculate w_vi and q_i using the same approach as LC, namely assuming LTE to equate the radiation and gas temperatures (i.e. T_gas = T_rad = T). The incoming intensity and the total flux are further set by the local gas temperature via the black-body intensity distribution:

$$I_{v_i}^{\text{core}}=\frac{2h{v}_i^3}{c^2}\frac{1}{{\text{e}}^{\raisebox{1ex}{$h{v}_i$}\!\left/ \!\raisebox{-1ex}{$k_{B}T$}\right.}-1},$$(16)

$$F={\sigma }_{\text{SB}}{T}^4,$$(17)

where h is the Planck constant, k_B the Boltzmann constant, T the gas temperature, and σ_SB the Stefan–Boltzmann constant. Combining specific intensity and flux from above in Eq. (12), we get the illumination function for an unattenuated Planck radiation (i.e., without any spectral features):

$$w_{v_i}=\frac{2\pi h}{{\sigma }_{\text{SB}}{c}^2}\frac{v_i^4}{T^4}\frac{1}{{\text{e}}^{\raisebox{1ex}{$h{v}_i$}\!\left/ \!\raisebox{-1ex}{$k_{B}T$}\right.}-1},$$(18)

where v_i is the line-centre frequency of line i.

Under this LTE assumption Eq. (5) simplifies to

$$k_{\text{L}}={\sigma }_{\text{cl}}{f}_{1\text{u}}{n}_1\left(1-{\text{e}}^{\raisebox{1ex}{$h{v}_i$}\!\left/ \!\raisebox{-1ex}{$k_{B}T$}\right.}\right),$$(19)

where f_lu= (gf)/g₁ is the oscillator strength for a given transition. Computing n₁ using the Saha and Boltzmann relations (see Appendix A for details), we also compute q_i using Eq. (11).

Setting T and ρ, we can now evaluate the sum in Eq. (15) for a given line list and an assumed range of τ_t. For T = 31 kK and ρ = 10⁻¹² g cm⁻³, we show the resulting force multiplier for the range log₁₀ τ_t = [−15, 5], as a triangle-marked green line in Fig. 1a. We note that τt depends on the ratio of density and velocity gradient rather than on each parameter individually. For a fixed density a range in τ_t is then equivalent to having a range in velocity gradient. Therefore, τ_t can be treated as an independent variable for the purposes of fitting the resulting M(τ_t) (see also LC). We can directly identify a similar behaviour of the force multiplier as LC, characterised by two distinct regimes roughly corresponding to (i) an optically thick (τ_t → ∞) repion and (ii) an optically thin (τ_t → 0)) region. In the optically thick region, log₁₀ M(τ_t) follows a nearly constant slope as a function of log₁₀ τ_t, which is a well-known result first described by CAK. Following their original notation, the force multiplier there scales with t as

$$M\left({\tau }_{\text{t}}\right)\propto {\tau }_{\text{t}}^{-\alpha },$$(20)

where α is a power-law exponent and fit parameter. Upon inspection of Fig. 1a, it is clear that such a power-law fit, shown as a blue dashed line, describes the optically thick part of the green curve well. However, as we approach the optically thin part, the force multiplier saturates to a constant value (shown in the Fig. 1a with a dashed red line) instead of continuing to increase following Eq. (20).

We can evaluate the saturation value by taking the optically thin limit τ_t → 0 of Eq. (15). Then following Gayley (1995), we define

$$\overline{Q}={\displaystyle \sum _{i=1}^{\text{all\hspace{0.17em}\hspace{0.17em}lines}}{w}_{v_i}{q}_i}.$$(21)

Using $\overline{Q}$, we can now quite accurately describe the force multiplier also in thee optically thin parts (neglected by CAK).

However, $\overline{Q}$ and α alone are still insufficient to constrain the force multiplier fully. For this, we also need to know where the optically thin-to-thick transition occurs, that is, τ_t = τ_{t trans}. To gain insight into what sets this transition, we re-examine Eqs. (15) and 21, noting that saturation can only take place for the majority of lines when

$$\frac{1-{\text{e}}^{-q_i{\tau }_{\text{t}}}}{q_i{\tau }_{\text{t}}}\lesssim 1.$$(22)

Upon closer inspection of Fig. 1a, we see that the saturation occurs already at some τ_{t trans} > 0. In other words, the line strength is limited to some effective maximum, which we define as Q₀ (see also Owocki et al. 1988). We can roughly estimate this maximum strength by locating the intersection of the power law at. high optical depth with the constant saturation level at low optical depth, such that Q₀ ≈ 1/τ_{t trans}; in Fig. 1a we highlight the intersection with a dashed black line.

Fig. 1 Force multiplier as computed by evaluating the numerical sum over the discrete set of spectral lines in the triangle-marked green curve. Panel a shows the saturation level indicated with the dashed red curve, a power-law region with the dashed blue line, and the locatton of a transition eotnl witt the dashed black line. Panel b shows the actual fit of two numerical results, where the triangle-marked green line has the same temperature and density as in (a) and the square-marked blue line assumes a different set of temperature and density. The derived fit parameters for each case are indicated in the corresponding legend.

2.3 Force Multiplier Fitted to an Exponentially Truncated Power-Law Line Distribution

Wee can now calculate the force multiplier tor various different combinattom of temperature and denstty. For each combination of (T, ρ) appropriate values of α, $\overline{Q}$, and Q₀ can be tabulated and later used to reproduce the value op the force multiplier without using the speceral line databases. This approach is quite analogous to LC, differing mainly by our choices of fitting parameters to the force multiplier However, a key advantage with the parameetrisation here it that it allows us to capture (and analyse) the underlying statistical distribution erf spectral lenes. As mentioned in Sect. 1, applications like LDI require this spectral line-distribution function to be specified in order to compute the corresponding non-Sobolev radiation line force (Owocki & Puls 1996), which as we discuss in Sect. 4.1 will now be possible using our tabulations of line strength parameters tabulated here.

To provide insights into the statistical properties of the spectral line distribution, it is instructive tostart by assuming an underlying distribution for the number of lines across the frequency and line-strength domains and then to investigate the resulting force multiplier. For the functional form of the line distribution, we use a parametrisation suggested by Gayley (1995), where the number of lines dN in some interval of strength (q, q + dq) and frequency interval (v, v + dv) is

$$\text{d}{N}_v=f\left(v\right)\frac{\overline{Q}}{\text{Γ}\left(\alpha \right)Q_0^2}{\left(\frac{q}{Q_0}\right)}^{\alpha -2}{\text{e}}^{-q/Q_0}\text{d}v\text{\hspace{0.17em}d}q.$$(23)

Here, Г(α) is the Gamma function, and f (v) describes the frequency distribution of line strength that we assume to follow f(v) = 1/v (see also Friend & Castor 1983; Puls et al. 2000). Combining this with Eq. (14), we integrate over the total number of lines in frequency range and line strength:

$$M_{\text{G}}\left({\tau }_{\text{t}}\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}{\displaystyle \int \text{d}{N}_v{w}_{v}q\frac{1-{\text{e}}^{-q{\tau }_{\text{t}}}}{q{\tau }_{\text{t}}}=\frac{\overline{Q}}{1-\alpha }\frac{{\left(1+Q_0{\tau }_{\text{t}}\right)}^{1-\alpha }-1}{Q_0{\tau }_{\text{t}}}},}$$(24)

where the subscript G refers to Gayley. This description has several advantages as compared to alternatively proposed parametrisations (see also Gayley 1995). First, it leads to a simple analytic form of the force multiplier that is sufficient to reproduce the observed asymptotic behaviour found for a discrete set of spectral lines. We can verify this by simply taking the optically thick and thin limits in Eq. (24). In the optically thick limit, when Q₀τ_t ≫ 1, we directly obtain the same power-law dependence as in Eq. (20)

$$M_{\text{G}}\left({\tau }_{\text{t}}\right)\approx \frac{\overline{Q}{Q}_0^{-\alpha }}{1-\alpha }{\tau }_{\text{t}}^{-\alpha },$$

such that when using the Gayley (1995) ansatz $Q_0=\overline{Q}$, we recover the well-known form of the CAK force multiplier:

$$M_{\text{G}}\left({\tau }_{\text{t}}\right)=\frac{{\overline{Q}}^{1-\alpha }}{1-\alpha }{\tau }_{\text{t}}^{-\alpha }.$$

In the opposite limit Q₀τ_t ≪ 1, a first-order Taylor expansion of Eq. (24) gives:

$$M_{\text{G}}\left({\tau }_{\text{t}}\right)\approx \overline{Q}.$$

This means that, for each combination of (T, ρ) used to compute $w_{v_i}$ and q_i for all lines in our database, Eq. (24) can be applied to fit the M(τ_t) that results from Eq. (15) across a very wide range of τ_t, covering both the optically thick and thin limits and using only the three basic fit parameters $\overline{Q}$, Q₀, and α. Second, this parametrisation offers insight into the underlying distribution of spectral lines through intuitive parameters. As in Sect. 2.2, $\overline{Q}$ and Q₀ correspond to the maximum force in the optically thin limit and the effective maximum strength of spectral lines, respectively. The values of the α, $\overline{Q}$, and Q₀ parameters can now be found by directly fitting Eq. (24) to the actual force multiplier as computed from the discrete line list. We note that here we do not consider the so-called δ parameter, introduced by Abbott (1982) to account for the shift of ionisation states over the wind. In cases when the line-distribution parameters are held fixed throughout the full wind, the δ parameter then accounts for the fact they should vary with position. This is typically applied in studies of 1D line-driven stellar winds, where the linedistribution parameters have often been tabulated as a function of stellar parameters (i.e. effective temperature, mass and radius of the star, Abbott 1982), and each considered wind model thus has been characterised by one single value for the line-distribution parameters. Those models then require the δ parameter to introduce back the radial variation neglected by the use of spatially fixed line-force parameters (see also Gormaz-Matamala et al. 2019). In our new method presented here, these effects are instead directly accounted for by allowing the line-distribution parameters to vary with the local density and temperature (thus, they also vary naturally with both position and ionisation; see a first benchmark application in Sect. 3).

Figure 1a shows the results of such a fit for two different combinations of temperature and density. Here, we show the same force multiplier as in Fig. 1a as well as the force multiplier for T = 10 kK and ρ = 10⁻¹⁷ g cm⁻³, along the corresponding best-fitting M_G(τ_t). Though we note that this parametrisation might not be sufficient to fully capture some of the fine details of the optically thick region, as well as perfectly reproduce the force multiplier in the transition region (see the discussion in LC), it is sufficient to represent the overall character of the force multiplier while also providing insight on the spectral line distribution. As shown explicitly in Sect. 3.3, the range in τ_t over which we compute M(τ_t) is sufficient to cover the full range of values present in our simple 1D RHD benchmark simulations of O-star winds. Moreover, since for a given (ρ, T) pair M(τ_t) saturates at τ_t ≪ τ_{t trans} and follows a single power law for τ_t ≫ τ_{t trans}, fitting over an even larger range of τ_t would in any case not yield meaningfully different values for the fit parameters. The specific values for density and temperature showcased in the triangle-marked green curve in Fig. 1a are rather good representations of the expected values for the already quite extensively analysed regime of early O-supergiants. Our best fit values α = 0.65, $\overline{Q}=Q_0=1104$ are indeed broadly consistent with previous studies of this region (Gayley 1995; Puls et al. 2000). In contrast, for the second pair of density and temperature, which is outside the early O-supergiant regime, the best fit parameters α = 0.48, $\overline{Q}=35$ and Q₀ = 2634 are very much outside these previously estimated ranges. This illustrates a general need for not assuming constant line-force parameters independent of conditions, which so far has often been the case in (at least time-dependent) simulations of line-driven flows. We will thus now expand the above calculations towards a much broader range of temperature and density combinations.

2.4 Grid of Line-Distribution Parameters

We apply the method described above to various combinations of temperature and density in the range log₁₀(ρ [g⁻¹cm³]) ∈ [−20, −10] and log₁₀(T [K⁻¹]) ∊ [4.0, 4.7], assuming a solar plasma composition as given by Asplund et al. (2009). Using colour maps, we present the results of log₁₀ $\overline{Q}$, ${\log }_{10}\left(Q_0\text{/}\overline{Q}\right)$, and α in Fig. 2.

Examining these figures, we focus on the region of interest for our initial study of O-type stars in Sect. 3. Typical surface conditions in O stars cover the region of log₁₀(ρ [g⁻¹ cm³]) ∈ [−14, −10] and log₁₀(T [K⁻¹]) ∈ [4.3, 4.7]. Inspecting this region in Figs. 2a, b, we find values of α ≈ 0.5–0.7 and $\overline{Q}\approx 1000–2000$ that again are in good agreement with previous estimates (e.g. Gayley 1995; Puls et al. 2000). Moreover, for the O-star regime, we also find that the ansatz $Q_0\approx \overline{Q}$ by Gayley (1995) holds reasonably well. However, we see that the tabulated values can largely deviate from these often used estimates outside this region. For example, also in the O-star temperature regime at moderately low densities (typical for the outer wind), we find $\overline{Q}<1000$ and α > 0.7. As further discussed in Sect. 3, such variation in the line-force parameters may then significantly affect the radiation force and alter the wind dynamics.

2.5 From Tables to Opacity

The tables presented here are now ready to be applied towards computations of the flux-weighted mean opacity, which can then be used in a broad range of applications where the radiation line force needs to be accounted for. In the Sobolev limit, we find a general expression for the flux-weighted opacity diagonal tensor by combining Eqs. (6) and (24):

$${\overline{\kappa }}_{\text{F}}=\frac{F}{F^2}\frac{{\kappa }_0\overline{Q}}{1-\alpha }{\displaystyle \oint \text{dΩ}\widehat{n}{I}^{\text{core}}\left(\widehat{n}\right)\frac{{\left[1+Q_0{\tau }_{\text{t}}\left(\widehat{n}\right)\right]}^{1-\alpha }-1}{Q_0{\tau }_{\text{t}}\left(\widehat{n}\right)}}.$$(25)

We note that in this general expression, we have allowed $I^{\text{core}}\left(\widehat{n}\right)={\displaystyle \int \text{d}v{I}_v^{\text{core}}\left(\widehat{n}\right)}$ and ${\tau }_{\text{t}}\left(\widehat{n}\right)$ to potentially be functions of the ray direction, whereas the parameters $\overline{Q}$, and Q₀ are all local properties of the gas, which for convenience have all been tabulated using the radial case (Sect. 2.1). As also noted by LC, since the different ray directions effectively serve to provide a scaling factor to the radial streaming case, calibrating our tables using the radial ray is justified.

In our numerical models, for each separate spatial cell at each time-step, we find appropriate values of α, $\overline{Q}$, Q₀, and τ_t based on the local density, temperature, and line-of-sight velocity gradient. In practice, values of α, $\overline{Q}$, and Q₀ are found by interpolations in our pre-described tables (i.e., very similar to how tables of Rosseland mean opacities often are used in various applications). In this first study, we further restrict ourselves to analysing the flux-weighted line opacity in the Sobolev limit.

However, in Sect. 4.2.2 we discuss a possible extension also towards non-Sobolev (observer’s frame) applications.

Fig. 2 Colour maps of the line-distribution statistical parameters (a) ${\log }_{10}\overline{Q}$, (b) α, and (c) ${\log }_{10}\left(Q_0\text{/}\overline{Q}\right)$ computed for a range of temperatures and densities using the Munich line database. The abscissa here measures the logarithm of the temperature in kK, and the ordinate gives the logarithm of the density in g cm−3.

3 First Applications

Generally, self-consistent massive-star wind calculations using detailed radiation transfer codes such as FASTWIND (Sundqvist et al. 2019). METUJE (Krtička & Kubát 2017), PoWR (Sander et al. 2017), or CMFGEN (Gormaz-Matamala et al. 2021) are highly time-consuming. Those computations are also exclusively 1D and only consider spherically symmetric and steady-state wind outflows. Compared to these steady-state models, relaxing a 1D O-star model by means of the approach introduced here increases the computational speed by a significant factor, > 100. Thus, our new method offers a direct way to straightforwardly take the full line statistics into account also within time-dependent 2D and 3D simulations (something that is computationally prohibited when using the time-consuming radiative transfer methods applied in the codes mentioned above). Furthermore, in comparison to previous time-dependent line-driven wind models (see e.g. Feldmeier & Thomas 2017; Dyda & Proga 2018), which typically rely on line-force parameters that are pre-described and fixed throughout the simulations (i.e. both in space and time), our new models allow these parameters to vary in space and time, self-consistently computing and updating them according to the local wind conditions.

In this initial study, our key objective is to benchmark our method and tables. To achieve this, we investigate the wind outflows of a calibrated O-star model grid used by Björklund et al. (2021). The values of stellar luminosity ℒ, stellar mass ℳ, and stellar radius ℛ of the considered models are summarised in Table 1. In addition to studying the behaviour of the line-force parameters, the O-star grid allows us to further validate our technique by comparing mass-loss rates found in our simulations to the values predicted for such stars by other studies using more detailed radiative transfer.

3.1 Numerical Simulations of Line-Driven O-star Winds

The radial variation in density and velocity of a spherically symmetric outflow are described by the hydrodynamical equations of mass and momentum conservation:

$$\frac{\partial \rho }{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\rho v\right)=0,$$(26)

$$\frac{\partial }{\partial t}\left(\rho v\right)+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\rho v\right)=-\frac{\partial P}{\partial r}-\rho g+\rho \frac{{\kappa }_{\text{e}}{F}_{\text{r}}}{c}+\rho \frac{{\kappa }_{\text{F}\left(r,r\right)}{F}_{\text{r}}}{c}.$$(27)

Here, υ is the radial velocity in cm s⁻¹ , g = ℳG/r² is the gravitational acceleration in cm s⁻², and P is the gas pressure in g cm⁻¹ s⁻². Using the ideal gas law, P = a²ρ, where a² = k_BT/(μm_H) is the isothermal sound speed, μ = 0.66 is the mean molecular weight of a fully ionised solar plasma, and m_H is the mass of the hydrogen atom. We further assume that the winds are isothermal, using the standard assumption (for O stars) T = 0.8T_eff (Puls et al. 2000), where $T_{\text{eff}}^4=ℒ\text{/}\left(4\pi {\sigma }_{\text{SB}}{R}^2\right)$ is the effective temperature of the star³. Finally, κ_eF_r/c and κ_Fr,rF_r/c are the two components of radiation acceleration from Thompson scattering and spectral lines, respectively, with F_r = ℒ/(4πr²) the radial radiation flux.

Assuming a uniformly bright star with I^core = F/π for μ ∈ [μ_*, 1], where ${\mu }_{*}^2=1-R^2\text{/}{r}^2$, Eq. (25) gives:

$${\kappa }_{\text{F}\left(r,r\right)}={\kappa }_0\eta \frac{\overline{Q}}{1-\alpha }\frac{{\left(1+Q_0{\tau }_{\text{t}}\right)}^{1-\alpha }-1}{Q_0{\tau }_{\text{t}}}{\text{cm}}^2{\text{g}}^{-1},$$(28)

where the Sobolev optical depth parameter from Eq. (13) reduces to τ_t = c ρκ₀/|dυ/dr|. Here η is a geometrical correction factor that accounts for the finite size of the star, the exact computation of which requires numerical evaluation of the solid angle integral given in Eq. (25). However, assuming Q₀τ_t ≫ 1 when evaluating η provides an analytic estimate (Pauldrach et al. 1986; see also Sect. 2.2.4 of Cranmer 1996):

$$\eta \left(r\right)=\frac{{\left(1+\sigma \right)}^{1+\alpha }-{\left(1+\sigma {\mu }_{*}^2\right)}^{1+\alpha }}{\left(1+\alpha \right)\left(1-{\mu }_{*}^2\right){\left(1+\sigma \right)}^{\alpha }\sigma },$$(29)

with homologous wind expansion parameter

$$\sigma =\frac{\text{d}\ln v}{\text{d}\ln r}-1.$$

For simplicity, we introduce the effective gravity g_eff(r) = ℳG[1 − Γ_e(r)]/r², where Γ_e(r) = κ_e(r)ℒ/(4πℳGc) is the Eddington ratio, so that Eq. (27) reads as

$$\frac{\partial }{\partial t}\left(\rho v\right)+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\rho v\right)=-\frac{\partial P}{\partial r}-\rho g_{\text{eff}}+\rho \frac{{\kappa }_{\text{F}\left(r,r\right)}{F}_{\text{r}}}{c}.$$(30)

To compute the effective gravity, we use the appropriate pre-tabulated values of the Thompson scattering opacity κ_e.

We used the midpoint time integration, with TVDLF Riemann solver, and a VANLEER flux limiter (van Leer 1979). As described in the previous section, computation of κ_F(r,r) requires knowledge of the local velocity gradient. We computed this by applying a simple finite differencing for non-constant intervals to the velocity field at each numerical time step (Sundqvist & Veronis 1970). We also limited the integration time step dt to avoid information propagating over multiple cells in a single step by

$$\text{d}t=\min \left(\text{d}{t}_{\text{CFL}},\text{CFL}\times \sqrt{\frac{c\text{Δ}r}{{\kappa }_{\text{F}\left(r,r\right)}{F}_{\text{r}}}}\right),$$

where dt_CFL is the time step set by the standard Courant–Friedrichs–Lewy (CFL) condition (Courant et al. 1928), using here CFL = 0.3, and ∆r is the width of one numerical cell.

Each simulation covers a radial range extending from a lower boundary at ℛ to an outer boundary at 6 ℛ. For the restricted parameter range investigated in this paper (essentially O stars), such an outer boundary is sufficient to ensure that the winds have reached velocities sufficiently close to their asymptotic values (i.e. they are not significantly accelerating at the outer boundary). In our numerical setup, we further apply a constant stretching factor ∆r_í+1/∆r_í = 1.002 between subsequent cells i, i + 1 to resolve the transonic region.

Fig. 3 Colour map showing the radial and temporal variation in the line-distribution parameters. From the left to right we show α, $\overline{Q}$, and Q0 radially normalised to their respective final profiles αfin, ${\overline{Q}}_{\text{fin}}$, and Q0 fin.

3.2 Initial and Boundary Conditions

We set the initial velocity using the so-called β-law:

$$v_{\beta }\left(r\right)=v_{\infty }{\left(1-b\frac{R}{r}\right)}^{\beta },$$(31)

with β = 0.5 and terminal velocity $v_{\infty }^2=2R{g}_{\text{eff}}\left(R\right)$. To avoid infinite density in the first computational cell, we chose b such that υ_β(ℛ) = 10 km s⁻¹. The initial density follows

$$\rho \left(r\right)={\rho }_{\text{b}}\frac{v_{\beta }\left(R\right)}{v_{\beta }\left(r\right)}{\left(\frac{R}{r}\right)}^2,$$(32)

where ρ_b is an input parameter. In principle, we could use an analytic solution of mass loss (see CAK for more details) for the initial condition, rather than setting the input value ρ_b. However, we chose not to do so to investigate the relaxation of the time-dependent solutions starting from different initial struc tures. Indeed, starting the same setup with different values of ρ_b (varying it by an order of magnitude), we find almost identical final structures, thus confirming that our results are quite independent of the specific initial conditions.

We set the upper boundary conditions by simple extrapolations of velocity and density (see e.g. Poniatowski et al. 2021). In contrast to our previous simulations, we here allow the lower boundary conditions to self-adjust) to the overlaying wind by updating the density of the first ghost сell ρ₁ when required. This is done if the following, condetion is met: $\left|{\rho }_1\text{/}\left(2{\overline{\rho }}_{\text{s}}\right)-1\right|\ge 20\%$, where ${\overline{\rho }}_{\text{s}}$ is the average density at the sonic point υ(R_s)² = a², with the averages taken over a dynamical timescale t_dyn = ℛ/υ_∞. When ever this condition is met, we set ${\rho }_1=2{\overline{\rho }}_{\text{s}}$. The remaining ghost celte at the lower boundary are then set by assuming quasi-hydrostatic equilibrium for the density, where after the velocity is determined from most conservation. In comparison to our previous models with ρ₁ fixed throughout the simulation duration, this procedure allows the wind to adjust itself somewhat better to initial conditions that may be rather far away from the final relaxed state.

3.3 Mod-0

Following the above setup, we first present a numerical model (Mod-0) with stellar parameters roughly corresponding to ζ Puppis, which is often referred to as a ‘prototypical’ O star.

The distinctive point of our simulation, as mentioned above, is that here the line-distribution parameters are interpulated⁴ from tables and thus consistent with the local conditions in the wind; therefore, they vary, in time and space. That is, as outlined in Secl. 2, we are computing α(ρ, T), $\overline{Q}\left(\rho ,T\right)$, Q₀(ρ, T) at each radius and each time-step in the simulation, and then we use Eq. (28) to obrain the corresponding κ_F(r, t) that is consistent with the local hydrodynamical structure. Since the local values of α, $\overline{Q}$, and Q₀ set the scale of the line force, it is natural that the radial and temporal variation in these parameter then also reflects upon the self-consistent mass loss and velocities of the models. To illustrate this variation, (in Fig. 3) we plot colour maps of all three parameters as a function of scaled radius x = 1 − ℛ/r and simulation time t in unite of dynamical timescale t_dyn. For better visualisation, we normalised each parameter to their radial profile from the final snapshot. As seen in the figure, after adjusting to the new conditions all three parameters relax to a quite steady behaviour in time. To ensure that our models are well retaxed, we typically run models for at least 90 t_dyn, although in the О star regime, the relaxation Is typically faster and a (quasi) steady-state ‘final’ configuration is reached much earlier than this. Because the relaxed radial profiles of α, $\overline{Q}$, and Q₀ are significantly different from the typically assumed constant profiles, we now discuss how this variation influences the wind dynamics.

We show the radial profile of α in Fig. 4a. Near the lower boundary, we observe that α starts with values reasonably close to the typically expected 2/3 (see e.g. Puls et al. 2008); however, when the radius increases, so does α. A possible interpretation of this rise comes from the distribution of optically thick and thin lines that contribute to the line force. Splitting Eq. (15) into two parts,

$$M\left({\tau }_{\text{t}}\right)\approx {\displaystyle \sum _i^{q_i{\tau }_{\text{t}}\ll 1}{w}_{v_i}{q}_i}+{\displaystyle \sum _j^{q_j{\tau }_{\text{t}}\gg 1}\frac{w_{v_j}}{{\tau }_{\text{t}}},}$$(33)

we notice how the optically thin part has no direct dependence on the velocity gradient (α = 0), whereas the optically thick part is directly proportional to it (α = 1). Following Puls et al. (2000), we may, thus interpret the value of α as the fraction of the line force that comes from optically thick lines to the totol line force. As such, we can associate the increase in α with radius as a change of the line force character from having relatively even contributions from weak and strong lines to being dominated by more optically thick lines. Indeed, driving the outer wind of О stars by only a few strong lines is consistent with previous models that focus exclusively on stationary winds (for example, see Fig. 9 in Krtička 2006 and the general discussion in Sect. 2 of Puls et al. 2008⁵). Moreover, this general interpretation seems so be supported by observations of O stars, where we see only a few strong spectral fines formed at a high wind velocity. Also, since differences in abundances or ionisation stages, or both can significantly influence these relatively few outer-wind driving lines, a significant scatter of terminal wind speeds can be expected (see also Puls et al. 2000), as is also observed (Garcia et al. 2014).

Inspecting Fig. 4b we further find a decrease in both $\overline{Q}$ and Q₀ towards the outer wind. Generally, the decrease of $\overline{Q}$ suggests either that the contributing lines in the outer wind are weaker compared to the inner wind or that fewer lines are contributing to the total line force. Both an increase in α and a decrease in $\overline{Q}$ lead to a lower radiation force, which in turn also reduce the wind velocity. However, a lower Q₀ can counteract; this, as this leads to the opposite effect, causing an increase in the radiation force (by effectively reducing me line optical depth).

In principle, as each ok the parameters depends on the local wind conditions while also affecting these locar conditions, they form feedback loops with the wind structure. To understand this, we consider the effects of varying each parameter separately in Eq. (28). Then an increase in α leads to a decrease in the line opacity, which affects the velocity and density profile. However, in the O-star wind regime, an increase in density reduces α (see Fig. 21b), thus forming positive feedback. Similarly, from Fig. 2a we see that increasing density also increases $\overline{Q}$, and thus the line opacity, again forming a positive feedback loop. By contrast, Q₀ and the line force form a negative feedback loop, where increasing Q₀ rather decreases the line opacity. A self-regulatory system forms through the interaction of these positive and negative feedback loops, which ultimately here determines the final mass-loss rate and velocity law of the wind.

In Fig. 5a, we show the relaxation of the optical depth parameter τ_t. Similar to the relaxation of the line-distribution parameters (Fig. 3), we observe a transient behaviour from the initial condition, which takes about 40 t_dyn to fully vanish. After this, a near-equilibrium configuration is reached. Importantly, the value of τ_t during the simulation never exceeds the range initially used to compute the line-distribution parameters. As such, our tabulations always cover the correct range in τ_t (as well as in density, velocity gradient and temperature)⁶. Inspecting the final, relaxed values of τ_t in Fig. 5b, we further note that the final saturation of M(τ_t) at τ_t ≪ τ_{t trans} is not reached for this particular simulation. Indeed, the minimum values still have Q₀ τ_t ≫ 1, and as such the full O-star wind resides in the power-law portion of M(τ_t), as also previously noted, for instance by Gormaz-Matamala et al. (2019). However, we emphasise that this is not guaranteed always to be the case, especially not for more complicated situations where the outflow might not be stationary and smooth, but rather time-dependent and clumpy with co-existing regions of (much) higher and lower densities. In such cases, very low values of τ_t would systematically be found in the low-density parts, and M(τ_t) might there readily reach the regime where it saturates to a maximum value $\overline{Q}$.

The resulting velocity structure can be seen in Fig. 6a (solid blue line), where we show the velocity in km s⁻¹ as a function of x, again for the final snapshot. For comparison, we show the initial velocity (dashed golden line) prescribed by Eq. (31). As is evident from the comparison, the final velocity profile does not resemble the initially assumed one very much, exhibiting a different acceleration profile and a higher terminal velocity of υ_∞ ≈ 1600 km s⁻¹. An increase in terminal velocity between the initial (CAK) and the final profile is typically expected because of the finite disc factor ç (Pauldrach et al. 1986; Friend & Abbott 1986). However, variation in α, $\overline{Q}$, and Q₀ will also affect the terminal velocity. To differentiate between effects introduced by the finite disc factor from the effects introduced by variation in line-distribution parameters, we also compute a ‘reference’ model, which uses the same stellar parameters as the Mod-0, but where we keep line-distribution parameters constant throughout the simulation. In this reference model, we use the near lower boundary conditions of the original model to set the values α = 0.66, $Q_0=\overline{Q}=1100$. This then results in a terminal velocity of $v_{\infty }^{\text{ref}}\approx 3000\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}$. Clearly, compared to the reference model, the terminal velocity is reduced by a factor of ~2. Since the only difference between the two models is how we treat the line-distribution parameters, this stems directly from their specific radial behaviour.

Further, to demonstrate that the wind velocity and density also relax as line-distribution parameters reach their final configuration, in Fig. 6b we plot a colour map of radial and temporal variation in the mass-loss rate zooming in on the first 30 t_dyn. The axes in this figure are the same as in Fig. 3. Like for the line-force parameters the relaxation takes around 20 t_dyn after which the model relaxes to a final self-consistent mass loss of Ṁ = 1.8 × 10⁻⁶M_⊙yr⁻¹. This mass loss then stays constant and only differs by about 20% from the value predicted by Björklund et al. (2021) and Krtička & Kubát (2017) for the same luminosity (see Table 1). Such an overall good agreement of our models with previous mass-loss results re-emphasises the general quality of the method. Indeed, as now shown, our computed mass-loss rates agree quite well with the predictions by Björklund et al. (2021) and Krtička & Kubát (2017) in the complete O-star domain.

Fig. 4 Radial profile of the line-distribution parameters of the final snapshot of Mod-0. Plot (a) shows the α parameter and (b) the $\overline{Q}$ parameter (solid red line) and Q0 parameter (dashed black line).

Fig. 5 Radial and temporal variation of the scaled Sobolev optical depth τt. (a) Colour map showing the radial and temporal variation in optical depth parameter τt. Extrema found in this map are max(log10 τt) = 3.95 and min(log10 τt) = 0.12. (b) Radial profile of τt from the final snapshot.

Fig. 6 Comparison of the initial and final velocity profiles, and the radial and temporal variation of the mass-loss rate, (a) Radial velocity profiles of the initial conditions with the dashed orange line and the final snapshot with the solid blue line, (b) Colour map showing the relaxation or the radial profile of the logarithm of the mass-loss rate spanning the range from the initial conditions up to 30 tdyn.

3.4 Grid of O Stars

Using the same setup as for Mod-0, we compute numerical models for the stellar parameters presented in Table 1. In the previous section, we considered in detail a specific model, pointing out that the discussed behaviours are common for all models within our grid. Specifically,in all models, α increases outwards in the wind, whereas both $\overline{Q}$ and Q₀ decrease. Near the lower boundary, we find values of all three line-force parameters within their previously estimated range. However, starting just slightly above the stellar surface, they begin to deviate from this range by significant amounts. Across all O-star models, we find $Q_0\approx \overline{Q}$ only in the near-surface region, whereas throughout the rest of the wind $Q_0<\overline{Q}$. In general, this means that studies using the standard ansatz $Q_0=\overline{Q}$ tend to overestimate the effective Sobolev optical depth Q₀τ_t of the line ensemble in low-density regions, which can then lead to an inconsistent radial line force.

Furthermore, although the typical conviction is that non-LTE (NLTE) effects are essential in setting the dynamics of O-star winds, our model-grid analysis seems to suggest the opposite, namely that an LTE-based line force (explicit assumption in Sect. 2.2) sufficiently well captures the fundamental dynamics of (at least) the near-star region, which essentially sets the mass-loss rate of the model. This suggestion stems from the (on average) tight agreement of the mass-loss rates found in our models with the predictions derived from the full NLTE models by Krtička & Kubát (2017) and Björklund et al. (2021).

To illustrate this point, we plot in Fig. 7 the resulting mass-loss rates as a function of luminosity. In addition, we performed a linear fit plotted and compare this mass-loss-luminosity relation to the similar relations by Krtička & Kubát (2017) and Björklund et al. (2021). We indeed observe exceptionally intimate agreement between our and the Krtička & Kubát (2017) relations and also a close correlation with the Björklund et al. (2021) results, with better agreement at the high luminosity end (see also below). For comparison, we also plot the mass-loss-luminosity relation by Vink et al. (2001). As compared to Vink et al. (2001), our models have a systematically lower mass loss by about a factor of approximately 3. Such a Tactor of 3’ difference in the O-star regime is (nowadays) quite well-established (see e.g. Najarro et al. 2011; Sundqvist et al. 2011; Šurlan et al. 2013; Cohen et al. 2014; Hawcroft et al. 2021), and our results thus suggest that this reduction is not a consequence of NLTE or co-moving frame transfer (CMF) effects. Rather our models seem to suggest that it is of prime importance to require local dynamical consistency (as in the models by Krtička & Kubát 2017; Björklund et al. 2021, and here), rather than relying on a global energy constraint (as in the models by Vink et al. 2001), when computing the mass-loss rates (for О stars). We also note that for lower luminosities, our rates start to become higher than those predicted by Björklund et al. (2021), which, most likely, is a direct result of the characteristic radiation-force dip in near-sonic regions observed in non-Sobolev models (Björklund et al. 2021).

Interestingly, our results thus suggest that this effect is more significant for lower luminosities and mass-loss rates, as already pointed out by Owocki & Puls (1999). Furthermore, Krtička & Kubát (2017) scale their computed non-Sobolev line force to the corresponding Sobolev one, which shifts their critical point (where the mass loss is set) upstream of the sonic point. As a result, their models might be considered as ‘effectively Sobolev-based’ and the nature of their steady-state solutions may thus be quite similar to those presented here, consequently giving a possible explanation to why we obtain such a tight agreement with these authors over the complete O-star domain. Regarding terminal wind speeds, we note that we generally obtain somewhat lower values than those by Björklund et al. (2021). This may stem from differences in the modelling techniques (e.g. CMF or our LTE assumption; for the latter, see the discussion in Sect. 4.2), and should be further investigated in future work.

In summary, the overall agreement with these previous works in the well-studied O-star domain provides good support for the general method developed here. We remind the reader, however, that our prime goal in this section has not been to provide new velocities and mass-loss rates for О stars, but rather to test our new (fast and efficient) method against models based on more detailed radiative transfer (but which are limited to ID steady-state systems). The key point of the analysis in this section is thus that, with the above benchmarks now in hand, our method is ready for exploitation also in more complicated situations requiring multi-dimensional RHD setups (e.g. Moens et al. 2022a).

Fig. 7 Comparison of the mass-los s rates from our grid of models (blue dots) to the mas s-los s values predicted in earlier works by Björklund et al. (2021), Krtička & Kubát (2017), and Vink et al. (2001). We plot the mass-loss-luminosity relations predicted by these authors with dashed orange, dotted-dashed maroon, and dotted black lines, respectively. We also over-plot the linear fit to our results (solid blue line).

4 Discussion

4.1 Comparison to Lattimer & Cranmer (2021)

Using a similar set of assumptions, LC recently computed force multipliers from an atomic database with approximately the same number of spectral lines as ours. However, a key distinction between theirs and our work is that they fit the force multipliers using an ad hoc functional form that seems to have no apparent relation to an underlying line-distribution function (their Eq. 37). By contrast, the functional form we fit here (Eq. (24)) can be directly derived from assuming a distribution of lines according to Eq. (23). As such, our study becomes much simpler to generalise both towards general multi-dimensional RHD simulations (see the discussion in Sect. 3) and towards situations wherein a non-Sobolev line force must be considered (see the discussion in Sect. 4.2.2).

Nonetheless, as our computations of the force multiplier are very similar, we may extract the corresponding $\overline{Q}$-values from the tabulations of LC in order to compare results stemming from different atomic line databases. To achieve this, we first re-normalise the results found by LC the same way as described in Sect. 2.1 and then, in Fig. 8a, we plot the colour map of ${\log }_{10}{\overline{Q}}_{\text{LC}}$ for a similar temperature and density grid as presented here . For a quantitative comparison to our results (seen in Fig. 2a), Fig. 8b then plots the ratio ${\log }_{10}\left({\overline{Q}}_{\text{LC}}/\overline{\text{Q}}\right)$.

Inspecting these figures, we immediately see that the two independent results are broadly in agreement, with a small-factor difference throughout most of the parameter space. However, an exception is the low-density and high-temperature region, where the difference is significant. Since the only principal difference between our studies regards the line databases, this mismatch should stem from differences in these databases.

We note that the line database compiled by LC for elements includes ionisation stages up to X (i.e. nine times ionised), while the Munich line database (see Appendix B) only extends to ionisation stages VIII. This difference in the available maximum ionisation stages becomes significant at low densities and high temperatures, as in this regime, it becomes increasingly easier to ionise elements. Therefore, the significant difference between $\overline{Q}$ and ${\overline{Q}}_{\text{LC}}$ in this region is most likely explained by the lack of spectral lines from high ionisation stages in our database. Nonetheless, we here focus on O-star applications, where the local conditions do not reach this corner of parameter space and thus should not pose an issue for our benchmark study. Indeed, examining the region of Fig. 8b corresponding to the density and temperatures of such O-star winds, the maximum difference between the two calculations is less than 20% in $\overline{Q}$.

Fig. 8 Colour maps of the ${\log }_{10}{\overline{Q}}_{\text{LC}}$ computed by LC (a) and differences between the ${\log }_{10}{\overline{Q}}_{\text{LC}}$ and ${\log }_{10}\overline{Q}$ presented in Fig. 2a (b). Maps are given as a function of temperature on the abscissa and density on the ordinate. The ${\log }_{10}{\overline{Q}}_{\text{LC}}$ is reconstructed from the data provided by LC.

4.2 Limitations of the Model

Having discussed the influence of the atomic data for our results, we now focus on some of the limitations of the model. These limitations stem from the various approximations introduced in Sect. 2, and here we focus on the three most important ones.

4.2.1 NLTE Effects and Non-Planckian Illumination

In our work, we assume LTE and set the radiation and gas temperatures to be equal. Overall, as also discussed in Sect. 3.3, this approximation seems to work reasonably well, resulting in O-star mass-loss rates that agree well with predictions derived from alternative models accounting for NLTE effects. Nonetheless, NLTE effects are generally known to affect the detailed ionisation/excitation balance in hot stars, so a natural follow-up investigation to this study will be to examine effects from this also upon the detailed line-force tables presented here. As such, in a forthcoming paper, we will extend our method towards an approximate NLTE treatment (allowing the gas and radiation temperatures to be different) following, for example, the basic method outlined in Puls et al. (2000).

In order to maintain generality when tabulating line-distribution parameters, we approximated the frequency distribution of the illuminating specific intensity with the Planck function. However, in a more realistic scenario, illumination may not be strictly Planckian, which might then affect the flux-weighted mean opacity. In particular, if lines overlap within the velocity domain of an outflow, this can lead to multiple line-resonances, which may then affect the illumination of the individual resonance zones (see e.g. Puls et al. 1993). Moreover, if the gas is exposed to an external strong ionising radiation field, as for example in active galactic nuclei and high-mass X-ray binaries, this can further alter both the ionisation balance and the illuminating intensities (see e.g. Dannen et al. 2019). Finally, by design, the straight line-force sum, Eq. (15), neglects all intrinsic line-overlap effects. This may act towards reducing the line force due to self-shadowing effects within intrinsically overlapping optically thick lines.

4.2.2 Non-Sobolev Line Force

In Sect. 2.3, we mentioned how a key advantage of the line-force parametrisation used here is that, it yields an explicit form for the underlying statistical distribution of lines. Using this statistical line-distribution function (i.e. Eq. (23)), one can then also compute a corresponding non-Sobolev line force based on our tabulations, for example, by following the observer’s frame escape-integral methods outlined in Owocki & Puls (1996). This would then allow us to investigate the effects from varying the line-force parameters upon the strong LDI (Owocki & Rybicki 1984), which is not captured within a Sobolev treatment. In terms of the general method developed here, a straightforward application would be to use our line-force tabulations directly in such LDI calculations. Then LDI models could be constructed without any (significant) further loss of computation time, and the effects of varying line-force parameters can be investigated in an analogous way as in this paper. Thus far, LDI simulations have always assumed prototypical line-force parameters (constant both in time and space), although it is clear that varying them will also impact the resulting LDI-generated structure (Driessen et al. 2019). First simulations of the LDI using self-consistent line-force parameters are already underway and will be presented in a forthcoming paper.

5 Conclusions

This paper has presented a comprehensive tabulation of line-distribution parameters that can be used in the fast calculation of flux-weighted line opacities and line forces in supersonic media (Sect. 2). In Sect. 3, we demonstrate that, when applied to the (rather well-calibrated) line-driven flows of О stars, our method, considering the limitations of the method outlined in Sect. 4.2, gives a very good agreement between our derived mass-loss rates and those obtained in the detailed steady-state studies by Krtička & Kubát (2017) and Björklund et al. (2021). Based on this benchmark study, we provide here some important conclusions.

First, the self-consistent variation in the line-strength distribution parameters (i.e. α, $\overline{Q}$, and Q₀) have critical feedback on the resulting line-driven flow structure, influencing the mass loss and velocity structure in the models. Thus, we conclude that the conventional approach of keeping these parameters constant in time and space within multi-dimensional, time-dependent line-driven flow simulations (e.g. Kee et al. 2016; Dyda & Proga 2018; Sundqvist et al. 2018; Schrøder et al. 2021) may lead to inconsistent results.

Second, assuming LTE when computing the line force seems to work reasonably well, at least in the near-star regions of О stars, where their mass-loss rates are determined. Nonetheless, a key follow-up study to this paper will present an extension of the basic method presented here towards (approximate) NLTE conditions, allowing the gas and radiation temperatures to be different.

Finally, we conclude that the method presented here for computing flux-weighted line opacities is a beneficial, versatile, and very fast tool for implementing more realistic treatments of radiation line forces in multi-dimensional, time-dependent RHD studies. For example, in a unified simulation encompassing both the optically thick, deep, subsonic regions and an overlying line-driven wind, the method here can be trivially applied within the hybrid-opacity model proposed by Poniatowski et al. (2021) (wherein ê_F is simply added to Rosseland mean opacities; see also the discussion in Castor 2004, ch. 6). Indeed, with this hybrid opacity model, our new tables have already been implemented into the RHD code presented in Moens et al. (2022b), and first simulations of 3D time-dependent WR wind outflows have been computed and will be presented in a parallel paper (Moens et al. 2022a).

Acknowledgements

The authors would like to thank all members of the KUL EQUATION group for fruitful discussion, comments, and suggestions. L.P., J.S., L.D., A.dK., and H.S. acknowledge support from the KU Leuven CI grant MAESTRO С16/17/007. JS further acknowledges additional support by the Belgian Research Foundation Flanders (FWO) Odysseus program under grant number G0H9218N. NDK acknowledges support from the National Solar Observatory (NSO), which is operated by the Association of Universities for Research in Astronomy, Inc. (AURA), under cooperative agreement with the National Science Foundation. We finally also thank the referee for their comments on the manuscript.

Appendix A Ionisation and Excitation Balance

In this section we describe how we find the ionisation and excitation balance required for the computation of the line force in Sect. 2.2. Under the LTE assumption, the ionisation balance is computed using the Saha-Boltzmann equations. We obtain the ionisation balance by first computing the number density n_z of a specific element with nuclear charge z, and atomic weight A_z (for hydrogen A_z = m_H). Taking the number abundance relative to hydrogen N_z (relative number abundances in this paper are obtained from the tabulated solar values in Asplund et al. 2009) and the mass fraction of hydrogen

$$X=\frac{m_{\text{H}}}{{\displaystyle {\sum }_z{A}_z{N}_z}},$$(A.1)

the number density of element z for a plasma with mass density ρ is:

$$n_z=\rho \frac{N_{z}X}{m_{\text{H}}}.$$(A.2)

Partition functions of every ionisation stage i of each element z are computed as

$$U_{z,i}={\displaystyle \sum _{l=1}^{I_{\max }}{g}_{\text{l}}\exp \left(\frac{{є}_{zi}-{\epsilon }_{zil}}{k_{\text{B}}T}\right)},$$(A.3)

where g₁ is the statistical weight of level l, ϵ_{z i} is the ionisation energy from the ground level, and ϵ_{z i l} is the excitation energy of the corresponding level. In general, sums are taken over all levels; however, as higher excitation levels generally have a small contribution, we only sum to some highest excitation level l_max (see Appendix В for details on the atomic data used in this paper, including information on maximum ionisation states and excitation levels for each ion). Since we have access to the extensive Munich atomic database, we opt to compute the partition functions using the above sum throughout our work, (rather than using fit functions such as those provided by Cardona et al. 2010). This is the same way as is being done within the well-calibrated, standard ID model atmosphere code FASTWIND (Puls et al. 2005), which also uses the same atomic database as here.

Now the Saha equation is used to compute the number densities of ionisation stage i + 1 relative to i for each element,

$$\frac{n_{zi+1}}{n_{zi}}=\frac{2U_{zi+1}}{n_{\text{e}}{U}_{zi}}{\left(\frac{2\pi m_{\text{e}}{k}_{\text{B}}T}{h^2}\right)}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\exp \left(\frac{{є}_{zi}-{є}_{zi+1}}{k_{\text{B}}T}\right).$$(A.4)

Here n_e is the number density of free electrons, and n_zi and n_zi+1 are the number densities of the i and i + 1 ionisation stages of element z. It is more convenient to express number densities for each ionisation stage with respect to the neutral stage of the element:

$$\frac{n_{zi}}{n_{z1}}={\displaystyle \prod _{j=1}^{i-1}\frac{n_{zj+1}}{n_{zj}},}$$(A.5)

where fractions on the right-hand side of expression are given by Eq. A.4. Knowing n_{z 1} the expression above gives us the number density of each ionisation stage. To find n_{z 1},we simply apply conservation of total number density:

$$n_{z1}=n_z{\left({\displaystyle \sum _{i=1}^{i_{\max }}\frac{n_{zi}}{n_{z1}}}\right)}^{-1}$$(A.6)

where the fractions in the sum are given by Eq. A.5.

Finally, it is important to note that these calculations depend on the number density of free electrons n_e, which in turn depends on the result of the same calculations through

$$n_{\text{e}}={\displaystyle \sum _{z=1}^{z_{\max }}{\displaystyle \sum _{i=1}^{i_{\max }}\left(i-1\right)n_{zi}}.}$$(A.7)

As such, iteration is generally required to determine the correct ionisation balance. In the first step of the iteration, n_e is set by assuming fully ionised hydrogen. This then allows the ionisation balance to be solved for and $n_{\text{e}}^{\text{new}}$ computed using Eq. (A.7). We then update the electron number density in the next iteration by taking the square root of the product of its values from the current and previous iterations (see also LC).

Now applying the Boltzmann excitation law, the occupation number density of element z in excitation level l and ionisation state i with respect to the ground-state occupation number density n_{z i 1} is

$$\frac{n_{zil}}{n_{zi1}}=\frac{g_{\text{l}}}{g_1}\exp \left(\frac{{є}_{zi}-{\epsilon }_{zil}}{k_{\text{B}}T}\right).$$(A.8)

Again applying conservation of number density, we get

$$n_{zil}=n_{zi}\frac{n_{zil}}{n_{zi1}}{\left({\displaystyle \sum _{l=1}^{l_{\max }}\frac{n_{zil}}{n_{zi1}}}\right)}^{-1},$$(A.9)

where the fraction terms come from Eq. A.8. In this way, we compute atomic occupation number densities for each excitation level of every ion contained in our database.

Appendix B Atomic Data

For the study presented in this paper, we used the Munich atomic database originally compiled by Pauldrach et al. (1998, 2001) (see also Puls et al. 2005). This database contains elements from hydrogen to zinc, with the exception of Li, Be, B, and Sc, which are rare and provide minimal effect. For the included elements, the database provides ionisation stages up to VIII. The total number of lines included is 4.2 · 10⁶, for which values of their lower levels are provided. We provide a detailed list of available ionisation stages for the element and the number of available lines for each ionisation stage in Table B.1.

References

All Tables

All Figures

Fig. 1 Force multiplier as computed by evaluating the numerical sum over the discrete set of spectral lines in the triangle-marked green curve. Panel a shows the saturation level indicated with the dashed red curve, a power-law region with the dashed blue line, and the locatton of a transition eotnl witt the dashed black line. Panel b shows the actual fit of two numerical results, where the triangle-marked green line has the same temperature and density as in (a) and the square-marked blue line assumes a different set of temperature and density. The derived fit parameters for each case are indicated in the corresponding legend.

In the text

	Fig. 2 Colour maps of the line-distribution statistical parameters (a) ${\log }_{10}\overline{Q}$, (b) α, and (c) ${\log }_{10}\left(Q_0\text{/}\overline{Q}\right)$ computed for a range of temperatures and densities using the Munich line database. The abscissa here measures the logarithm of the temperature in kK, and the ordinate gives the logarithm of the density in g cm−3.
In the text

	Fig. 3 Colour map showing the radial and temporal variation in the line-distribution parameters. From the left to right we show α, $\overline{Q}$, and Q0 radially normalised to their respective final profiles αfin, ${\overline{Q}}_{\text{fin}}$, and Q0 fin.
In the text

	Fig. 4 Radial profile of the line-distribution parameters of the final snapshot of Mod-0. Plot (a) shows the α parameter and (b) the $\overline{Q}$ parameter (solid red line) and Q0 parameter (dashed black line).
In the text

	Fig. 5 Radial and temporal variation of the scaled Sobolev optical depth τt. (a) Colour map showing the radial and temporal variation in optical depth parameter τt. Extrema found in this map are max(log10 τt) = 3.95 and min(log10 τt) = 0.12. (b) Radial profile of τt from the final snapshot.
In the text

Fig. 6 Comparison of the initial and final velocity profiles, and the radial and temporal variation of the mass-loss rate, (a) Radial velocity profiles of the initial conditions with the dashed orange line and the final snapshot with the solid blue line, (b) Colour map showing the relaxation or the radial profile of the logarithm of the mass-loss rate spanning the range from the initial conditions up to 30 tdyn.

In the text

Fig. 7 Comparison of the mass-los s rates from our grid of models (blue dots) to the mas s-los s values predicted in earlier works by Björklund et al. (2021), Krtička & Kubát (2017), and Vink et al. (2001). We plot the mass-loss-luminosity relations predicted by these authors with dashed orange, dotted-dashed maroon, and dotted black lines, respectively. We also over-plot the linear fit to our results (solid blue line).

In the text

	Fig. 8 Colour maps of the ${\log }_{10}{\overline{Q}}_{\text{LC}}$ computed by LC (a) and differences between the ${\log }_{10}{\overline{Q}}_{\text{LC}}$ and ${\log }_{10}\overline{Q}$ presented in Fig. 2a (b). Maps are given as a function of temperature on the abscissa and density on the ordinate. The ${\log }_{10}{\overline{Q}}_{\text{LC}}$ is reconstructed from the data provided by LC.
In the text