Issue 
A&A
Volume 667, November 2022



Article Number  A113  
Number of page(s)  15  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/202142888  
Published online  15 November 2022 
Method and new tabulations for fluxweighted line opacity and radiation line force in supersonic media^{★}
^{1}
Institute of Astronomy, KU Leuven,
Celestijnenlaan 200D,
3001
Leuven, Belgium
email: luka.poniatowski@kuleuven.be
^{2}
National Solar Observatory,
22 Ohi’a Ku St,
Makawao, HI
96768, USA
^{3}
Bartol Research Institute, Department of Physics and Astronomy, University of Delaware,
Newark, DE
19716, USA
^{4}
University of Iowa,
Iowa City, IA, USA
^{5}
Astronomical Institute Anton Pannekoek, Amsterdam University,
Science Park 904,
1098 XH
Amsterdam, The Netherlands
Received:
11
December
2021
Accepted:
21
April
2022
Context. In accelerating and supersonic media, understanding the interaction of photons with spectral lines can be of utmost importance, especially in an accelerating flow. However, fully accounting for such line forces is computationally expensive and challenging, as it involves complicated solutions of the radiative transfer problem for millions of contributing lines. This currently can only be done by specialised codes in 1D steadystate flows. More general cases and higher dimensions require alternative approaches.
Aims. We present a comprehensive and fast method for computing the radiation line force using tables of spectrallinestrength distribution parameters, which can be applied in arbitrary (multiD, timedependent) simulations, including those that account for the linedeshadowing instability, to compute the appropriate opacities.
Methods. We assume local thermodynamic equilibrium to compute a fluxweighted line opacity from ~4 million spectral lines. We fit the opacity computed from the line list with an analytic result derived for an assumed distribution of the spectral line strength and found the corresponding linedistribution parameters, which we tabulate here for a range of assumed input densities ρ ∈ [10^{−20}, 10^{−10}] g cm^{−3} and temperatures T ∊ [10^{4}, 10^{47}] K.
Results. We find that the variation in the linedistribution parameters plays an essential role in setting the wind dynamics in our models. In our benchmark study, we also find a good overall agreement between the Ostar massloss rates of our models and those derived from steadystate studies that use a more detailed radiative transfer.
Conclusions. Our models reinforce the idea that selfconsistent variation in the linedistribution parameters is important for the dynamics of linedriven flows. Within a wellcalibrated Ostar regime, our results support the proposed methodology. In practice, utilising the provided tables, yielded a factor >100 speedup in computational time compared to specialised 1D modelatmosphere codes of linedriven winds, which constitutes an important step towards efficient multidimensional simulations. We conclude that our method and tables are ready to be exploited in various radiationhydrodynamic simulations where the line force is important.
Key words: stars: earlytype / stars: atmospheres / stars: winds / outflows / stars: massloss / radiative transfer / hydrodynamics
Data used to create Fig. 2 are available in electronic form at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/vizbin/cat/J/A+A/667/A113, and also at https://zenodo.org/record/6505135
© L. G. Poniatowski et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
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 massivestar atmospheres, the outward push from radiation causes the star to lose mass through powerful radiationdriven stellar winds. Prominent spectral features such as PCygni lines provide observational evidence of radiationdriven 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, stateoftheart radiationhydrodynamic (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 multiline radiative transfer problem in a reference frame), have exclusively focused on 1D steadystate cases (e.g. Sander et al. 2017; Sundqvist et al. 2019). On the other hand, multidimensional, timedependent 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 ‘velocitystretch’ 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 socalled Sobolev approximation (Sobolev 1960), assuming that hydrodynamic quantities^{1} are constant over the photonline 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 steadystate) simulations of hotstar winds (see also GormazMatamala et al. 2019).
In timedependent numerical RHD simulations based on a CAK parameterisation, lineforce parameters are typically assumed to take some predescribed values; for example, in Poniatowski et al. (2021) an ad hoc radial stratification for the CAK powerlaw index parameter α was assumed. However, in a more realistic treatment, the parameters describing the line force should rather be selfconsistently 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 (fluxweighted) 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 multidimensional extensions of the 1D WolfRayet (WR) models presented in Moens et al. (2022b) and an extension of the current method to corresponding nonSobolev, observer’s frame lineforce calculations (essentially following Owocki & Puls 1996). This latter will allow us, for the first time, to investigate effects of selfconsistent lineforce parameters on the strong ‘line deshadowing instability’ (LDI) (Owocki & Rybicki 1984); thus far, all such LDI simulations have simply used predescribed, prototypical lineforce 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 fluxweighted line opacity (Sect. 2). We then fit our results from the complete line list to an expression derived from the assumed linestrength distribution function, tabulating them by means of our (three) fitting parameters (Sect. 2). To benchmark our method, we apply our tables in a 1D timedependent RHD showcase study of the calibrated grid of linedriven Ostar 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 FluxWeighted Mean Opacity from Spectral Lines in an Accelerating Supersonic Medium
As mentioned above, using a statistical approach, one may approximate the fluxweighted mean opacity in both the Sobolev and nonSobolev limits (Owocki & Puls 1996). In this paper, we use the Sobolev approximation to evaluate the fluxweighted 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 FluxWeighted Mean Opacity from a Spectral Line in the Sobolev Approximation
Assuming isotropic emission, the radiation acceleration on a medium illuminated by a frequencydependent 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^{−3}, and c is the speed of the light in km s^{−1}. 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 fluxweighted 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}}\xb7F}{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 fluxweighted 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_{0} 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^{2}/(m_{e}c) is the classical frequencyintegrated line crosssection, with e and m_{e} being the charge and mass of the electron, and ϕ_{v} = ϕ(v − v_{0}) is the normalised line profile function, with v_{0} the linecentre frequency; k_{L} then represents the frequencyintegrated lineextinction 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 opticaldepth equation above from physical to frequency space. Then, under the Sobolev approximation (Sobolev 1960), where the lineofsight 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 directionaldependent 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}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\nabla \left(\widehat{n}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}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}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\nabla \left(\widehat{n}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}v\right)\right},$$(13)
where κ_{0} = 0.34 cm^{2} g^{−1} is a nominal mass absorption coefficient set to a constant for convenience of the normalisation^{2}, 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 fluxweighted mean opacity then gives the socalled 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 fluxweighted 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^{6} lines (for more information, see Pauldrach et al. 1998, 2001). For a discrete list of intrinsically nonoverlapping 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 blackbody 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 linecentre 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_{1} is the oscillator strength for a given transition. Computing n_{1} 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^{−12} g cm^{−3}, we show the resulting force multiplier for the range log_{10} τ_{t} = [−15, 5], as a trianglemarked 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_{10} M(τ_{t}) follows a nearly constant slope as a function of log_{10} τ_{t}, which is a wellknown 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 powerlaw exponent and fit parameter. Upon inspection of Fig. 1a, it is clear that such a powerlaw 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 thintothick transition occurs, that is, τ_{t} = τ_{t trans}. To gain insight into what sets this transition, we reexamine 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_{0} (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_{0} ≈ 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 trianglemarked green curve. Panel a shows the saturation level indicated with the dashed red curve, a powerlaw 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 trianglemarked green line has the same temperature and density as in (a) and the squaremarked 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 PowerLaw 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_{0} 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 linedistribution function to be specified in order to compute the corresponding nonSobolev 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 linestrength 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{\Gamma}\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_{0}τ_{t} ≫ 1, we directly obtain the same powerlaw 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 wellknown 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_{0}τ_{t} ≪ 1, a firstorder 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_{0}, 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_{0} 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_{0} 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 socalled δ parameter, introduced by Abbott (1982) to account for the shift of ionisation states over the wind. In cases when the linedistribution 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 linedriven 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 linedistribution parameters. Those models then require the δ parameter to introduce back the radial variation neglected by the use of spatially fixed lineforce parameters (see also GormazMatamala et al. 2019). In our new method presented here, these effects are instead directly accounted for by allowing the linedistribution 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^{−17} g cm^{−3}, along the corresponding bestfitting 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 Ostar 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 trianglemarked green curve in Fig. 1a are rather good representations of the expected values for the already quite extensively analysed regime of early Osupergiants. 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 Osupergiant regime, the best fit parameters α = 0.48, $\overline{Q}=35$ and Q_{0} = 2634 are very much outside these previously estimated ranges. This illustrates a general need for not assuming constant lineforce parameters independent of conditions, which so far has often been the case in (at least timedependent) simulations of linedriven flows. We will thus now expand the above calculations towards a much broader range of temperature and density combinations.
2.4 Grid of LineDistribution Parameters
We apply the method described above to various combinations of temperature and density in the range log_{10}(ρ [g^{−1}cm^{3}]) ∈ [−20, −10] and log_{10}(T [K^{−1}]) ∊ [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_{10} $\overline{Q}$, ${\mathrm{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 Otype stars in Sect. 3. Typical surface conditions in O stars cover the region of log_{10}(ρ [g^{−1} cm^{3}]) ∈ [−14, −10] and log_{10}(T [K^{−1}]) ∈ [4.3, 4.7]. Inspecting this region in Figs. 2a, b, we find values of α ≈ 0.5–0.7 and $\overline{Q}\approx 1000\u20132000$ that again are in good agreement with previous estimates (e.g. Gayley 1995; Puls et al. 2000). Moreover, for the Ostar 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 Ostar 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 lineforce 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 fluxweighted 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 fluxweighted 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\Omega}\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_{0} 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 timestep, we find appropriate values of α, $\overline{Q}$, Q_{0}, and τ_{t} based on the local density, temperature, and lineofsight velocity gradient. In practice, values of α, $\overline{Q}$, and Q_{0} are found by interpolations in our predescribed 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 fluxweighted line opacity in the Sobolev limit.
However, in Sect. 4.2.2 we discuss a possible extension also towards nonSobolev (observer’s frame) applications.
Fig. 2 Colour maps of the linedistribution statistical parameters (a) ${\mathrm{log}}_{10}\overline{Q}$, (b) α, and (c) ${\mathrm{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, selfconsistent massivestar 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 (GormazMatamala et al. 2021) are highly timeconsuming. Those computations are also exclusively 1D and only consider spherically symmetric and steadystate wind outflows. Compared to these steadystate models, relaxing a 1D Ostar 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 timedependent 2D and 3D simulations (something that is computationally prohibited when using the timeconsuming radiative transfer methods applied in the codes mentioned above). Furthermore, in comparison to previous timedependent linedriven wind models (see e.g. Feldmeier & Thomas 2017; Dyda & Proga 2018), which typically rely on lineforce parameters that are predescribed and fixed throughout the simulations (i.e. both in space and time), our new models allow these parameters to vary in space and time, selfconsistently 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 Ostar 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 lineforce parameters, the Ostar grid allows us to further validate our technique by comparing massloss rates found in our simulations to the values predicted for such stars by other studies using more detailed radiative transfer.
Resulting massloss rates Ṁ from our hydrodynamic models for a grid of Ostar models with stellar parameters ℒ, ℳ, ℛ, and T_{eff}.
3.1 Numerical Simulations of LineDriven Ostar 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^{−1} , g = ℳG/r^{2} is the gravitational acceleration in cm s^{−2}, and P is the gas pressure in g cm^{−1} s^{−2}. Using the ideal gas law, P = a^{2}ρ, where a^{2} = k_{B}T/(μ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}=\mathcal{L}\text{/}\left(4\pi {\sigma}_{\text{SB}}{R}^{2}\right)$ is the effective temperature of the star^{3}. Finally, κ_{e}F_{r}/c and κ_{Fr,r}F_{r}/c are the two components of radiation acceleration from Thompson scattering and spectral lines, respectively, with F_{r} = ℒ/(4πr^{2}) 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 ρκ_{0}/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_{0}τ_{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}\mathrm{ln}v}{\text{d}\mathrm{ln}r}1.$$
For simplicity, we introduce the effective gravity g_{eff}(r) = ℳG[1 − Γ_{e}(r)]/r^{2}, 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 pretabulated 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 nonconstant 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=\mathrm{min}\left(\text{d}{t}_{\text{CFL}},\text{CFL}\times \sqrt{\frac{c\text{\Delta}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 linedistribution parameters. From the left to right we show α, $\overline{Q}$, and Q_{0} radially normalised to their respective final profiles α_{fin}, ${\overline{Q}}_{\text{fin}}$, and Q_{0 fin}. 
3.2 Initial and Boundary Conditions
We set the initial velocity using the socalled βlaw:
$${v}_{\beta}\left(r\right)={v}_{\infty}{\left(1b\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^{−1}. 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 timedependent 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 selfadjust) to the overlaying wind by updating the density of the first ghost сell ρ_{1} 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})^{2} = a^{2}, 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 quasihydrostatic equilibrium for the density, where after the velocity is determined from most conservation. In comparison to our previous models with ρ_{1} 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 Mod0
Following the above setup, we first present a numerical model (Mod0) 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 linedistribution parameters are interpulated^{4} 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_{0}(ρ, T) at each radius and each timestep 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_{0} set the scale of the line force, it is natural that the radial and temporal variation in these parameter then also reflects upon the selfconsistent 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) steadystate ‘final’ configuration is reached much earlier than this. Because the relaxed radial profiles of α, $\overline{Q}$, and Q_{0} 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^{5}). 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 outerwind 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_{0} 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_{0} 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 Ostar 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_{0} and the line force form a negative feedback loop, where increasing Q_{0} rather decreases the line opacity. A selfregulatory system forms through the interaction of these positive and negative feedback loops, which ultimately here determines the final massloss 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 linedistribution parameters (Fig. 3), we observe a transient behaviour from the initial condition, which takes about 40 t_{dyn} to fully vanish. After this, a nearequilibrium configuration is reached. Importantly, the value of τ_{t} during the simulation never exceeds the range initially used to compute the linedistribution parameters. As such, our tabulations always cover the correct range in τ_{t} (as well as in density, velocity gradient and temperature)^{6}. 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_{0} τ_{t} ≫ 1, and as such the full Ostar wind resides in the powerlaw portion of M(τ_{t}), as also previously noted, for instance by GormazMatamala 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 timedependent and clumpy with coexisting regions of (much) higher and lower densities. In such cases, very low values of τ_{t} would systematically be found in the lowdensity 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^{−1} 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^{−1}. 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_{0} will also affect the terminal velocity. To differentiate between effects introduced by the finite disc factor from the effects introduced by variation in linedistribution parameters, we also compute a ‘reference’ model, which uses the same stellar parameters as the Mod0, but where we keep linedistribution 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 linedistribution parameters, this stems directly from their specific radial behaviour.
Further, to demonstrate that the wind velocity and density also relax as linedistribution parameters reach their final configuration, in Fig. 6b we plot a colour map of radial and temporal variation in the massloss rate zooming in on the first 30 t_{dyn}. The axes in this figure are the same as in Fig. 3. Like for the lineforce parameters the relaxation takes around 20 t_{dyn} after which the model relaxes to a final selfconsistent mass loss of Ṁ = 1.8 × 10^{−6}M_{⊙}yr^{−1}. 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 massloss results reemphasises the general quality of the method. Indeed, as now shown, our computed massloss rates agree quite well with the predictions by Björklund et al. (2021) and Krtička & Kubát (2017) in the complete Ostar domain.
Fig. 4 Radial profile of the linedistribution parameters of the final snapshot of Mod0. Plot (a) shows the α parameter and (b) the $\overline{Q}$ parameter (solid red line) and Q_{0} 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(log_{10} τ_{t}) = 3.95 and min(log_{10} τ_{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 massloss 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 massloss rate spanning the range from the initial conditions up to 30 t_{dyn}. 
3.4 Grid of O Stars
Using the same setup as for Mod0, 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_{0} decrease. Near the lower boundary, we find values of all three lineforce 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 Ostar models, we find ${Q}_{0}\approx \overline{Q}$ only in the nearsurface 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_{0}τ_{t} of the line ensemble in lowdensity regions, which can then lead to an inconsistent radial line force.
Furthermore, although the typical conviction is that nonLTE (NLTE) effects are essential in setting the dynamics of Ostar winds, our modelgrid analysis seems to suggest the opposite, namely that an LTEbased line force (explicit assumption in Sect. 2.2) sufficiently well captures the fundamental dynamics of (at least) the nearstar region, which essentially sets the massloss rate of the model. This suggestion stems from the (on average) tight agreement of the massloss 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 massloss rates as a function of luminosity. In addition, we performed a linear fit plotted and compare this masslossluminosity 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 masslossluminosity 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 Ostar regime is (nowadays) quite wellestablished (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 comoving 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 massloss 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 radiationforce dip in nearsonic regions observed in nonSobolev models (Björklund et al. 2021).
Interestingly, our results thus suggest that this effect is more significant for lower luminosities and massloss rates, as already pointed out by Owocki & Puls (1999). Furthermore, Krtička & Kubát (2017) scale their computed nonSobolev 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 Sobolevbased’ and the nature of their steadystate 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 Ostar 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 wellstudied Ostar 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 massloss 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 steadystate 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 multidimensional RHD setups (e.g. Moens et al. 2022a).
Fig. 7 Comparison of the masslos s rates from our grid of models (blue dots) to the mas slos 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 masslossluminosity relations predicted by these authors with dashed orange, dotteddashed maroon, and dotted black lines, respectively. We also overplot 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 linedistribution 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 multidimensional RHD simulations (see the discussion in Sect. 3) and towards situations wherein a nonSobolev 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 renormalise 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 ${\mathrm{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 ${\mathrm{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 smallfactor difference throughout most of the parameter space. However, an exception is the lowdensity and hightemperature 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 Ostar 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 Ostar winds, the maximum difference between the two calculations is less than 20% in $\overline{Q}$.
Fig. 8 Colour maps of the ${\mathrm{log}}_{10}{\overline{Q}}_{\text{LC}}$ computed by LC (a) and differences between the ${\mathrm{log}}_{10}{\overline{Q}}_{\text{LC}}$ and ${\mathrm{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 ${\mathrm{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 NonPlanckian 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 Ostar massloss 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 followup investigation to this study will be to examine effects from this also upon the detailed lineforce 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 linedistribution 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 fluxweighted mean opacity. In particular, if lines overlap within the velocity domain of an outflow, this can lead to multiple lineresonances, 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 highmass Xray 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 lineforce sum, Eq. (15), neglects all intrinsic lineoverlap effects. This may act towards reducing the line force due to selfshadowing effects within intrinsically overlapping optically thick lines.
4.2.2 NonSobolev Line Force
In Sect. 2.3, we mentioned how a key advantage of the lineforce parametrisation used here is that, it yields an explicit form for the underlying statistical distribution of lines. Using this statistical linedistribution function (i.e. Eq. (23)), one can then also compute a corresponding nonSobolev line force based on our tabulations, for example, by following the observer’s frame escapeintegral methods outlined in Owocki & Puls (1996). This would then allow us to investigate the effects from varying the lineforce 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 lineforce 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 lineforce parameters can be investigated in an analogous way as in this paper. Thus far, LDI simulations have always assumed prototypical lineforce parameters (constant both in time and space), although it is clear that varying them will also impact the resulting LDIgenerated structure (Driessen et al. 2019). First simulations of the LDI using selfconsistent lineforce parameters are already underway and will be presented in a forthcoming paper.
5 Conclusions
This paper has presented a comprehensive tabulation of linedistribution parameters that can be used in the fast calculation of fluxweighted line opacities and line forces in supersonic media (Sect. 2). In Sect. 3, we demonstrate that, when applied to the (rather wellcalibrated) linedriven flows of О stars, our method, considering the limitations of the method outlined in Sect. 4.2, gives a very good agreement between our derived massloss rates and those obtained in the detailed steadystate 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 selfconsistent variation in the linestrength distribution parameters (i.e. α, $\overline{Q}$, and Q_{0}) have critical feedback on the resulting linedriven 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 multidimensional, timedependent linedriven 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 nearstar regions of О stars, where their massloss rates are determined. Nonetheless, a key followup 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 fluxweighted line opacities is a beneficial, versatile, and very fast tool for implementing more realistic treatments of radiation line forces in multidimensional, timedependent RHD studies. For example, in a unified simulation encompassing both the optically thick, deep, subsonic regions and an overlying linedriven wind, the method here can be trivially applied within the hybridopacity 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 timedependent 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 SahaBoltzmann 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}_{\mathrm{max}}}{g}_{\text{l}}\mathrm{exp}\left(\frac{{\u0454}_{z\text{\hspace{0.17em}}i}{\epsilon}_{z\text{\hspace{0.17em}}i\text{\hspace{0.17em}}l}}{{k}_{\text{B}}T}\right)},$$(A.3)
where g_{1} 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 wellcalibrated, 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}_{z\text{\hspace{0.17em}}i+1}}{{n}_{z\text{\hspace{0.17em}}i}}=\frac{2{U}_{z\text{\hspace{0.17em}}i+1}}{{n}_{\text{e}}{U}_{z\text{\hspace{0.17em}}i}}{\left(\frac{2\pi {m}_{\text{e}}{k}_{\text{B}}T}{{h}^{2}}\right)}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{1ex}{$2$}\right.}\mathrm{exp}\left(\frac{{\u0454}_{z\text{\hspace{0.17em}}i}{\u0454}_{z\text{\hspace{0.17em}}i+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}_{z\text{\hspace{0.17em}}i}}{{n}_{z\text{\hspace{0.17em}}1}}={\displaystyle \prod _{j=1}^{i1}\frac{{n}_{z\text{\hspace{0.17em}}j+1}}{{n}_{z\text{\hspace{0.17em}}j}},}$$(A.5)
where fractions on the righthand 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}_{z\text{\hspace{0.17em}}1}={n}_{z}{\left({\displaystyle \sum _{i=1}^{{i}_{\mathrm{max}}}\frac{{n}_{z\text{\hspace{0.17em}}i}}{{n}_{z\text{\hspace{0.17em}}1}}}\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}_{\mathrm{max}}}{\displaystyle \sum _{i=1}^{{i}_{\mathrm{max}}}\left(i1\right){n}_{z\text{\hspace{0.17em}}i}}.}$$(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 groundstate occupation number density n_{z i 1} is
$$\frac{{n}_{z\text{\hspace{0.17em}}i\text{\hspace{0.17em}}l}}{{n}_{z\text{\hspace{0.17em}}i\text{\hspace{0.17em}}1}}=\frac{{g}_{\text{l}}}{{g}_{1}}\mathrm{exp}\left(\frac{{\u0454}_{z\text{\hspace{0.17em}}i}{\epsilon}_{z\text{\hspace{0.17em}}i\text{\hspace{0.17em}}l}}{{k}_{\text{B}}T}\right).$$(A.8)
Again applying conservation of number density, we get
$${n}_{z\text{\hspace{0.17em}}i\text{\hspace{0.17em}}l}={n}_{z\text{\hspace{0.17em}}i}\frac{{n}_{z\text{\hspace{0.17em}}i\text{\hspace{0.17em}}l}}{{n}_{z\text{\hspace{0.17em}}i\text{\hspace{0.17em}}1}}{\left({\displaystyle \sum _{l=1}^{{l}_{\mathrm{max}}}\frac{{n}_{z\text{\hspace{0.17em}}i\text{\hspace{0.17em}}l}}{{n}_{z\text{\hspace{0.17em}}i\text{\hspace{0.17em}}1}}}\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^{6}, 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.
Available ionisation stages for elements from H to Zn, and the number of available lines for each ionisation stage included in the database used in our study.
References
 Abbott, D. C. 1982, ApJ, 259, 282 [Google Scholar]
 Asplund, M., Grevesse, N., Sauvai, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Björklund, R., Sundqvist, J. O., Puls, J., & Najarro, F. 2021, A&A, 648, A36 [EDP Sciences] [Google Scholar]
 Cardona, O., MartínezArroyo, M., & LópezCastillo, M. A. 2010, ApJ, 711, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I. 2004, Radiation Hydrodynamics (Cambridge, UK: Cambridge University Press) [CrossRef] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [Google Scholar]
 Cohen, D. H., Wollman, E. E., Leutenegger, M. A., et al. 2014, MNRAS, 439, 908 [NASA ADS] [CrossRef] [Google Scholar]
 Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32 [Google Scholar]
 Cranmer, S. R. 1996, PhD thesis, Bartol Research Institute, University of Delaware, USA [Google Scholar]
 Dannen, R. C., Proga, D., Kallman, T. R., & Waters, T. 2019, ApJ, 882, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, S. W., Stone, J. M., & Jiang, Y.F. 2012, ApJS, 199, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Driessen, F. A., Sundqvist, J. O., & Kee, N. D. 2019, A&A, 631, A172 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Driessen, F. A., Kee, N. D., & Sundqvist, J. O. 2021, A&A, 656, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dyda, S., & Proga, D. 2018, MNRAS, 481, 5263 [NASA ADS] [CrossRef] [Google Scholar]
 Feldmeier, A., & Thomas, T. 2017, MNRAS, 469, 3102 [NASA ADS] [CrossRef] [Google Scholar]
 Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701 [Google Scholar]
 Friend, D. B., & Castor, J. I. 1983, ApJ, 272, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Garcia, M., Herrero, A., Najarro, F., Lennon, D. J., & Alejandro Urbaneja, M. 2014, ApJ, 788, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [Google Scholar]
 Goldberg, J. A., Jiang, Y.F., & Bildsten, L. 2021, ApJ, 929, 156 [Google Scholar]
 GormazMatamala, A. C., Curé, M., Cidale, L. S., & Venero, R. O. J. 2019, ApJ, 873, 131 [NASA ADS] [CrossRef] [Google Scholar]
 GormazMatamala, A. C., Curé, M., Hillier, D. J., et al. 2021, ApJ, 920, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Hawcroft, C., Sana, H., Mahy, L., et al. 2021, A&A, 655, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Higginbottom, N., Proga, D., Knigge, C., et al. 2014, ApJ, 789, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, Y.F., Cantiello, M., Bildsten, L., et al. 2018, Nature, 561, 498 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, Y.F., Stone, J. M., & Davis, S. W. 2019, ApJ, 880, 67 [Google Scholar]
 Kee, N. D., Owocki, S., & Sundqvist, J. O. 2016, MNRAS, 458, 2323 [Google Scholar]
 Krticka, J. 2006, MNRAS, 367, 1282 [NASA ADS] [CrossRef] [Google Scholar]
 Krticka, J., & Kubát, J. 2017, A&A, 606, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lattimer, A. S., & Cranmer, S. R. 2021, ApJ, 910, 48 [CrossRef] [Google Scholar]
 Lucy, L. B., & Solomon, P M. 1970, ApJ, 159, 879 [Google Scholar]
 Moens, N., Poniatowski, L. G., Hennicker, L., et al. 2022a, A&A, 665, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moens, N., Sundqvist, J. O., El Mellah, I., et al. 2022b, A&A, 657, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Najarro, F., Hanson, M. M., & Puls, J. 2011, A&A, 535, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owocki, S. P. 2015, in Astrophysics and Space Science Library, Very Massive Stars in the Local Universe, ed. J. S. Vink, 412, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & Puls, J. 1996, ApJ, 462, 894 [Google Scholar]
 Owocki, S. P., & Puls, J. 1999, ApJ, 510, 355 [Google Scholar]
 Owocki, S. P., & Rybicki, G. B. 1984, ApJ, 284, 337 [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] [Google Scholar]
 Pauldrach, A. W. A., Lennon, M., Hoffmann, T. L., et al. 1998, in Astronomical Society of the Pacific Conference Series, Properties of Hot Luminous Stars, ed. I. Howarth, 131, 258 [NASA ADS] [Google Scholar]
 Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Petrenz, P., & Puls, J. 2000, A&A, 358, 956 [NASA ADS] [Google Scholar]
 Poniatowski, L. G., Sundqvist, J. O., Kee, N. D., et al. 2021, A&A, 647, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Proga, D., Stone, J. M., & Drew, J. E. 1998, MNRAS, 295, 595 [Google Scholar]
 Puls, J., Owocki, S. P., & Fullerton, A. W. 1993, A&A, 279, 457 [NASA ADS] [Google Scholar]
 Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Vink, J. S., & Najarro, F. 2008, A&A Rev., 16, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Sander, A. A. C., Hamann, W. R., Todt, H., Hainich, R., & Shenar, T. 2017, A&A, 603, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sander, A. A. C., Vink, J. S., & Hamann, W. R. 2020, MNRAS, 491, 4406 [Google Scholar]
 Schaerer, D., & Schmutz, W. 1994, A&A, 288, 231 [NASA ADS] [Google Scholar]
 Schrøder, S. L., MacLeod, M., RamirezRuiz, E., et al. 2021, ApJ, submitted, [arXiv:2107.09675] [Google Scholar]
 Sobolev, V. V. 1960, Moving Envelopes of Stars (Harvard University Press) [Google Scholar]
 Sundqvist, H., & Veronis, G. 1970, Tellus, 22, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., Puls, J., Feldmeier, A., & Owocki, S. P. 2011, A&A, 528, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Owocki, S. P., & Puls, J. 2018, A&A, 611, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Björklund, R., Puls, J., & Najarro, F. 2019, A&A, 632, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Šurlan, B., Hamann, W. R., Aret, A., et al. 2013, A&A, 559, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Leer, B. 1979, J. Comput. Phys., 32, 101 [Google Scholar]
 Vink, J. S. 2022, ARA&A, 60, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
In principle, this difference between T and T_{eff} introduces a small inconsistency with the T = T_{rad} assumption that went into the calculation of the tables. However, from test calculations that also use T = T_{eff}, we have verified that this does not cause any significant differences for the results presented in this section.
Although we also note that an earlier study by Schaerer & Schmutz (1994) does seem to rather imply a radial decline of α within a stationary wind; the way of practically deriving the actual values for α in that study is, however, quite different from here, and as such it is difficult to assess where this difference in results might come from.
Data were obtained from the online repository provided by the authors: https://github.com/aslyv2/RadWinds
All Tables
Resulting massloss rates Ṁ from our hydrodynamic models for a grid of Ostar models with stellar parameters ℒ, ℳ, ℛ, and T_{eff}.
Available ionisation stages for elements from H to Zn, and the number of available lines for each ionisation stage included in the database used in our study.
All Figures
Fig. 1 Force multiplier as computed by evaluating the numerical sum over the discrete set of spectral lines in the trianglemarked green curve. Panel a shows the saturation level indicated with the dashed red curve, a powerlaw 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 trianglemarked green line has the same temperature and density as in (a) and the squaremarked 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 linedistribution statistical parameters (a) ${\mathrm{log}}_{10}\overline{Q}$, (b) α, and (c) ${\mathrm{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 linedistribution parameters. From the left to right we show α, $\overline{Q}$, and Q_{0} radially normalised to their respective final profiles α_{fin}, ${\overline{Q}}_{\text{fin}}$, and Q_{0 fin}. 

In the text 
Fig. 4 Radial profile of the linedistribution parameters of the final snapshot of Mod0. Plot (a) shows the α parameter and (b) the $\overline{Q}$ parameter (solid red line) and Q_{0} 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(log_{10} τ_{t}) = 3.95 and min(log_{10} τ_{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 massloss 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 massloss rate spanning the range from the initial conditions up to 30 t_{dyn}. 

In the text 
Fig. 7 Comparison of the masslos s rates from our grid of models (blue dots) to the mas slos 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 masslossluminosity relations predicted by these authors with dashed orange, dotteddashed maroon, and dotted black lines, respectively. We also overplot the linear fit to our results (solid blue line). 

In the text 
Fig. 8 Colour maps of the ${\mathrm{log}}_{10}{\overline{Q}}_{\text{LC}}$ computed by LC (a) and differences between the ${\mathrm{log}}_{10}{\overline{Q}}_{\text{LC}}$ and ${\mathrm{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 ${\mathrm{log}}_{10}{\overline{Q}}_{\text{LC}}$ is reconstructed from the data provided by LC. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.