Issue 
A&A
Volume 684, April 2024



Article Number  A177  
Number of page(s)  21  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/202348206  
Published online  23 April 2024 
2D unified atmosphere and wind simulations of Otype stars
Instituut voor Sterrenkunde, KU Leuven,
Celestijnenlaan 200D,
3001
Leuven,
Belgium
email: dwaipayan.debnath@kuleuven.be
Received:
9
October
2023
Accepted:
13
January
2024
Context. Massive and luminous Otype star (O star) atmospheres with winds have been studied primarily using onedimensional (1D), spherically symmetric, and stationary models. However, observations and theory have suggested that O star atmospheres are highly structured, turbulent, and timedependent. As such, when making comparisons to observations, presentday 1D modeling tools require the introduction of ad hoc quantities such as photospheric macro and microturbulence, wind clumping, and other relevant properties.
Aims. We present a series of multidimensional, timedependent, radiationhydrodynamical (RHD) simulations for O stars that encapsulate the deeper subsurface envelope (down to T ~ 450 kK), as well as the supersonic linedriven wind outflow in one unified approach. Our overarching aim is to develop a framework that is free from the adhoc prescriptions that plague presentday 1D models. Here, we start with an analysis of a small set of such multidimensional simulations and then compare them to atmospheric structures predicted by their 1D counterparts.
Methods. We performed timedependent, twodimensional (2D) simulations of O star atmospheres with winds using a fluxlimiting RHD finite volume modelling technique. Opacities are computed using a hybrid approach combining tabulated Rosseland means with calculations (based on the Sobolev approximation) of the enhanced line opacities expected for supersonic flows. The initial conditions and comparison models were derived using similar procedures as those applied in standard 1D stationary model atmosphere with wind codes.
Results. Structure starts appearing in our simulations just below the ironopacity peak at ~200 kK. Local pockets of gas with radiative accelerations that exceed gravity then shoot up from these deep layers into the upper atmosphere, where they interact with the linedriven wind outflow initiated around or beyond the variable photosphere. This complex interplay creates large turbulent velocities in the photospheric layers of our simulations, on the order of ~30–100km s^{−1}, with higher values for models with higher luminositytomass ratios. This offers a generally good agreement with observations of large photospheric ‘macroturbulence’ in O stars. When compared to 1D models, the average structures in the 2D simulations display less envelope expansion and no sharp densityinversions, along with density and temperature profiles that are significantly less steep around the photosphere, and a strong anticorrelation between velocity and density in the supersonic wind. Although the wind initiation region is complex and highly variable in our simulations, our average massloss rates agree well with stationary wind models computed by means of full comoving frame radiative transfer solutions.
Conclusions. The different atmospheric structures found in 2D and 1D simulations are likely to affect the spectroscopic determination of fundamental stellar and wind parameters for O stars as well as the empirical derivation of their chemical abundance patterns. To qualitatively match the different density and temperature profiles seen in our multidimensional and 1D models, we need to add a modest amount of convective energy transport in the deep subsurface layers and a large turbulent pressure around the photosphere to the 1D models.
Key words: hydrodynamics / instabilities / methods: numerical / stars: atmospheres / stars: massive / stars: winds, outflows
© The Authors 2024
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 Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Due to the influence of a strong linedriven outflow (Castor et al. 1975), atmospheric models of massive Otype stars (O stars) are required to treat the subsurface layers and the overlying wind simultaneously (Gabler et al. 1989). While such ‘unified’ model atmosphere tools have since reached a high level of maturity to date, they have all been computed under the assumption of a 1D, spherically symmetric, and stationary atmosphere (e.g. FASTWIND, SantolayaRey et al. 1997; Puls et al. 2005, 2020; Sundqvist & Puls 2018; CMFGEN, Hillier & Miller 1998; Hillier & Lanz 2001; POWR, Gräfener et al. 2002; Hamann & Gräfener 2004; Sander et al. 2012). On the other hand, observations and theory strongly suggest that the coupled envelopes, atmospheres, and winds of massive O stars are highly structured and variable, involving convective and radiative envelope and wind instabilities (Hearn 1972; Owocki & Rybicki 1984; Glatzel 1994; Blaes & Socrates 2003; Cantiello et al. 2009; Jiang et al. 2015) as well as various observed phenomena such as photospheric ‘macroturbulence’ (Conti & Ebbets 1977; Howarth et al. 1997; SimónDíaz et al. 2017) and ‘wind clumping’ (Eversberg et al. 1998; Puls et al. 2006; Oskinova et al. 2007; Sundqvist et al. 2010; Šurlan et al. 2012).
Recent timedependent, radiationhydrodynamical (RHD) simulations of the socalled ‘ironopacity peak’ region (Rogers & Iglesias 1992; Stothers & Chin 1993) below the surface of massive stars have indeed displayed such structured and turbulent stellar envelopes (Jiang et al. 2015, 2018). However, since these models did not treat the effects of enhanced lineopacity in a medium with significant velocity gradient (the ‘linedriving’ effect, Castor et al. 1975), they could not simulate the launch of the supersonic wind outflow from around the stellar surface. In Poniatowski et al. (2022) we developed an opacityformalism able to treat the subsurface layers and the linedriven wind in a unified way and then applied this in multidimensional simulations of the atmospheres and outflows of classical WolfRayet stars (Moens et al. 2022a).
In this work, we present the first results from our subsequent attempts to build a unified multidimensional atmosphere and wind modeling framework for O stars. To do so, we started from typical 1D and stationary atmospheric and wind models, extended these downwards to well below the critical ironopacity peak region, and then perturbed them and used them as initial conditions to solve the timedependent RHD equations. We use the MPIAMRVAC tool (Xia et al. 2018; Keppens et al. 2023), which solves partial differential equations (PDE) using a finite volume solver under the parallelised MPI framework. We applied the RHD module of MPIAMRVAC (Moens et al. 2022b), and used the ‘Munich Atomic Database’ (Pauldrach et al. 1998; Puls et al. 2000) to calculate selfconsistent lineopacities following Poniatowski et al. (2022). In this paper, our analysis is focused on three 2D simulations of prototypical O stars in the Galaxy, including first comparisons to the subsurface and atmospheric structures of standard 1D models. In followup work, we plan to use our new framework for extension toward full 3D simulations, to make various improvements to our basic modeling tools, for calculation of synthetic spectra, to further benchmark and improve corresponding 1D tools, and to perform direct comparisons to observations.
In Sect. 2, we go through the basic RHD equations and opacityformalism that have been used in our modeling as well as describe the 1D initial conditions along with our prescription of perturbations. In Sect. 3, we present our simulation results. We explore the formation of structures in our prototypical O star models, discuss their general atmospheric and wind properties, as well as make comparisons between the different simulations. In Sect. 4, we then compare the average properties of our 2D simulations with comparable 1D O star models. In Sect. 5, we discuss our main results and further improvements to our methods. Section 6 concludes our work with a summary and a future outlook.
2 Simulation methods
2.1 Radiationhydrodynamics
The previous studies by Moens et al. (2022b,a) developed radiationhydrodynamic (RHD) modules of the general hydrodynamics code MPIAMRVAC (Xia et al. 2018) and used these to study multidimensional WolfRayet stellar atmospheres and winds. Following their model and study, we also solve numerically the timedependent RHD equations on a finite volume mesh with a ‘boxinastar’ setup, including corrections for spherical divergence terms (Sundqvist et al. 2018; Moens et al. 2022b,a). We provide the further below. This process involves solving the mass, momentum, and energy equations for the gas: $${\partial}_{t}\rho +\nabla \cdot (\rho v)=0,$$(1) $${\partial}_{t}(v\rho )+\nabla \cdot \left(v\rho v+{p}_{\text{g}}\mathbb{I}\right)={f}_{\text{grav}}+{f}_{\text{rad}},$$(2) $${\partial}_{t}{e}_{\text{g}}+\nabla \cdot \left({e}_{\text{g}}v+{p}_{\text{g}}v\right)={f}_{\text{grav}}\cdot v+{f}_{\text{rad}}\cdot v+\dot{q}.$$(3)
Here, v is the velocity of the gas, ρ is the gas density, p_{g} is the gas pressure, 𝕀 is a unit tensor, and e_{g} is the total gas energy which comprises internal and kinetic components given by, $${e}_{\text{g}}=\frac{{p}_{\text{g}}}{{\gamma}_{\text{g}}1}+\frac{\rho {v}^{2}}{2},$$(4)
where we neglect any ionisation effects in the equation of state of the gas by assuming a constant adiabatic index γ_{g} = 5/3. The source terms on the righthand side of Eq. (2), ƒ_{grav} and ƒ_{rad}, are the gravitational and radiation force respectively. The heating and cooling of the gas due to its interaction with radiation is given by $\dot{q}$. Since $\dot{q}$ and ƒ_{rad} are dependent on the radiation field, we need to treat radiation and its coupling with the gas formally. To do so, we write the energy equation for the radiation field^{1} as: $${\partial}_{t}{E}_{\text{rad}}+\nabla \cdot \left({E}_{\text{rad}}v+{F}_{\text{diff}}\right)=\dot{q}\nabla v\text{\hspace{0.17em}}:{P}_{\text{rad}}.$$(5)
Here, $E\equiv (4\sigma /c){T}_{\text{rad}}^{4}={a}_{r}{T}_{\text{rad}}^{4}$ represents the frequencyintegrated radiation energy density for a radiation temperature T_{rad}. Also, P_{rad} is the frequencyintegrated radiation pressure tensor and a_{r} is the radiation constant; E is related to the mean intensity J as E = (4π/c)J and to the radiation pressure via the Eddington tensor ƒ = P_{rad}/E. Additionally, ∇v : P_{rad} on the righthand side of Eq. (5) is the dyadic product between the gradient of the velocity vector and P_{rad}. In our formulation, the conservation equations above are solved on a Cartesian grid as the multigrid method implemented to solve the diffusive part of the radiation equation is not adapted for spherical meshes (Teunissen & Keppens 2019; Moens et al. 2022a). As such, in the present extended simulations, we need to modify certain terms to account for sphericity effects on the fluxes. Specifically, a direction (here x) in the Cartesian setup is taken to be the corresponding radial (r) direction, and terms like the divergence operator in this direction are then modified to account for sphericity by adding appropriate source terms to the conservation equations. The extent of the tangential directions (here only one) is then assumed to be small so that effects of curvature can be neglected and fluxes in the tangential directions do not need to be corrected. The method thus effectively becomes equivalent to solving the equations on a spherical mesh, but neglecting curvature effects within a small interval in the tangential direction. The exact source terms that are added for individual equations are given in Appendix A of Moens et al. (2022a).
The radiative flux F_{diff} is derived from $${\mathit{F}}_{\text{diff}}=\frac{c\lambda}{\kappa \rho}\nabla E=D\nabla E,$$(6)
where κ is the fluxweighted opacity, c the speed of light, and λ a fluxlimiter obtained from an approximation that bridges the analytic and opposite optically thick and freestreaming limits, also providing the required Eddington tensor, and the second equality introduces the radiative diffusion coefficient D (see Eqs. (13)– (20) and discussion in Moens et al. 2022b). The radiative force density is $${f}_{\text{rad}}=\rho \frac{\kappa {\mathit{F}}_{\text{diff}}}{c}.$$(7)
As in Moens et al. (2022b,a), throughout this paper, we assume the flux, Planck and energyweighted mean opacities to be the same, given by Eq. (10). Thus $$\dot{q}=c\kappa \rho E4\pi \kappa \rho B,$$(8)
where $B\equiv \sigma {T}_{\text{g}}^{4}/\pi $ is the frequencyintegrated Planck function for gas temperature T_{g} = p_{g}µm_{H}/(k_{b}ρ). We assume a constant mean molecular weight µ = 0.61 and our method allows for different radiation and gas temperatures. k_{b} is the Boltzmann constant and m_{H} is the mass of a proton. Here, the gravitational acceleration is described by: $${f}_{\text{grav}}=\rho \frac{G{M}_{\star}}{{r}^{2}}\widehat{r}={g}_{\text{g}}\widehat{r}\text{,}$$(9)
meaning we have neglected a small modeled envelope and wind mass in comparison to the underlying stellar mass M_{★} (a good assumption for the O stars considered here).
2.2 Hybrid opacity and finite disk correction factor
To account properly for the ‘linedriving’ (Castor et al. 1975) effect in the moving parts of O star atmospheres and winds, we use the ‘hybrid opacity’ technique introduced by Poniatowski et al. (2021; see also Appendix A for connection to Castor (2004)’s ‘expansion opacity’). This approximates the total opacity of the medium as the sum of: $$\kappa ={\kappa}^{\text{Ross}}+{\kappa}^{\text{line}}.$$(10)
Here, κ is the total opacity of the medium, κ^{line} is the opacity due to lines in a supersonic medium (linedriving effect), and κ^{Ross} are the Rosseland mean opacities (in this paper κ^{Ross} = κ^{OPAL} as taken from tabulations by the ‘OPAL’ project, Iglesias & Rogers 1996).
Following Poniatowski et al. (2022), we compute κ^{line} under the Sobolev (1960) and local thermodynamical equilibrium (LTE) approximations, summing over all lines present in the ‘Munich Atomic Database’ (in total ~4 × 10^{6} lines, Pauldrach et al. 1998, 2001), such that $${\kappa}^{\text{line}}={f}^{*}{\kappa}_{0}{\displaystyle \sum _{\text{i}}{w}_{v,\text{i}}}{q}_{\text{i}}\left(\frac{1{e}^{{q}_{\text{i}}\tau}}{{q}_{\text{i}}\tau}\right)={f}^{*}{\kappa}_{0}M(\tau ),$$(11)
where q_{i} is the linestrength, w_{ν,i} the illumination function (here as defined by Eqs. (11) and (12) in Poniatowski et al. 2022) and the index i runs over all lines in the database. Then we have: $$\tau ={\kappa}_{0}c\rho {\left\frac{\text{d}{v}_{r}}{\text{d}r}\right}^{1}$$(12)
which is a characteristic Sobolev optical depth, with κ_{0} as the normalisation constant. In contrast to the previous Wolf–Rayet simulations of Moens et al. (2022b) we here also introduce the socalled finite disk correction factor, ƒ^{*} (see below).
It is important to realize that while the sum in Eq. (10) typically reproduces the correct optically thick and thin limits, it is only an approximation in between. If the medium becomes very dense or the velocity gradient approaches zero or both, τ (∝ ρ/dv_{r}/dr) becomes very large and κ^{line} → 0, so that κ → κ^{Ross}. On the other hand, in the limit that all lines would be optically thin, the sum in Eq. (11) yields κ^{line}(0 < q_{i}τ ≪ 1) → ${\displaystyle {\sum}_{i}}{w}_{v,i}}{q}_{i}\equiv \overline{Q$, where we have introduced Gayley (1995)’s $\overline{Q}$ (also see below). Following Poniatowski et al. (2022) this becomes $\overline{Q}={\displaystyle {\sum}_{i}{B}_{v,i}}{\alpha}_{l,i}/(B\rho )$, where α_{l} is the frequencyintegrated line extinction coefficient. Additionally, in such diluted parts of our simulations, we typically find that the OPAL opacities approach the Thomson opacity κ^{Ross} → κ^{Th}. That is, in this limit κ approaches an appropriate Planck mean for an ensemble of (intrinsically nonoverlapping) lines and a Thomson scattering continuum. In an ‘intermediate’ region, however, lines could in principle be ‘doublecounted’ as they may contribute in both κ^{line} and κ^{Ross} (see also Fig. 2). This is potentially a significant weakness in our current opacity approximation, which should be further investigated in future work. On the other hand, we indeed find good agreement between average massloss rates of the simulations presented here and those computed for O stars by means of full frequencydependent comovingframe radiative transfer (in 1D and assuming a steadystate), suggesting that the issue may not be critical, at least not for the parameter range explored in this paper (see further discussion in Sect. 4.3). Moreover, as outlined in Appendix A our approximation is also equivalent to the modification of Sobolevbased ‘expansion opacity’ formulae suggested by Castor (2004) based on calculations of Fe III lines.
The results from the complete lineensemble summation are fitted using (Gayley 1995): $${\kappa}^{\text{line}}={f}^{*}{\kappa}_{0}\frac{\overline{Q}}{1\alpha}\frac{\left({\left(1{Q}_{0}\tau \right)}^{1\alpha}1\right)}{{Q}_{0}\tau}={f}^{*}{\kappa}_{0}{M}_{\text{G}}(\tau )\text{.}$$(13)
Here, we choose the normalisation κ_{0} = 0.34 cm^{2} g^{−1}, which reflects the Thomson scattering opacity in a fully ionised plasma with a chemical composition similar to that of the Sun. Also, $\overline{Q}$, Q_{0}, and α are the line force fit parameters; $\overline{Q}$ sets the maximum value of the line opacity in the limit where all lines are optically thin (see above), whereas Q_{0} and α give an effective maximum line strength and powerlaw index for a mixture of optically thick and thin lines, respectively. For this work, we have computed a large table of κ^{line} values appropriate to the conditions for the O stars under consideration. For each density and temperature pair, we use the three line force fit parameters to fit M_{G}(τ) to the integrated M(τ) calculated directly from our atomic database, namely, we obtain and tabulate α(ρ, T), $\overline{Q}(\rho ,T),$, and Q_{0}(ρ, T); an example is given in Fig. 1. These fitted and tabulated parameters are then used within the timedependent simulations, providing varying and locally consistent values for κ^{line}. We note that since the line force parameters thus vary in space and time depending on the ionisation state of the model, there is no need to introduce further fitparameters (such as Abbott’s ‘δ', see discussion in Poniatowski et al. 2022). The technique is similar to how OPAL opacity tables are constructed from precalculations of the Rosseland mean opacity. Similarly to the OPAL tables, κ^{line} can be calculated for various chemical compositions. To obtain relative abundances we here use the default scale (Grevesse & Noels 1993) from the OPAL tables, thereby, hydrogen massfraction X = 0.70, helium massfraction Y = 1 − Z_{⊙} = 0.28 for metallicity Z_{⊙} = 0.02. We note that this differs slightly from the oftenused chemical abundance scale by Asplund et al. (2009; with Z_{⊙} = 0.013). Figure 2 displays ‘average’ opacity curves for one of the 2D O star models presented in the next section. The figure illustrates a characteristic opacity behaviour that is dominated by the Rosseland mean in deep optically thick layers, and by the lineopacity component in the wind outflow above the stellar surface. To a large extent, the diminishing contribution of linedriving to the total opacity in deep atmospheric layers is simply due to the inverse dependence with density seen in Eq. (13).
As shown above, ƒ^{*} is the standard finitedisk correction factor to the line force (Pauldrach et al. 1986; Friend & Abbott 1986), as given by, for instance, Eq. (29) in Poniatowski et al. (2022), which corrects the radial streaming expression by taking into account nonradial light rays and velocity gradients. It is here numerically implemented as in Owocki & Puls (1996). Above the stellar surface, this produces a varying ƒ^{*} depending on the local wind conditions. In the deeper regions of the atmosphere, we set a constant ƒ^{*} = 1/(1 + α_{★}) (the limiting value at the stellar surface), where α_{★} is the line force parameter α at the photosphere R_{★} (see definition below). Since R_{★} is not constant in our models (instead it varies with time depending on the local conditions in the atmosphere), in practice, we let our models relax until the overall variation of the stellar radius is small and then we use a constant R_{★} in the evaluation of ƒ^{*} for the remainder of the simulation. As can also be seen in Fig. 2, the values of κ^{line} are mostly small as compared to κ^{Ross} in the layers below R_{★}. A fixed R_{★} for the computation of ƒ^{*} is nevertheless used here in order to avoid additional initiation of structures due to small variations of ƒ^{*} at the stellar surface, as our implementation is approximate anyway and used here primarily to avoid issues related to ‘overloaded’ linedriven wind solutions in the radial streaming case. Indeed, strictly speaking, ƒ^{*} as implemented here is valid only for spherically symmetric systems; a full multidimensional evaluation would require costly numerical integration across the stellar disk (see, e.g. Petrenz & Puls 2000; Kee et al. 2016). Notwithstanding this caveat, this 1D finitedisk correction factor has been quite successfully applied in a range of other nonspherical linedriven systems like magnetically channeled O star winds (ud Doula & Owocki 2002). Nevertheless, this approximate 1D treatment should be improved upon in future work.
Fig. 1 Line force multiplier M(τ) as a function of the characteristic Sobolev optical depth τ. The blackdashed line is plotted by calculating M(τ) at temperature 30 kK and density 10^{−10} g cm^{−3}, whereas the red line is calculated by fitting M(τ) using M_{G}(τ) in Eq. (13). The best fit for the line force parameters for this particular temperature and density are $\overline{Q}=2278$, Q_{0} = 1418, and α = 0.59. 
Fig. 2 Averaged opacity 〈κ〉 as a function of scaled radius x = 1 − R_{★}/r, with R_{★} as the 2D averaged optical photosphere (reddashed line). In the figure, the brown curve is the Rosseland mean opacity and the purple curve is the finitedisk corrected lineopacity. The black dashdotted curve is the total opacity using the hybrid opacity scheme. 
2.3 Initial conditions
To compute the initial conditions for our timedependent simulations, we used a procedure similar to that outlined in Sect. 2 of SantolayaRey et al. (1997; 1st version of the FASTWIND unified model atmosphere code). This method assumes a spherically symmetric, stationary, and analytic wind structure described by a radial ‘β’ velocity law v_{r}(r) = v_{∞}(1 − br_{T}/r)^{β} and mass density ρ = Ṁ/(4πr^{2}v_{r}) on top of an atmosphere in (quasi)hydrostatic and (here) local thermal equilibrium. The two components of the atmosphere are smoothly connected at a transition velocity υ_{T} = ƒc_{iso,gas}(T_{eff}), where ${c}_{\text{iso,gas}}=\sqrt{\frac{{p}_{\text{g}}}{\rho}}$ is the isothermal gas sound speed, the parameter ƒ ≈ 0.1–1, and b is derived from v_{T}; this analytic wind structure provides an intermediate boundary condition for inward integration of the static gas and radiative momentum equations on a suitable grid of column mass m. For a given set of input parameters L_{★}, M_{★}, R_{★}, elemental abundances, Ṁ, v_{∞}, β (=1 throughout this work) then, this procedure is iterated until a unified structure that fulfills $${R}_{\star}\equiv r\left({\tau}_{\text{R}}=2/3\right)$$(14)
is found, where τ_{R} is the spherically modified Rosseland optical depth (see below). In practice, for opacity, κ, and gas pressure, p_{g}, we solve numerically the following coupled differential equations in the layers below the intermediate boundary given by r_{T}: $$\frac{\text{d}{p}_{\text{g}}}{\text{d}m}={g}_{\star}(1\text{\Gamma}){\left(\frac{{R}_{\star}}{r}\right)}^{2},$$(15) $$\frac{\text{d}T}{\text{d}m}=\frac{3\kappa {T}_{\text{eff}}^{4}}{16{T}^{3}}{\left(\frac{{R}_{\star}}{r}\right)}^{2},$$(16) $$\frac{\text{d}{\tau}_{\text{R}}}{\text{d}m}=\kappa {\left(\frac{{R}_{\star}}{r}\right)}^{2},$$(17) $$\frac{\text{d}r}{\text{d}m}=\frac{1}{\rho},$$(18)
along with the definitions and relations $$\sigma {T}_{\text{eff}}^{4}\equiv {F}_{\star}\equiv \frac{{L}_{\star}}{4\pi {R}_{\star}^{2}},$$(19) $${g}_{\star}\equiv \frac{G{M}_{\star}}{{R}_{\star}^{2}},$$(20) $$\Gamma \equiv \frac{{g}_{r}}{{g}_{g}}=\frac{{L}_{\star}\kappa}{4\pi G{M}_{\star}c},$$(21) $${p}_{\text{g}}=\frac{\rho {k}_{\text{b}}T}{\mu {m}_{\text{H}}},$$(22)
where g_{r} is the radial component of the radiation force as seen in Eq. (7). Also, for simplicity, we assume a constant µ = 0.61 (set by assuming fully ionised hydrogen and helium). Above, it is further implicitly assumed that the gas and radiation temperatures are equal at T. We also note that since here dP_{rad}/dm = g_{r}/ρ = κF/c = κL_{★}/(4πr^{2}c), Eq. (15) above may equivalently be formulated as a hydrostatic equation for the total pressure P_{tot} = P_{rad} + p_{g} of the radiating gas. To classify models, it is further useful to introduce the classical Eddington parameter, Γ_{e}, and luminosity, L_{edd}: $${\Gamma}_{\text{e}}\equiv \frac{{L}_{\star}{\kappa}_{\text{e}}}{4\pi G{M}_{\star}c}\equiv \frac{{L}_{\star}}{{L}_{\text{edd}}},$$(23)
where κ_{e} is the electron scattering opacity for a fully ionised medium of a certain chemical composition.
Two differences between our procedure here to compute initial conditions and the method outlined in SantolayaRey et al. (1997) are that we (i) do not account for any variations in ‘Hopfparameters’, but instead assume the Eddington approximation 3P_{rad} = E_{rad} = a_{r}T^{4} is valid throughout the atmosphere and (ii) use opacities κ = κ^{Ross} from the OPAL database (instead of fitting a Kramer’s like law) for a given set of input chemical abundances. The latter point allows us to probe much deeper atmospheric layers than FASTWIND so that we are able to consider the critically important regions around the socalled ‘ironopacity peak’ at T ~ 2 × 10^{5} K; specifically, while a typical FASTWIND model only extends down to column masses m ~ 100 we here terminate the downward integration only when a high temperature T ~ 4–5 × 10^{5} K has been reached (which in the O star regime occurs at m ≫ 100).
The black curve in Fig. 3 indicates the density structure for a template ‘early’ O star model (with small initial perturbations around the convective instability region, see Sect. 2.5). Additionally, the cyan curve corresponds to the FASTWIND model computed for the same set of input parameters. This model indeed agrees very well with our simplified 1D prescription, however the FASTWIND calculation only extends down to layers with ~100 kK, as explained above^{2}. The input parameters T_{eff}, g_{★} and Ṁ for the FASTWIND model are the same as those listed in Table 1 (O4 model). We note that the table excludes the terminal wind speed v_{∞} which for the 1D model is set to 2100 km s^{−1}, made to correspond to an extrapolation of the results of the 2D ‘average’ structure, and which is also typical for observed stars at this temperature range (Hawcroft et al. 2024).
From Fig. 3 (and later in Fig. 13), we can further see that in the region around the abovementioned ironopacity peak, there is a clear densityinversion. This occurs because in this region of a high Rosseland mean opacity, Γ > 1, such that the gradient in the hydrostatic Eq. (15) changes sign (see also e.g. Gräfener et al. 2012; Jiang et al. 2015; Köhler et al. 2015; Sanyal et al. 2015). To overcome this region of enhanced opacity the radiative stellar envelope reacts by expanding significantly, similar to what has been observed in envelopemodels of more evolved massive stars (Petrovic et al. 2006; Gräfener et al. 2012; Grassitelli et al. 2018; Poniatowski et al. 2021). Moreover, the transition region to the analytic wind outflow is also clearly visible in the figure, in particular through the characteristic change from an essentially exponential atmosphere to a more modest ~r^{−2} dropoff in density.
Fig. 3 Initial 1D gas density (in g cm^{−3}) structure as a function of modified radius x_{0} = 1 − R_{0}/r, with R_{0} as the lower boundary radius. The cyan curve is corresponding FASTWIND density structure as described in Sect. 2.3. The reddashed vertical line is the optical photosphere R_{★}. The shadedpurple region is the instability region where we provide the initial perturbations (see zoomedin part). 
2.4 Numerical grid
The lower boundary of our simulation is situated inside the stellar envelope, at a fixed R_{0} given by the initial condition calculation. The outer boundary extends well into the supersonic regions of the outflow, we set it here to r = 3R_{0}. This allows us to take a simultaneous and ‘unified’ approach to study the deep and optically thick layers of the atmosphere, the wind launching region, and the supersonic outflow up until a point where a majority of the gas particles have reached escape velocities.
The simulations that we present in this paper have been run on a Cartesian mesh, with four levels of refinement. Although MPIAMRVAC is capable of adaptive mesh refinement (AMR), in our simulations, we use a fixedradius refinement since otherwise, we run into methodological problems with smallscale ‘clumped’ regions in particularly the wind outflow. At our lowest resolution (Level 1), the numerical domain is 32 cells in the xdirection and 256 cells in the ydirection covering 0.2R_{0} and 3R_{0}, respectively. The ydirection is here identified with the radial coordinate r, applying the spherical correctionterms discussed in Moens et al. (2022b,a), and the xdirection is thus the lateral coordinate. Three times doubled to the highest resolution means we have 256 and 2048 cells in the lateral and radial directions, respectively. The highest level resolution (Level 4) cell covers 9.76 × 10^{−4} R_{0} in length. This highest resolution is set from the lower boundary R_{0} of the simulation up to the stellar surface R_{★} of the 1D model used as an initial condition.
Such a refinement is very useful for achieving proper resolution in the deeper optically thick atmospheric regions without net radial outflow while still keeping the total number of cells in the outer supersonic outflow regions computationally viable. To this end, we ensured that we resolve a few pressure scale heights in our simulations, especially in the deeper regions up to ~R_{★}. For a radiating gas in hydrostatic equilibrium, this pressure scale height is: $${H}_{\text{p}}^{\text{tot}}=\frac{{c}_{\text{iso,gas}}^{2}}{g{\beta}_{\text{p}}},$$(24)
where β = p_{g}/P_{tot} = p_{g}/ (p_{g} + P_{rad}). Figure 4 confirms a posteriori that we do indeed resolve several total pressure scale heights in layers deeper than and around R_{★}.
Fundamental parameters for the 〈2D〉 O star models studied in this paper.
Fig. 4 Number of total pressure height scale ${H}_{\text{p}}^{\text{tot}}$ resolved per cell in our simulation. Level 4 corresponds to the highest level of resolution in our models. Levels 3–1 are not shown since it is zoomed in to indicate the deep layers and the photosphere. See text for details. 
2.5 Perturbations
To initiate structure formation in our simulations, the initial conditions are perturbed in the predicted subsurface convective region. This region may be identified following the linear analysis in Blaes & Socrates (2003), where instabilities in an optically thick radiationpressure dominated gas were studied. The instability criterion for convection then becomes: $$\frac{{\gamma}_{\text{g}}{p}_{\text{g}}}{4\left({\gamma}_{\text{g}}1\right){E}_{\text{rad}}}{N}_{\text{g}}^{2}+\frac{1}{3}{N}_{\text{r}}^{2}0,$$(25)
where $$\begin{array}{l}{N}_{\text{g}}^{2}:={f}_{\text{grav}}\cdot \left(\frac{1}{\rho {c}_{\text{s},\text{gas}}^{2}}\nabla {p}_{\text{g}}\nabla \mathrm{ln}\rho \right),\\ {N}_{\text{r}}^{2}:={f}_{\text{grav}}\cdot \left(\frac{1}{3\rho {c}_{\text{s},\text{rad}}^{2}}\nabla {E}_{\text{rad}}\nabla \mathrm{ln}\rho \right).\end{array}$$
These are the BruntVaisala frequencies for the gas and radiation parts, respectively. This instability criterion may be interpreted as a modified Schwarzschild (1958) criterion for media with significant radiation pressure. c_{s,rad} and c_{s,gas} are the adiabatic radiation and gas sound speeds defined below in Sect. 2.7. The predicted convective instability region for the 1D initial condition model of a prototypical O star is shown in Fig. 3 as purpleshaded parts.
In this convective region, we add small perturbations (max 20%) to the initial density profile (see the zoomedin part in Fig. 3). Such perturbations have a plane wave dependence on r, namely, $$\delta \rho ~A\mathrm{exp}\{ik\cdot r\},$$(26)
with A a constant set here to 0.2ρ and k the wavevector of the perturbation. The final perturbations are a superposition of several wave modes with different wave numbers k, and we assume the perturbations have a Gaussian distribution in Fourier space. We consider a range of wavenumbers k limited by the total pressure scale height and the cell size of our simulations, namely, $${k}_{x}\in \left[\frac{2\pi}{{H}_{\text{p}}^{\text{tot}}},\frac{2\pi}{{L}_{\text{cell},x}}\right],{k}_{y}\in \left[\frac{2\pi}{{H}_{\text{p}}^{\text{tot}}},\frac{2\pi}{{L}_{\text{cell},y}}\right].$$(27)
Since we follow Blaes & Socrates (2003), where results are limited to the short wavelength regime, we set the scale height as a lower limit for the wavenumbers (and therefore also the wavelength). The number of wave modes N(k) with a certain wavenumber k is then a linear combination of Gaussians centered around (k_{x,c}, k_{y,c}), (k_{x,c}, − k_{y,c}), (−k_{x,c}, k_{y,c}) and (k_{x,c}, k_{y,c}), respectively, with $${k}_{x,c}=\left({k}_{x,\mathrm{min}}+{k}_{x,\mathrm{max}}\right)/2,{k}_{y,c}=\left({k}_{y,\mathrm{min}}+{k}_{y,\mathrm{max}}\right)/2.$$(28)
Finally, perturbations are calculated by Fourier transforming N(k), so that $$\delta \rho =0.2\rho {\displaystyle \int \mathrm{exp}}\{ik\cdot r\}N(k)\text{d}k$$(29)
Additionally, in this region, we add perturbations to the lateral velocity which are sinusoidal in r and y.
The analysis by Blaes & Socrates (2003) further shows that acoustic waves can potentially be unstable around opacity bumps in a massive stellar envelope (see also Owocki 2015). As in Jiang et al. (2015), we have focused here on the convective instability, however, we note that it is possible that such acoustic waves, which essentially are local versions of ‘strange mode’ instabilities (Glatzel 1994), could also be driven unstable in the present O star simulations.
2.6 Boundary conditions
In the lateral direction, the boundaries are periodic for all conserved quantities solved in Eqs. (1)–(3), (5). The mass density at the first ghost cell at the lower boundary is fixed at the value given by the initial conditions and is then extrapolated to the remaining ghost cells assuming hydrostatic equilibrium. For the gas momentum, we let it vary, extrapolating from the first active cell to the ghost cells assuming mass conservation. We set the radial component of the radiation energy E_{rad} by a fixed input radiative luminosity at the bottom boundary, using a finite difference formulation of Eq. (6) to extract E_{rad} from its gradient.
The gas energy e_{g} is then set by assuming T_{rad} = T_{g} at the lower boundary, using the ideal gas law and Eq. (4).
In the case of the outer boundary, the outflow is highly supersonic, and density, momentum, and gas energy are linearly extrapolated. For the radiation energy E_{rad}, we analytically solve Eq. (6) from the outer boundary radius r_{o} to r → ∞ by assuming constant opacity and fluxlimiter, and a diffusive flux and mass density that varies according to r^{−2} outside r_{o}, yielding ${E}_{\text{rad},\text{o}}=\frac{{F}_{\text{diff,o}}{r}_{\text{o}}}{3{D}_{\text{o}}}$ where subscripts o in the parameters signify that values are taken at the radial outer boundary.
Fig. 5 Mass loss rates (top) and effective temperature (bottom) for the O4 model as varying with time (shown here in 10^{3} seconds). The black curve displays the laterally averaged quantity varying across time written as 〈2D〉_{L}, and the blackdashed line is the lateraltemporal averaged quantity, here 〈2D〉. 
2.7 Timescales
The dynamics of the deep, optically thick atmosphere and the outflowing wind may in general be governed by different timescales. For the deep atmosphere without a net radial outflow, we may choose a characteristic dynamical timescale, t_{d}, as: $${t}_{\text{d},\text{a}}=\frac{{H}_{\text{p}}^{\text{tot}}}{{c}_{\text{s}}},$$(30)
where c_{s} is the total adiabatic sound speed in the radiating gas as such, $${c}_{\text{s}}^{2}={c}_{\text{s},\text{rad}}^{2}+{c}_{\text{s},\text{gas}}^{2},$$(31)
for adiabatic radiation and gas sound speeds of $${c}_{\text{s},\text{gas}}^{2}=\frac{5{p}_{\text{g}}}{3\rho},$$(32) $${c}_{s,\text{rad}}^{2}=\frac{4{P}_{\text{rad}}}{3\rho}.$$(33)
On the other hand, the dynamics of the outflowing wind, starting around an ‘average’ sonic point R_{sonic} ≃ r [v_{r} > c_{iso,gas}], which in our simulations vary in space and time depending on the local conditions, is typically better characterised by, $${t}_{\text{d},\text{w}}=\frac{{R}_{\star}}{v},$$(34)
with v a characteristic wind speed; here such a typical wind speed simply has been chosen from inferences of O star winds and some visual inspection of our simulations. For the O star regime considered in this paper, typical numbers are v ~ 1000 km s^{−1}, R_{★} ~ 17 R_{⊙}, ${H}_{\text{p}}^{\text{tot}}~0.03{R}_{\star}$, and c_{s} ~ 150 km s^{−1}, yielding t_{d,a} ~ 1000 s and t_{d,w} ~ 10 000 s.
Finally, thermal relaxation of the deeper atmospheric layers by radiative diffusion proceeds on a timescale related to the total thermal energy content per unit area divided by the total energy flux. Following Freytag et al. (2012) this may be approximated for the case here of radiationdominated O star atmospheres by ${t}_{\text{th}}~{E}_{\text{rad}}{H}_{\text{p}}^{\text{tot}}/{F}_{\text{diff}}~{10}^{3.\dots 4}\hspace{0.17em}\text{s}$. Alternatively, a thermal timescale may be estimated for the entire atmosphere following Grassitelli et al. (2016); Moens et al. (2022b), t_{th} = GM_{★} M_{atm}/ (R_{★}L_{★}), where M_{atm} represents the mass of the stellar atmosphere. Integrating over the average radial density profile of our simulations from the ironbump to the outer boundary, this gives typical numbers ≈ 40 000 s ≈ 4t_{d,w}. This estimate yields higher characteristic numbers than t_{d,w}, and in practice, we thus also examine whether our simulations have reached an approximate energetic steadystate by inspecting the temporal behaviour in total luminosity, see Sect. 3.2.
3 Simulation results
3.1 Classification of models
Generally, O stars are classified according to their spectral type, which is related to the effective temperature, T_{eff}, as defined at the photospheric radius, R_{★} (see Eqs. (14) and (19)). In our 2D simulations, however, it is the mass, radius, and radiative luminosity at the lower boundary that are specified as input parameters (see Sect. 2.6). As such, the photospheric radii and effective temperatures of the simulations are emergent properties that typically vary both in time and space, and in general, they do not have the same values as the corresponding ones of the 1D initial conditions (see detailed discussion in the next section). As such, we calculate lateral and temporal averages of fundamental stellar parameters and classify our models according to them; that is, hereafter, 〈X〉 for quantity X means a lateral and temporal average taken at some radial cell. In the case of 〈L_{★}〉 and 〈T_{eff}〉 the radial cell corresponds to R_{★} = 〈r(τ_{Ross} = 2/3)〉, whereas for 〈Ṁ〉 we take an average of a few radial cells close to the outer boundary of our simulation (where gas particles have become unbound). To ensure that transients stemming from the initial conditions do not influence our results we make sure that such averages are taken several dynamical time steps after our models have dynamically relaxed. Additionally, we also make sure that we take averages only after turning off a varying R_{★} when computing the finite disk correction factor (see previous section).
That our models are dynamically well relaxed is also evident from Fig. 5, showing that for advanced simulationtimes timefluctuations in 〈Ṁ〉 and 〈T_{eff}〉 are rather stochastic around the mean value. In the context of this Fig. 5, it is important to recall that we are only modeling 0.2R_{0} of the star in one of the two lateral directions, namely, only a fraction of the overall stellar and wind volume. As such, the timevariance displayed in this figure will notably not be representative of the predicted timevariability for observations of these global stellar and wind parameters. Specifically, 〈T_{eff}〉 is obtained by computing the diffusive radiative flux at the average photosphere 〈R_{★}〉. This allows us to organise the first O star simulations presented here according to their 〈T_{eff}〉, and thereby to also name them according to their approximate ‘spectral type’. Mainly, this represents a convenient organisation of our models and, thus, it does not necessarily reflect the ‘true’ spectral type the modeled stars would be classified as according to their spectral appearance. Table 1 provides averaged fundamental stellar parameters of our models, along with their averaged emergent mass loss rates.
Fig. 6 Colour maps of relative radiation temperature (top panel), relative density (second panel), and the radial velocity (third panel) and Γ (bottom panel), for several snapshots at different times (shown on top in 10^{3} seconds) at the beginning of our simulation for the O4 model. The figures are zoomed in to show the inner regions of our simulations. In the lateral direction, the figure extends to 0.2 R_{0} and in the radial direction to 1.5 R_{0}, with R_{0} as the lower boundary radius of our simulation. 
3.2 Structure formation, relaxation, and initial conditions
In Figs. 6 and 7, we show snapshots from the temporal evolution of our prototypical ‘O4’ simulation during the initial and dynamically relaxed phase. As can be seen in Fig. 6, structure starts developing within the perturbed ironopacity peak region, which also moves spatially inward in the simulation because of the readjusting temperature and density scales. Simultaneously, a selfconsistent linedriven wind outflow is initiated from the regions slightly above the optical photosphere at around 〈R_{sonic}〉. Since the stellar flux is not constant in the wind initiation region, this outflow also starts to develop structure, leading to a complex interplay between structure formation deep in the atmosphere and in the overlying wind. After about ~10 t_{d.w} (or, ~90 ks) the simulation starts to reach a more dynamically relaxed state^{3}.
As mentioned in the previous section, thermal relaxation of the deeper atmospheric layers may occur on different characteristic time scales. We examine this by taking lateral timeaverages of the total model luminosity 〈L_{tot}〉_{L} along the radius, where L_{tot} is obtained from the stationary radial energy equation including full contributions from diffusive radiation, radiative and gas enthalpy (convection), as well as kinetic and gravitational energies (see Eq. (35) and the corresponding discussion below). If timeaverages of 〈L_{tot}〉 are sufficiently constant over radius, we may consider the simulation to be energetically relaxed (see also Goldberg et al. 2022). The colored curves in Fig. 8 show laterally averaged total luminosities in the O4 simulation, for a timeaverage of 4 snapshots within ~t_{d,w}. By contrast, the black solid line shows a much longer average over the whole time range for which the colored lines were constructed, specifically for ~100t_{d,w}. As can be seen in the plot, the long timeaverage shows a fairly constant luminosity in radius, whereas the shorter averages display variations on the order of a few tens percent. Overall, this suggests that, in a statistical sense, our models are also energetically well relaxed.
After the initial adjustment period, Fig. 7 next displays timesnapshots from more advanced and dynamically relaxed stages of the O4 simulation. This shows how the Γ < 1 subEddington layers close to the lower boundary have become quasistable, characterised by low variation in temperature and velocity; the radial velocity dispersion is here about ~1 km s^{−1} (see Fig. 11) and the temperature dispersion is ~1.5%. As the temperature then decreases outwards from about 450 kK at the lower boundary, the opacity increases and creates local superEddington regions (below the photosphere) with Γ > 1. In accordance with previous simulations of this opacitydriven unstable zone (Jiang et al. 2015), the net result is a highly structured and turbulent subsurface atmosphere characterised by temperature fluctuations, strong density contrasts, and large turbulent velocities that ‘overshoot’ into the overlying cooler surface layers with lower opacities. In contrast to previous work, however, in the simulations presented here the turbulent gas stemming from the deeper atmosphere naturally interacts with the linedriven layers around and above the photosphere. This complex interplay creates a situation wherein the linedriven wind acts like a suction mechanism, preventing some gas parcels that otherwise would have turned over from falling back into the deeper stellar envelope. Simultaneously, the turbulence arising from the deeper atmosphere introduces a natural variation in flux and lineopacity in the wind launching region. This leads to further structure formation and prevents the development of a stationary linedriven outflow as seen in 1D O star simulations with a similar Sobolevbased κ^{lıne} as here but assuming a static stellar surface and constant stellar luminosity (Poniatowski et al. 2022). The overall trend in behaviour is quite similar in all three models considered in this work.
To examine how sensitive our models are to initial conditions, we rerun the O4 simulation but starting now instead from an essentially hydrostatic envelope (achieved in practice by assuming a negligible massloss rate ~10^{−10} M_{⊙} yr^{−1} and v _{∞}~ 10 km s^{−1} as boundary conditions in the 1D procedure outlined in the Sect. 2.3); that is, we start from a hydrostatic initial condition with effectively no wind to investigate if a linedriven outflow is still launched. Figure 9 demonstrates the evolution of suitable running averages for radial velocity as a function of time, taken specifically for 10 snapshots over a few t_{w,d} starting from the beginning of the simulation. Additionally, the figure also displays the initial conditions and a longer timeaverage taken after the wind has dynamically relaxed. The top panel shows the new ‘O4^{*}’ model starting from the effectively nowind condition (note the nearzero velocity curve for the initial condition); the bottom panel displays our standard O4 model. As discussed elsewhere in the paper, for this O4 model the adjustment of the atmosphere from the initial conditions involves inward movement of the photosphere. This leads to a similar adjustment of the linedriven wind, so that while the final average wind outflow is initiated from approximately the same optical depth layers as the initial condition this now occurs closer to the lower boundary radius. Without a significant initial velocity gradient, however, linedriving is at first quite inefficient in the O4^{*} simulation and the whole upper atmosphere begins to fall inwards. Since this then gives rise to velocity gradients, linedriving starts to kick in and an outflow is subsequently developed during later snapshots. Comparing the final wind velocity curves of the O4^{*} and O4 models, the only significant difference in the end is that the O4^{*} model, as expected, takes much longer to develop a statistically relaxed outflow. This shows that the resulting linedriven winds in our O star simulations are quite robust against specific details in the initial wind conditions.
Fig. 7 Similar to Fig. 6, but now at times after the models have dynamically relaxed from their initial conditions. The blackdashed line is the averaged 2D optical photosphere 〈R_{★}〉. 
Fig. 8 Running averaged total luminosity 〈L_{tot}〉 as a function of scaled radius x = 1 − R_{★}/r. The thin lines represent lateral averages taken over a few time snapshots (with lighter colours representing earlier snapshots), whereas the black line represents a long timeaverage of the total luminosity. 
3.3 Model O4
We go on to examine in greater detail our prototypical model O4 for an early Osupergiant in the Galaxy. As discussed above, after the initial relaxation, the O4 simulation instigates a turbulent deeper atmosphere driven by the ironopacity peak, with linedriving taking over slightly above the optical photosphere to initiate an overall outflow. That is, although many gas parcels do indeed reach Γ ≳ 1 already around the subsurface ironopacity peak, the O stars considered in this paper are not able to drive a net radial wind outflow from these regions. This is a key difference between the O star models of this paper and the optically thick winds arising from the simulations of WRstars in Moens et al. (2022a).
In Fig. 10, we illustrate probability density profiles (analogous to normalised 2D histograms) for radial velocity, radiation temperature, gas density, and Γ for the O4 model (middle panel). The colour coding for the plots represents the likelihood of finding a cell with a specific value for the given quantity (v_{r}, ρ, etc.) at a given radial coordinate. Thereby, at any radius, values that are coloured yellow are more likely to appear in the simulation than values that are coloured blue. The overplotted red curve is then the laterally and time ‘averaged’ quantity. As evident from Fig. 10, at the ironopacity peak 〈Γ〉 indeed grazes unity before again decreasing as we move outwards in the atmosphere toward R_{★}. Since there is no netoutflow from these subsurface layers, this in turn produces a very shallow 〈ρ〉 profile in the region where 〈Γ〉 ~ 1 , however in contrast to the initial conditions (see Fig. 3) there is no sharp inversion present but rather a plateaulike feature. As 〈Γ〉 then decreases,〈ρ〉 starts to decrease more rapidly again until linedriving eventually becomes efficient enough to produce higher 〈Γ〉 (recall that essentially κ^{line} ~ ρ^{−α}, Eqs. (12), (13)). When this occurs, 〈Γ〉 quickly reaches values well above unity, thereby launching a linedriven supersonic outflow. This general behaviour is then also reflected in the 〈v_{r}〉 profile. Namely, around the ironopacity peak 〈v_{r}〉 first increases slightly to around ~20 km s^{−1}. The following decrease in 〈Γ〉 then produces a region with a netinfall of material before linedriving takes over above the average photosphere. As such, the net effect is that at R_{★} we interestingly enough observe a negative 〈v_{r}〉 40 km s^{−1} in our simulation, indicating that the average O star surface is actually slightly infalling.
This general situation further produces large velocity dispersions, as illustrated in Fig. 11. Focusing again on the O4 model, we see how the radial velocity dispersion 〈v_{dısp}〉 indeed is close to zero in the deepest layers of our simulations but then increases rapidly once instabilities connected to the ironopacity peak become effective. The moving gas particles in this region then shoot into the cooler upper atmosphere, interacting with the linedriven outflow launched from the turbulent surface layers.
For the O4model, 〈v_{dısp}〉 ~ 80–90 km s^{−1} around R_{★}, after which it increases even further in the wind. As further discussed below, this general behaviour is overall in good agreement with O star observations of photospheric optical absorption lines as well as UV resonance lines formed in the wind (Hawcroft et al. 2021). Let us also point out here that although the characteristic (sub)photospheric 〈v_{dısp}〉 values in our simulations are above the ‘gas’ sound speed, c_{s,gas}, they are still below the ‘total’ sound speed, c_{s}, of the radiating atmosphere (see Eqs. (31)–(33)).
Figure 12 illustrates the average relative density for each radial velocity at different radii for the O4 mode. Blue here represents highdensity clumps and yellow represents the underdense structures. Although the underdense structures reach their local escape velocities very fast and close to R_{★}, on average the wind only escapes much further out. This is because the underdense structures move at much higher velocities than their dense counterparts. As discussed above, linedriving is the main driver of the wind outflow here. However, since lineopacity and gas density essentially are inversely proportional to each other (see above), underdense structures experience much stronger linedriving than dense clumps, resulting in significantly higher velocities for the former. The upshot is a strong anticorrelation between density and velocity in the O star winds studied here; that is, dense clumps in the wind move significantly slower than the rarefied gas in between them (see also Moens et al. 2022a for a similar result for WRstars). By contrast, below R* the gas is quite dense, linedriving is not efficient, and κ thus does not have this direct inverse dependence on density. This then disrupts the above mechanism, resulting in a turbulent atmosphere without any such clear (anti) correlation between velocity and density.
Fig. 9 Running averages of the radial velocities for the O4^{*} model (top) and O4 model (bottom), taken over 10 snapshots starting from the beginning of the simulation, with lighter colours representing earlier times. The red dashdotted lines are the initial velocity profiles. The green dashdotted line is the final average for the O4^{*} model and the orange dashed line is for the O4 model. See text. 
Fig. 10 Probability density maps for different quantities, from top to bottom radial velocity, gas density, Eddington Γ, and radiation temperature for different models O8 (left), O4 (middle), and O2 (right). 40 snapshots have been used to create these probability density maps. The red curves represent the lateraltemporal averages. The vertical orangedashed line is 〈R_{★}〉, whereas the green dashdotted vertical lines in the top panel are 〈R_{sonic}〉. The cyan horizontal lines in the Γ panels are where its value is unity. The colours signify the probability of finding a quantity at a certain cell, with yellow having a higher probability and blue having a lower one. 
Fig. 11 2D average radial velocity (dashed lines) and its dispersion (solid lines) as a function of the scaled radius x = 1 − R_{★}/r. The colours correspond to the different models. 
3.4 Models O2 and O8
Here, we discuss two additional models, one (O2) with a higher luminosity and thus higher classical Γ_{e} = 〈L_{★}〉/L_{edd}, resulting in a higher 〈T_{eff}〉, and one (O8) with a lower mass and luminosity, thus resulting in a lower Γ_{e} and 〈T_{eff}〉. We computed these models by means of the same procedure; however, by varying the input parameters of the initial condition calculations, we can obtain different combinations of (lower boundary and stellar) radius, luminosity, and so on, for the simulations; the resulting 〈2D〉 fundamental parameters of these models are provided in Table 1.
Inspecting once again the probability density clouds of Fig. 10, we observe that all three models exhibit similar qualitative behaviour. The lower boundaries for all models are (quasi)stable, significant structures start developing around the ironopacity peak, and linedriven winds are initiated around and slightly above the stellar surface. However, due to its lower baseline classical Eddington ratio, 〈Γ〉 of the cooler O8 simulation on average stays slightly below unity throughout the ironopacity peak opacity region, leading to somewhat lower (but still significant) dispersion in the atmosphere. Vice versa, the hotter O2 simulation has a higher Eddington ratio, and as such 〈Γ〉 again grazes (even slightly exceeds) unity already in the subsurface layers. But since the atmosphere there still is too dense for efficient linedriving, 〈Γ〉 drops again once the ironopacity peak is overcome so that a deep subsurface netoutflow is not launched, mimicking the O4 model. We note, however, that for this very luminous model, it is somewhat ‘easier’ for linedriving to become efficient, meaning a lower lineopacity forcemultiplier boost is needed to overcome the effective gravity and launch an outflow. As such, in the O2 simulation, the average sonic point 〈R_{sonic}〉 lies at 〈τ_{Ross}〉 = 0.62, meaning the average wind is on the verge of becoming (marginally) optically thick, and also that the characteristic region with negative average velocity now is located slightly beneath the stellar surface. Again this demonstrates how the simulations presented here, and previously in Moens et al. (2022a), are able to very naturally capture the transition from lower luminosity stars with optically thin linedriven winds (O stars, hot subdwarfs) to higherluminosity objects with optically thick winds (‘slash’ or WNhstars, classical WRstars).
The variability in the subsurface layers of the three simulations also follows simple qualitative trends with increasing classical Eddington ratio, Γ_{e}; the dispersion in density, temperature, and velocity is lowest for the O8 simulation and highest for the O2 model. Specifically, Fig. 11 illustrates how 〈v_{dısp}〉 around the photosphere exceeds 100 km s^{−1} in the O2 model, is below 50 km s^{−1} in the O8 model, and (as discussed already above) lies in between for the O4 model.
Fig. 12 Average relative density for each radial velocity at different radii for the O4 model. 〈ρ_{r}〉 is the average density at every radial cell. Colours here represent relative density (ρ/〈ρ_{r}〉) with yellow displaying underdense regions and blue representing overdense clumps. The red dashdotted line indicates the local escape velocity, the black dashed vertical line is the average photosphere 〈R_{★}〉 and the black dashdotted vertical line is the average sonic point 〈R_{sonic}〉. The zoomedin part highlights the regions below the optical photosphere. 
4 Comparison to 1D models
The previous section demonstrated the development of a very turbulent, timedependent O star atmosphere coupled with a wind outflow. We go on to examine how the ‘average’ properties of these atmospheres may compare to standard 1D, stationary models.
Figure 13 compares the 〈2D〉 density and temperature structure of the O4 model to those of 1D models computed by means of the procedures outlined in Sect. 2.3. We note that (for visualisation purposes) the abscissa in this figure displays a scaled radial coordinate x0 and that the 〈2D〉 and 1D models here have been calibrated such that they have the same R_{0}. As we do not know in advance the emergent average parameters of the stellar wind, such as mass loss rates and terminal wind speeds, the outer part of the 1D initial conditions will not generally correspond to our 2D model. Hence, we reran the 1D model with the same stellar parameters (L_{★}, M_{★}, and R_{0}) but with the updated emergent wind properties, thus allowing us to make a fair comparison. The figure clearly shows how the 〈2D〉 model has significantly less envelope expansion than the 1D model, and thereby also a lower photospheric radius R_{★} and higher T_{eff}. Specifically, for the 〈2D〉 and 1D O4 models we obtain R_{★} = 16.7R_{⊙} and 18.7R_{⊙}; T_{eff} = 39.6 kK and 37.2 kK, respectively. Moreover, the sharp density inversion seen in the 1D initial conditions is no longer present in the 〈2D〉 density structure, and the slope of the density profile is significantly shallower in the 〈2D〉 model around the photosphere (i.e. the characteristic photospheric scale height is larger). This indeed suggests that processes typically not included in standard 1D unified model atmosphere with wind codes may be important for also the average structure of these key atmospheric layers.
Fig. 13 Initial 1D and resulting 〈2D〉 gas density [in g cm^{−3}] (top) and temperature [in K] (bottom) structure for the O4 model as a function of modified radius x = 1 − R_{0}/r, with R_{0} as the lower boundary radius. The ‘red dashdotted’ line is the 〈2D〉 photosphere and the ‘blue dashdotted’ is the 1D photosphere. The cyan dashdotted curve is produced using FASTWIND as described with specifications in Sect. 2.3. 
4.1 Convective energy transport and turbulent pressure
Modified 1D models. Spurred by the findings above and based on an interpretation of the structures observed in the deep atmosphere of our 2D simulations as primarily a result of turbulence due to the convective instability (see also discussion in Jiang et al. 2015), we have incorporated a simple treatment of convection and turbulence in our 1D code for computing initial conditions, following Kippenhahn et al. (2013) their Chapter 7^{4}. This treats energy transport within the standard mixinglength theory (MLT) framework for nonadiabatic convection, accounting for the effects of radiation pressure and cooling. The free input parameter, α_{MLT} = ℓ/H_{p}, is the standard MLT parameter for mixing length, ℓ, and ‘total’ pressure scale height, H_{p} = P/(dP/dr), and, for simplicity, we have applied the standard Schwarzchild criterion for when convective energy transport occurs (rather than the more rigorous criterion discussed in Sect. 2.5). Additionally, we have also included a turbulent pressure term ${P}_{\text{turb}}=\rho {v}_{\text{turb}}^{2}$ in the hydrostatic Eq. (15), where v_{turb} is an input turbulent velocity parameter, such that now dP/dm = GM_{*}/r^{2} for total pressure P = p_{g} + P_{rad} + P_{turb}. The connection of the subsurface layers to the overlying analytic wind outflow then follows the same procedure as before.
To make suitable comparisons, we first calibrated the new 1D models such that their T_{eff} and R_{★} now agree with the corresponding 〈2D〉 simulations. If this calibration is done without the addition of convective energy and turbulent momentum transport, the result is again that the 1D models have a much more inflated envelope, a density inversion around the ironopacity peak, and a steeper densityprofile around the photosphere.
First adding convective energy transport by setting α_{MLT} = 1.0, but still keeping v_{turb} = 0 km s^{−1}, gets rid of the densityinversion also in the 1D structure. However, around the photosphere (above the convectively unstable region) the scaleheight is still too small, producing the same steep densityprofile.
In the next step, we first add a constant v_{turb} ≃ 90 km s^{−1} throughout the atmosphere in our 1D model. This then indeed produces a shallower densityslope around the photosphere, but due to the constant turbulent velocity throughout the complete atmosphere, we (re)introduce the envelope expansion of the lowermost layers. Finally, then, we introduce a simple variation of v_{turb} such that at the photosphere we have a maximum v_{turb} ≃ 90 km s^{−1} which then gradually decreases below the photosphere to essentially 0 km s^{−1} in the lowermost layers. This produces a reasonable agreement between the 1D and 〈2D〉 results, as illustrated in Fig. 14. It is important to here note that we have not aimed for a perfect match of the 〈2D〉 and 1D structures. Rather, our goal has merely been to demonstrate that the overall effects seen in the subsurface layers of the 2D simulations may, in principle, be captured by these modifications of presentday 1D models. More quantitative comparisons and benchmarks will, however, require more indepth analysis of the detailed structures seen in a larger array of simulations. Moreover, we have here not touched upon how the rather complicated atmospheric 〈2D〉 average velocity field might be captured in such 1D stationary models (e.g. the negative average radial velocity around the photosphere and the subsequent rather complex wind initiation region).
Nonetheless, it is quite striking how much steeper the density and temperature slopes around the photosphere are in the FASTWINDlike 1D comparison models. We may understand this large difference by comparing the characteristic density scale height in a model accounting for turbulent pressure to that in a purely radiative model, ${H}_{\rho}^{\text{turb}}\approx {H}_{\rho}^{\text{rad}}\left(1+{v}_{\text{turb}}^{2}/{a}_{\text{g},\text{iso}}^{2}\right).$ For v_{turb} ~ 90 km s^{−1} and α_{g,iso} ~ 30 km s^{−1} (appropriate for the photospheric layers) then, we get thus ${H}_{\rho}^{\text{turb}}\approx 10{H}_{\rho}^{\text{rad}}$. That is to say, the density scale height is increased by a full order of magnitude, explaining the significantly different slopes seen in the models displayed in Figs. 13 and 14. Moreover, since the turbulent velocities in our 2D simulations decrease rapidly below the iron opacity peak, the effective density scale height there becomes significantly reduced, leading to less envelope expansion than in 1D comparison models with high v_{turb} also in the lowermost atmosphere. As further discussed in Sect. 5, these differences may have a significant effect on spectral line formation and thus spectroscopic determination of fundamental parameters in the O star regime.
2D model averages. The parameters introduced and values used above to (at least qualitatively) reproduce the average density and temperature profiles of the 〈2D〉 atmosphere may further be directly compared to and interpreted through the radiationhydrodynamic equations. Taking the stationary, radial components of the combined gas and radiation momentum and energy equations give (Mihalas & Mihalas 1984; Turner & Stone 2001; Goldberg et al. 2022; Jiang 2023) $$\frac{\partial}{\partial r}\left(\rho {v}_{r}^{2}+{p}_{\text{g}}\right)=\rho \frac{G{M}_{*}}{{r}^{2}}+\rho \frac{\kappa {F}_{\text{diff}}}{c},$$(35) $${v}_{r}\left({p}_{\text{g}}+{P}_{\text{rad}}+{e}_{\text{g}}+{E}_{\text{rad}}+\frac{\rho {v}_{r}^{2}}{2}\rho \frac{G{M}_{*}}{r}\right)+{F}_{\text{diff}}={F}_{\text{tot}}=\frac{{L}_{\text{tot}}}{4\pi {r}^{2}},$$(36)
where L_{tot} = const. is the total luminosity, and where we have assumed an isotropic velocity field in the momentum equation. We note further that for the case γ_{g} = 5/3 and 3P_{rad} = E_{rad}, which is a good assumption in the deep optically thick layers, the stationary energy equation above may be equivalently written as: $$\dot{M}\left({h}_{\text{g}}+{h}_{\text{r}}+\frac{{v}^{2}}{2}\frac{G{M}_{*}}{r}\right)+{L}_{\text{diff}}={L}_{\text{tot}},$$(37)
for gas and radiation enthalpies h_{g} = (5/2)p_{g}/ρ and h_{r} = 4P_{rad}/ρ, respectively, and diffusive luminosity L_{diff} = 4πr^{2} F_{diff} (e.g. Owocki et al. 2017). These stationary relations invite us to identify the following suitable lateral and timeaverages as corresponding to ‘turbulent pressure’ and ‘convective flux’ in our simulations: $$\overline{{P}_{\text{turb}}}=\langle \rho {v}_{r}^{2}\rangle ,$$(38) $$\overline{{F}_{\text{conv}}}=\langle {v}_{r}\left({p}_{\text{g}}+{P}_{\text{rad}}+{e}_{\text{g}}+{E}_{\text{rad}}\right)\rangle ,$$(39)
with associated turbulent and convective velocities $$\overline{{v}_{\text{turb}}}=\sqrt{\frac{\overline{{P}_{\text{turb}}}}{\langle \rho \rangle}},$$(40) $$\overline{{v}_{\text{conv}}}=\frac{\overline{{F}_{\text{conv}}}}{\langle {p}_{\text{g}}+{P}_{\text{rad}}+{e}_{\text{g}}+{E}_{\text{rad}}\rangle}.$$(41)
That is, the turbulent velocity is here simply identified as the densityweighted average root mean square (rms) velocity (see also Goldberg et al. 2022). As mentioned above, this expression for $\overline{{P}_{\text{turb}}}$ essentially assumes that the turbulent velocity field is isotropic in the relevant layers (specifically for our case that $\langle \rho {v}_{r}^{2}\rangle \approx \langle \rho {v}_{t}^{2}\rangle )$, or alternatively (as pointed out in Jiang 2023) that the scale height is much smaller than the radius; these conditions are generally true in the subsurface layers of our simulations that do not experience a net radial outflow. On the other hand, the appropriate convective velocity is rather identified from the energy equation, and as such in general $\overline{{v}_{\text{turb}}}\ne \overline{{v}_{\text{conv}}}\text{.}$. For the O star simulations here the additional kinetic and gravitational energy fluxes are small, and $\overline{{F}_{\text{conv}}}\text{}$ is now essentially what MLT tries to estimate (see also discussion in Jiang 2023). As can be seen in Fig. 15, the maximum amount of energy transported by $\overline{{F}_{\text{conv}}}\text{}$ in our O4 model is ~ 10% of the total flux. It is dominated by the radiative enthalpy component, with gas enthalpy providing only ~0.1% of the total flux. $\overline{{F}_{\text{conv}}}\text{}$ is further centered around the ironopacity peak regions, but with a soft upper boundary that reaches significantly cooler layers than the corresponding 1D model. For this 1D model with α_{MLT} = 1 (calibrated to approximately match the 〈2D〉 density and temperature profiles) the maximum energy transported via F_{conv} is ~20%. We note that to obtain similar energy transport by the convective flux in the 〈2D〉 and 1D O4 models one only needs α_{MLT} = 0.5. However, this is not enough to wash away the density inversion in our 1D model. This difference will require more analysis in a future work, but this effect may stem from the simple fact that our 1D model applies a simple Schwarschild criterion without ‘overshooting’ for the boundaries of the convective zone. As such, the region of convective transport is much narrower than in the 〈2D〉 model (see Fig. 15), which then may require somewhat more efficient transport to obtain a similar net effect.
For the O2 and O8 models, we find a similar behaviour of $\overline{{F}_{\text{conv}}}\text{}$, but with maximum levels of $\overline{{F}_{\text{conv}}}\text{}$ that reach ~11% and ~4% of the total flux from 2D respectively. Similar to the O4 model, gas enthalpy is quite insignificant in comparison to radiation enthalpy for both models. As evident from Fig. 16, for the O4 model, in the layers leading up to the photosphere, radiation and turbulent pressures dominate the total pressure balance. As we approach the deeper layers of the stellar envelope, the turbulent pressure becomes lower so that we have a clear radiationdominated atmosphere around the ironopacity peak. As the opacity then decreases again, the turbulent pressure declines rapidly and the gas pressure, p_{g}, starts to become somewhat more significant, approaching (though remaining below) the level of radiation pressure at the lower boundary. In the bottom panel of Fig. 16, we display the turbulent velocity $\overline{{v}_{\text{turb}}}$ for the O4 model. The qualitative behaviour of $\overline{{v}_{\text{turb}}}$is quite similar to the straight average v_{dısp} discussed in the previous section, with very similar values in the deep subsonic layers. Once we approach the upper atmosphere, we notice a difference in that the densityweighted $\overline{{v}_{\text{turb}}}$ is somewhat lower than v_{dısp}. This is most likely related to the influence of the linedriven wind and the densityvelocity anticorrelation discussed in the previous section. For the O4 model in the layers around the optical photosphere, we observe the turbulent velocities to be ~60–80 km s^{−1}. The general behaviour seen in the simulation is, thus, quite consistent with the modification of the pressure balance introduced above in the corresponding 1D model. Interestingly, $\overline{{v}_{\text{conv}}}$ only reaches a peak value (in the ironopacity peak region) of ~20 km s^{−1}, which is much lower than the turbulent velocity discussed above; effectively this means that the turbulent momentum transport in our simulations is much more efficient than the advective energy transport in the ironopacity peak region. Moreover, below the optical photosphere, although the gas sound speed c_{Sıgas} is around ~30 km s^{−1}, the total sound speed c_{s} is above 100 km s^{−1}. As such, both $\overline{{v}_{\text{turb}}}\text{and}\overline{{v}_{\text{conv}}}$ stay well below the local total sound speed in our simulations.
Broadly, we also observe the same trends for $\overline{{v}_{\text{turb}}}$ in the O2 and O8 models as for the O4 simulation, with $\overline{{v}_{\text{turb}}}$ at the photosphere ~100 km s^{−1} for the O2 simulation and ~30 km s^{−1} for the O8 model. In the case of the O8 simulation, p_{g} and P_{rad} are quite similar in the deeper subsurface layers, with pg even slightly above at the lower boundary. On the other hand, for the O2 the qualitative behaviour of p_{g} and P_{rad} are similar to the O4 model, although $\overline{{P}_{\text{turb}}}$ approaches P_{rad} much quicker than in the O4 model.
Fig. 14 〈2D〉 gas density (in g cm^{−3}; top panel) and temperature (in K; bottom panel) as a function of the scaled radius x = 1 − R_{★}/r for the O4 model. The 1D structure is calibrated to the 〈T_{eff}〉 and 〈R_{★}〉 obtained from the 2D model. The orange line displays the calibrated 1D model with an α_{MLT} = 1 and a turbulent velocity of 90 km s^{−1} at the photosphere and reduced gradually in the deeper regions close to the hydrostatic envelopes. The cyan dashdotted curve is produced using FASTWIND, using the specifications outlined in Table 1 and Sect. 2.3. 
Fig. 15 Convective flux to the total flux as a function of temperature for the O4 model. The black curve is the 〈2D〉, and the orange curve is the calibrated 1D model to match the 2D with α_{MLT} = 1. 
Fig. 16 2D average gas pressure, radiation pressure, and turbulent pressure (top panel) and turbulent velocity, dispersion of radial velocity, and radial velocity (bottom panel) for the O4 model as a function of the scaled radius x = 1 − R_{★}/r with the upper axis as radiation temperature in 10^{3} K. 
4.2 Porosity reduction of radiation force
Because of potential (anti)correlations between density, radiative diffuse flux, and opacity, the averaged radiation force can be different in multidimensional simulations as compared to 1D models based on their individual averages. That is, generally 〈ρκF_{diff}〉≠ 〈ρ〉〈κ〉〈F_{diff}〉 (Shaviv 1998, see also discussion in Jiang 2023). Figure 17 displays this difference for our O4 simulation, illustrating that while the quasisteady lowermost atmosphere shows no significant modification, the averaged radiative force becomes reduced in our 2D simulations when approaching the photosphere. This is in qualitative agreement with the analysis by Schultz et al. (2020; based on the simulations by Jiang et al. 2015, 2018), although the reduction here is quantitatively somewhat lower. We have not included this effect in the 1D comparison models above, also because in layers where the reduction is most prominent we follow the customary approach for 1D unified model atmospheres used for spectroscopy and apply an analytic βtype velocity law (see previous sections). The effect could be significant also for deeper quasistatic layers, however, and thereby modify the overall pressure gradient balance discussed above; this should be investigated in detail in future work when trying to further calibrate 1D atmospheric models (see also discussion in Schultz et al. 2020). Moreover, we note that these types of ‘porosity effects’ may be important for wind launching in other domains (Owocki et al. 2004) as well as for capturing spectrum synthesis effects stemming from the clumpy atmosphere and wind within stationary 1D models (see discussion in Sect. 5.3).
Fig. 17 Fluxporosity effect as a function of scaled radius for our O4 model. 
4.3 Mass loss rates
Additionally, we compared the averaged mass loss rates 〈Ṁ〉 (see Table 1) of our models to mass loss rates derived by dedicated 1D, stationary theoretical modeling (here Krticka & Kubát 2017; Björklund et al. 2021). Since we are using chemical abundances prescribed by Grevesse & Noels (1993) we applied a simple metallicity scaling 0.02/0.013 in the O star recipe by Björklund et al. (2021), to reflect the different baseline solar abundance scale. Similarly, for the Krticka & Kubát (2017) recipe, we also applied this metallicityscaling based on Björklund et al. (2021). Although this recalibration is not perfect (since changes in abundances of important driving ions do not necessarily follow the baseline metallicity scaling, see discussion in Sundqvist et al. 2019), it nonetheless captures the principal effect. Additionally, we also compared to the recent empirical study of luminous Galactic Osupergiants by Hawcroft et al. (2021). The overall good agreement between the average 〈Ṁ〉 computed here and the rates predicted by these (albeit 1D, stationary) more focused studies (see Fig. 18) provide further support that we indeed are able to capture the linedriving effect rather well within our new modeling framework (see also Poniatowski et al. 2022).
In this context, it is further important to recall that comparisons of ‘predicted’ 〈Ṁ〉 values cannot typically be carried out against standard 1D unified model atmosphere codes used for spectroscopic studies; this is because in these cases, Ṁ is most often (though see, e.g. Sander et al. 2017) treated as a free inputparameter accompanied by a parameterised analytic velocity field (see details previous sections). By contrast, in the simulations presented here 〈Ṁ〉 and the atmospheric velocity field are intrinsic and selfconsistent emergent properties of the models.
Fig. 18 Mass loss rates in M_{⊙} yr^{−1} as a function of solar luminosity. The ★markers are the mass loss rates calculated from the models considered in this work, whereas the dotdashed line using Krtička & Kubát (2017), and the dashed line is calculated using Björklund et al. (2021). The triangle markers are plotted from the empirical mass loss rate study by Hawcroft et al. (2021; their optically thick GA rates). 
5 Discussion
5.1 Convective energy transport
As shown in Sect. 4.1, radiation enthalpy is able to carry only a small part of the total flux through the subsurface iron opacity peak region of our model O stars. This means our simulations are generally located in the region of socalled ‘inefficient’ convection, which we may interpret through some simple scaling relations. For the radiationdominated atmospheres studied here, an upper limit for convective flux would scale as (e.g. Gräfener et al. 2012) F_{Conv} ~ C_{s}E_{rad}, so that F_{Conv}/F_{diff} ~ (c_{s}/c)(T/T_{eff,0})^{4} ~ (c_{s} /c)τ, where ${T}_{\text{eff},0}^{4}={L}_{\star}/\left(\sigma 4\pi {R}_{0}^{2}\right)$ is an effective (flux) temperature evaluated below the ironopacity peak at ~R0 and where the generic scaling τ ~ (T/T_{eff,0})^{4} has been applied further. Alternatively, we may view this scaling through the characteristic dynamical and radiative diffusion timescales for the region, ${t}_{\text{d},\text{a}}~{H}_{\text{p}}^{\text{tot}}/{c}_{\text{s}}$ and ${t}_{\text{diff}}~(\tau /c){H}_{\text{p}}^{\text{tot}}$, such that t_{diff}/t_{d,a} ~ (c_{s}/c)τ. That is, when the diffusion timescale is shorter than the dynamical timescale, gas particles quickly adjust to their surroundings such that energy transport by means of convection (advection) is no longer efficient for carrying a significant fraction of the total flux entering fr1om the quasistable radiative zone below. Setting t_{diff} = t_{d,a} then also gives us the ‘critical optical depth’ τ_{c} = c/c_{s} as defined and used in Jiang et al. (2015) to distinguish between efficient and nonefficient convective energy transport. Inspection gives τ on the order of hundreds for our simulations (see also corresponding figures displaying a mean optical depth scale), so that F_{conv}/F_{diff} ~ t_{diff}/t_{d,a} ~ 0.1 indeed is a quite good characteristic scaling number, as found in our O star simulations (see previous section). We note that for WRstars, this number will likely be even lower owing to their higher core effective temperatures (Gräfener et al. 2012; Moens et al. 2022a) whereas for more extended and cooler LBVlike massive stars, the ratio of convective to diffusive fluxes may be significantly higher for the ironopacity peak zone (Jiang et al. 2015, 2018).
Since the majority of the radiative flux thus cannot be carried by convective transport through the ironopacity peak region, it means the radiative flux is only marginally lowered and that many gas particles there effectively have Γ > 1 . As such, they can be accelerated to above the local gas sound speed, thereby producing the high values of turbulent velocities generally seen in our simulations. This is overall consistent with the results found by Jiang et al. (2015), whose simulation StarTop is the one that is closest to the parameters of our O star models here. For this simulation, also they find inefficient convective energy transport, with radiative enthalpy carrying less than a percent of the total flux. That number is even lower than found here, but likely this is largely caused by the different choices of the stellar input parameters defining the simulations. Moreover, similar to the results of this paper, they too find large turbulent velocities that exceed the local gas sound speed in their simulations.
5.2 Photospheric macroturbulence and microturbulence
As discussed, in our models we find large atmospheric turbulent velocities, with the prototypical O4 model exhibiting ~60 km s^{−1} around the optical photosphere; the O2 model has higher values than this, and the O8 simulation somewhat lower, thus marking a clear trend with Γ_{e} (see Fig. 11). This trend is in overall good qualitative agreement with observations of photospheric absorption lines of O stars, where it is typically necessary to apply a varying amount of ‘extra linebroadening’ to match the line profiles computed by means of 1D atmospheric models to highresolution observations. Interpreted through the (ad hoc) fit parameter ‘macroturbulence’, such studies do indeed find photospheric ‘macroturbulent velocities’ that are on the same order as the turbulent velocities seen in our simulations (~30– 100 km s^{−1}), and moreover also find the same overall trend as here, that is higher macroturbulence for objects closer to the classical Eddington limit (e.g. SimónDíaz et al. 2017).
Additionally, 1D models often apply a moderate amount of (also adhoc) photospheric ‘microturbulence’ (~ 10–20 km s^{−1} in standard FASTWINDbased analysis of O stars, e.g. Hawcroft et al. 2021), both when solving for the ionisation and excitation balance and when computing synthetic spectra. A distinction between micro and macroturbulence is usually motivated by means of optical depth arguments, as essentially borrowed from spectral analysis of cool, lowmass stars (see detailed overview in Gray 2008); in practice, this just means that microturbulence is explicitly included in the lineprofile functions for computing number densities and synthetic spectra whereas macroturbulence is only applied afterward as a convolution that further broadens the spectral lines but preserves their equivalent widths (thereby also causing degeneracyproblems between rotational and macroturbulent linebroadening, see Sundqvist et al. 2013).
Overall, our simulations thus provide a natural explanation for the need for this ad hoc ‘extra linebroadening’ mechanism necessary to include in the 1D model atmosphere and spectrum synthesis codes (which in a sense mirrors results found already for lowmass sunlike stars, e.g. Asplund et al. 2000). It remains to be seen, however, how the large turbulent velocities observed in our simulations may (or may not) be captured in such standard 1D methods for performing quantitative spectroscopic analysis of O stars with winds. Certainly, quantitative comparisons would require spectral line synthesis calculations directly from our simulations, investigating not only the overall broadening of lines stemming from our models but also their predicted strengths, shapes, and variability due to complex dynamics. Some promising first results in this respect have been published by Schultz et al. (2022), based on the simulations by Jiang et al. (2015, 2018), thus neglecting the influence of the wind outflow.
5.3 Wind turbulence and clumping
The characteristic turbulent velocities in our simulations increase even further in the supersonic wind outflow. Again this qualitatively agrees with observations of UV resonance wind lines in O stars (socalled ‘P Cygni’ lines), which typically require turbulent velocities (usually included in the lineprofile function as microturbulence, see above) that increase with wind velocity when fitted by means of 1D, stationary models, reaching terminal values on order ~(0.1–0.2)v_{∞} (e.g, Haser et al. 1995; Brands et al. 2022; Hawcroft et al. 2024).
Clumping factors ƒ_{cl} = 〈ρ^{2}〉〈ρ〉^{2} range between ~2–10 in the outflowing wind parts of our simulations. While these are quite reasonable numbers, in this context it is important to recall that our method neglects the influence of the strong linedeshadowing instability (LDI, Owocki & Rybicki 1984; Owocki et al. 1988; Rybicki et al. 1990, see further discussion in Sect. 5.5), which is believed to cause additional strong structure formation in O star winds. Moreover, for both Sobolevbased (Moens et al. 2022a) and LDIbased (Sundqvist et al. 2018; Driessen et al. 2022b) calculations it has been found that increasing the spatial dimensionality of the simulation typically decreases the quantitative values of ƒ_{cl} in the wind. As such, we defer further more detailed quantifications of clumping factors in our simulations to future 3D modeling.
Nevertheless, it is likely that the most important finding regarding clumping in our simulations concerns the ‘twocomponent medium’ ansatz applied in (to our knowledge) all major spectroscopic 1D modeling attempts to account for wind clumping (Hillier 1991; Oskinova et al. 2007; Sundqvist et al. 2010; Sander et al. 2017; Sundqvist & Puls 2018). Essentially, all such methods had previously relied on the assumption that the wind consists of smallscale density clumps embedded in a very rarefied (almost void) medium so that (almost) the whole wind mass is contained within the clumps (and the wind excitation and ionisation balance then can be derived only for the clumped medium). The recent studies by Hawcroft et al. (2021), Brands et al. (2022), and Verhamme et al. (in prep.), however, all find strong empirical indications that the ‘interclump medium’ seems to be much less rarefied than previously thought, and that a significant fraction of the total wind mass might not reside in overdense, smallscale clumps. Indeed, inspection of the probability colour maps in this paper really does not seem to find much basic support for a twocomponent medium ansatz; rather the most likely winddensity is close to the mean, with a large dispersion around it, and further with a strong anticorrelation between density and wind velocity. This was also found in the 3D WR models by Moens et al. (2022a; see discussion therein) and similar distributions that are rather centered around the mean have actually also been observed in recent 2D LDI simulations (Driessen et al. 2022a).
While these preliminary notions will need further verification from more detailed and dedicated studies, in view of the above it is certainly the case that current assumptions underlying the treatment of wind clumping in presentday 1D model atmosphere codes are in need of a careful reinvestigation. Owocki & Sundqvist (2018) recently made a first such attempt to relax the standard twocomponent medium assumption, but further work is required since that study focused solely on continuum absorption, whereas the most important effects typically are observed for spectral line formation.
5.4 General implication for 1D stellar atmosphere with wind models and spectral synthesis
The most significant difference in the average photospheric structure of our 2D simulations and those of the 1D models are most probably related to the markedly different slopes of gas density and temperature around the optical photosphere (see Fig. 14). In our simulations, this effect scales with Eddington parameter Γe, with larger differences for the O2 and O4 models than for the less luminous O8 model. Additionally, the 〈2D〉 models also display a significantly more complex average velocity structure around the photosphere and in the wind initiation region than usually assumed in the 1D models used for quantitative spectroscopy, as well as wind clumping properties that generally do not seem to support the basic assumptions normally adopted in 1D codes to account for such effects (see above).
These differences will likely affect the formation process of essential spectral lines, thereby also affecting empirical spectroscopic derivation of fundamental stellar parameters (T_{eff}, g_{★}) as well as determinations of chemical abundances in the O star regime. Moreover, general calibrations and comparisons between evolution predictions and spectroscopic results may also be significantly affected, such as the socalled ‘massdiscrepancy problem’ for O stars (e.g. Herrero et al. 1992).
As shown previously, in order to reasonably match the 〈2D〉 and 1D density and temperature structures, we need to in the latter add a large turbulent pressure $\overline{{P}_{\text{turb}}}$ in the hydrostatic equation, and, moreover, a moderate $\overline{{F}_{\text{conv}}}\text{}$ in the energy equation to get rid of density inversions around the iron opacity peak. It thus remains to be tested if these effects can be captured also without a downward extension of the lower boundary in current 1D O star atmospheric models such as FASTWIND, PoWR, CMFGEN, amongst others (which typically lies well above this region, see Sect. 2.3).
5.5 Limitations of our simulations
To make these first multidimensional unified atmosphere and wind simulations for O stars possible, there are a number of simplifications we have made in our simulation techniques. Although we believe that none of them will drastically change the general nature of the results presented here, it will be an important focus of future work to try and improve upon some of them.
First, as mentioned previously, we assume that the flux, Planck, and energyweighted mean opacities are the same. In particular, this may affect $\dot{q}$, namely, the balance between radiative heating and cooling, since really only the pure absorption part of the total opacity should be accounted for (see also discussion in Moens et al. 2022b). While similar approximations (i.e. using the full κ^{line} for energy and Planck means) have been quite successfully applied to approximate effects of ‘line blanketing’ in stationary 1D modeling of the temperature structure in optically thick WR winds (Lucy & Abbott 1993), it is likely that our current approach overestimates $\dot{q}$ in the windparts of the O stars studied here. Essentially, our present methodology introduces a strong wind blanketing effect that forces the radiation and gas temperatures to be almost identical; that is, our simulations effectively behave like LTE models. However, while this approximation has provided good estimates for the line force even in the O star regime (e.g. Poniatowski et al. 2022), it is more questionable for computing $\dot{q}$ in such O star winds. As such, we have very recently attempted to extend our method for computing κ^{line} to (approximate) nonLTE conditions (following Abbott & Lucy 1985; Lucy & Abbott 1993; Puls et al. 2000), separating the flux, energy, and Planck mean opacities. This method, and the first results regarding effects on the predicted energy balance and structure of our simulations, will be presented in an upcoming paper.
Secondly, as already discussed in Sect. 2.2 our hybrid opacity formalism may doublecount line opacities in an ‘intermediate’ region wherein the Sobolevbased κ^{line} has already started to become effective but κ^{Ross} still have contributions from line opacity. Since such situation may indeed occur around the regions where the wind is launched, it is natural to ask if this might significantly affect our predicted wind characteristics. In this respect, however, the fact that our average massloss rates overall agree well with those predicted by detailed 1D models based on full nonLTE and multifrequency comoving frame (CMF) radiative transfer (Fig. 18) suggests that the impact is not severe, at least not in an average sense for the regime investigated in this paper. Specifically, as seen in Fig. 18 we obtain good agreement with the massloss rates predicted by Björklund et al. (2021). These models are based on the iterative method by Sundqvist et al. (2019), using CMF solutions over the entire relevant frequencyrange (typically between ~200 and 10 000 Å, see Puls et al. 2020 for a detailed description) when computing the radiatvie acceleration, and also use the same line list as employed in this work. That is, if either our Sobolevbased linedriving calculations or the issues related to adding opacities had been severe, we would likely not find such good average agreements with these very detailed (albeit 1D and stationary) models. We also find equally good agreement with the alternative CMF linedriven wind models by Krtička & Kubát (2017), which have been developed completely independently of ours and also use different line lists. Nevertheless, it will be an important task in future work to try and further improve the hybrid opacity approximation employed here.
Thirdly, in our calculations for the κ^{line} we assume the Sobolev approximation, which inevitably neglects the influence of the LDI. While it is believed that this LDI might be strongly damped in optically thick atmospheric regions (Gayley & Owocki 1995), it will likely further increase structure formation in parts of the wind that are optically thin. However, a proper inclusion of the LDI would require that we replace our formulation of κ^{line} in terms of the local velocity gradient and instead perform (very costly) nonlocal integrations to obtain the corresponding optical depths. Typically, all LDI simulations carried out thus far have been based on some variant of the escapeintegral formalism of Owocki & Puls (1996), which is different than the method presented and applied here; as such, it remains to be seen if and how our current modeling framework could be adapted to include the LDI effect.
Fourthly, to solve the RHD equations in our simulations we use analytical closure relations between radiation energy density, pressure, and flux. While this analytic closure reproduces the correct optically thick and thin limits, it is again in the ‘intermediate’ regime, where it may be problematic. To evaluate the potential errors introduced by our analytical approximation, Moens et al. (2022a) calculated the radiation quantities by solving the 3D radiative transfer equation using a short characteristics (SC) method (Hennicker et al. 2020), which illustrated that (except in the outer winds) the analytical closure approximations overall actually yielded very similar results to such full 3D radiative transfer solutions (see further discussion in Moens et al. 2022a). While this suggests that our current approximate ‘bridging laws’ to close the radiation equations might not be too bad, it will nonetheless be important in future work to try and replace these (or at least further calibrate them) with full solutions to the 3D radiative transfer equation (for example via a variable Eddington tensor closure, see, e.g. Jiang et al. 2015). How to effectively account for the accumulative effect of linedriving stemming from thousands of spectral lines within such an approach remains a challenging question for multidimensional simulations though.
Finally, in this study, we have only performed 2D timedependent simulations. Although we do not believe this affects our overall results and conclusions, some quantitative properties may change when introducing the third spatial dimension, as has been already shown in Moens et al. (2022a) for 3D WR stars. In a future study, we will extend also our O star simulations to full 3D. Similarly, in such upcoming studies, we also aim to consider the effects of rotation and magnetic fields, which are both neglected here.
6 Summary and outlook
In this work, we have presented the first multidimensional, timedependent unified atmosphere with wind models for O stars. To this end, we use a ‘boxinastar’ approach (including spherical correction terms), while solving the RHD equations using an analytical fluxlimiting closure approximation in the radiation energy equation, and accounting properly for linedriving effects. The initial conditions for our simulations are derived using procedures very similar to those applied in the standard 1D stationary model atmosphere with wind codes; however, here these are extended to deeper layers of the stellar envelope to cover the unstable subsurface ironopacity peak region. Perturbations derived from linear stability analysis are then added to the 1D gas density and lateral velocity in this region to initiate structure formation.
In our simulations, we find that structure starts appearing in the previously perturbed ironopacity peak region. Simultaneously, around and beyond the optical photosphere, the densities are low enough so that line driving becomes efficient. This initiates a net supersonic wind outflow from the variable stellar surface, introducing more structure formation and leading to a complex interplay between the deeper atmospheric layers and the overlaying wind. As the models dynamically relax, we find a (quasi) stable O star envelope beneath the iron opacity peak zone and a very turbulent atmosphere with wind above it. We find large turbulent velocities in the photospheric and wind layers of our simulations, in overall good qualitative agreement with results derived from spectroscopic observations of O stars.
The next part of our analysis compared the lateral and time ‘averaged’ (2D) atmospheric quantities with their 1D counterparts, focusing on the prototypical O4 model. Here, we find that envelope expansion in the 1D model is stronger than in 〈2D〉, leading to a lower effective temperature and larger stellar radius in the former. Moreover, the density inversion present in the 1D profile gets washed away in the 〈2D〉 profile, and the slope of the density and temperature profiles around the optical photosphere is much shallower in the 〈2D〉 simulation (indicating a much larger characteristic photospheric scale height). As such, we introduced a simplified treatment of convection and turbulence in our 1D model, treating energy transport using standard mixing length theory, accounting for the effects of radiation pressure and cooling, and introducing a turbulent pressure term in the hydrostatic momentum equation. With this, we found good qualitative agreement in the adjusted 1D and 〈2D〉 O4 model structures for parameters α_{MLT} ~ 1.0 and a v_{turb} that increased from zero in the lowermost parts to ~90 km s^{−1} in the photospheric layers. The 〈2D〉 velocity field further shows a rather complex structure in particularly the wind initiation region, going slightly negative around the photosphere before rapidly (though slightly less rapidly than in corresponding 1D stationary models) accelerating to high supersonic speeds when linedriving becomes efficient; this is due to the inverse relation between density and linedriving opacity, the outflowing wind parts also display a clear anticorrelation between density and velocity. Finally, we derived the ‘average’ mass loss rates for our models and found good agreement with the detailed (albeit 1D and stationary) theoretical mass loss calculations by Krtička & Kubát (2017) and Björklund et al. (2021); this lends further support to our basic approach for including the effects of linedriving into our simulations.
A key followup study to this work will regard the computation of synthetic spectra and spectroscopic comparison to observations. This is work underway in our group, using the 3D radiative transfer code framework by Hennicker et al. (2020). Since quantities such as the mass loss rates, wind and turbulent velocities, and ‘atmospheric clumping’ in our multidimensional simulations are emergent properties (rather than adjustable free input parameters), these studies will provide key basic tests regarding the realism of our new O star unified model atmosphere calculations. Furthermore, they will directly test whether the observed strength and broadness of photospheric and wind lines are also quantitatively matched by our models, how this may affect spectroscopic derivation of fundamental stellar parameters and, as such, provide further insights into how presentday 1D codes might have to be calibrated to account for basic multidimensional effects. Moreover, since the photosphere and wind are variable in our O star simulations, this will likely cause both spectral line and stochastic photometric variability (see also Schultz et al. 2022).
As discussed in Sect. 5.5, in another upcoming work we will also separate the flux, Planck, and energyweighted mean opacities, using an approximate nonLTE method based primarily on Puls et al. (2000). Finally, the methods developed in this paper (also, Moens et al. 2022a) are generally applicable also to other interesting regions of the upper Hertzsprung–Russel (HR) diagram, for instance: Bsupergiants and LBVs.
Acknowledgements
The computational resources used for this work were provided by Vlaams Supercomputer Centrum (VSC) funded by the Research FoundationFlanders (FWO) and the Flemish Government. D.D., C.S., N.M., L.P., and J.S. acknowledge the support of the European Research Council (ERC) Horizon Europe under grant agreement number 101044048. O.V., N.M., L.P., and J.S. acknowledge support from KU Leuven C1 grant MAESTRO C16/17/007, the Belgian Research Foundation Flanders (FWO) Odysseus program under grant number G0H9218N. FWO grant G077822N. J.S. would further like to extend his thanks to Stan Owocki and Jo Puls for many fruitful discussions over the years, and to the latter also for providing his easytofollow (albeit in handwritten German, of course) koch rezept for MLTbased convective energy transport including radiative pressure and cooling. We thank the referee for their constructive comments to the manuscript. Finally, the authors would like to thank all members of the KULEQUATION group for fruitful discussion, comments, and suggestions We made significant use of the following packages to analyze our data: NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), matplotlib (Hunter 2007), Python amrvac_reader (Keppens et al. 2021).
Appendix A Castor’s expansion opacity
The expansion opacity model uses the escape probability concept from Sobolev (1960) in a situation when a spectrum is populated with multiple lines and the medium has a significant velocity gradient (Karp et al. 1977; Castor 2004). Under these conditions, for lines with a frequency separation corresponding to a velocity interval Δv < v_{∞}, the photon meanfreepath between interactions with different lines is ℓ_{mfp} = 1/(κ_{eff}ρ), where κ_{eff} is the effective opacity. Here l_{mfp} = Δl/P for Δl the physical mean distance between two lines and $P=1{e}^{{\tau}_{s}}$ the line transition probability. By means of the Doppler shift and Sobolev line optical depth, we obtain: $${\tau}_{s}=\frac{{\alpha}_{l}c}{{v}_{0}}{\left\frac{d{v}_{l}}{dl}\right}^{1},$$(A.1)
and the effective line opacity is thus: $${\kappa}_{\text{eff}}=\frac{1}{\rho}\frac{{\alpha}_{l}}{\Delta v}\frac{1{\text{e}}^{{\tau}_{s}}}{{\tau}_{s}},$$(A.2)
where α_{l} is the frequencyintegrated line extinction coefficient. Considering a frequency interval Δv within which there is i number of lines (in a line list), and adding a continuum opacity source given by Thompson scattering κ_{TH}, the total effective opacity is $${\kappa}_{\text{eff}}=\frac{1}{\rho}{\displaystyle \sum _{i}^{\text{linesin}\Delta v}\frac{{\alpha}_{l}^{i}}{\Delta v}}\left(\frac{1{\text{e}}^{{\tau}_{s}^{i}}}{{\tau}_{s}^{i}}\right)+{\kappa}_{\text{TH}}$$(A.3) $$=\frac{d{v}_{l}}{dl}\frac{1}{c\rho}{\displaystyle \sum _{i}^{\text{linesin}\Delta v}\frac{{v}^{i}}{\Delta v}}\left(1{\text{e}}^{{\tau}_{s}^{i}}\right)+{\kappa}_{\text{TH}}.$$(A.4)
The latter expression is equivalent to Castor (2004)’s Eq. 6.130 (also Eq. 9 in Friend & Castor 1983) for the effective opacity.
Furthermore, Castor (2004) performs explicit opacity calculations using a Fe III line list in LTE (his Sect. 6.9.4). He finds that for large values of τ = τ_{S} /q_{i} (τ = s in his notation), the above formula considerably underestimates the effective opacity since the Sobolev line opacity tends to zero when the velocity gradient becomes very small, due to the neglect of the intrinsic line width. To make the expansion opacity formula agree better with explicit calculations, thus, the author suggests that the Thomson continuum be replaced with the actual harmonic Rosseland mean opacity calculated without the velocity gradient (but including the intrinsic line widths); namely, $${\kappa}_{\text{eff}}=\frac{1}{\rho}{\displaystyle \sum _{i}^{\text{linesin}\Delta v}\frac{{\alpha}_{l}^{i}}{\Delta v}}\left(\frac{1{\text{e}}^{{\tau}_{s}^{i}}}{{\tau}_{s}^{i}}\right)+{\kappa}^{\text{Ross}}$$(A.5)
Connection to hybrid opacity model. The fluxweighted opacity of the above expansion opacity is: $$\frac{{\kappa}_{\text{eff}}{F}_{v}\Delta v}{F}$$(A.6)
when taken over the complete frequency range, this becomes: $$\kappa ={\kappa}_{0}{\displaystyle \sum _{i}^{\text{alllines}}{q}_{i}}{w}_{i}\left(\frac{1{\text{e}}^{{\tau}_{s}^{i}}}{{\tau}_{s}^{i}}\right)+{\kappa}^{\text{Ross}}$$(A.7)
which we have translated to notation used in the present paper; this illustrates how the hybrid opacity model suggested by Poniatowski et al. (2022) in effect is equivalent to the suggestion in Castor (2004) for modification of the Sobolev expansion opacity to account for the intrinsic widths of the lines in the limit of large values for τ (∝ ρ/dv_{l}/dl).
References
 Abbott, D. C., & Lucy, L. B. 1985, ApJ, 288, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Nordlund, A., Trampedach, R., Allende Prieto, C., & Stein, R. F. 2000, A&A, 359, 729 [NASA ADS] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, 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]
 Blaes, O., & Socrates, A. 2003, ApJ, 596, 509 [CrossRef] [Google Scholar]
 Brands, S. A., de Koter, A., Bestenlehner, J. M., et al. 2022, A&A, 663, A36 [Google Scholar]
 Cantiello, M., Langer, N., Brott, I., et al. 2009, A&A, 499, 279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Castor, J. I. 2004, Radiation Hydrodynamics (Cambridge University Press) [CrossRef] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [Google Scholar]
 Conti, P. S., & Ebbets, D. 1977, ApJ, 213, 438 [CrossRef] [Google Scholar]
 Driessen, F., Sundqvist, J. O., & Kee, N. D. 2022a, PhD Thesis, KU Leuven, Belgium [Google Scholar]
 Driessen, F. A., Sundqvist, J. O., & Dagore, A. 2022b, A&A, 663, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eversberg, T., Lépine, S., & Moffat, A. F. J. 1998, ApJ, 494, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Freytag, B., Steffen, M., Ludwig, H. G., et al. 2012, J. Comput. Phys., 231, 919 [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]
 Gabler, R., Gabler, A., Kudritzki, R. P., Puls, J., & Pauldrach, A. 1989, A&A, 226, 162 [NASA ADS] [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [Google Scholar]
 Gayley, K. G., & Owocki, S. P. 1995, ApJ, 446, 801 [NASA ADS] [CrossRef] [Google Scholar]
 Glatzel, W. 1994, MNRAS, 271, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Goldberg, J. A., Jiang, Y.F., & Bildsten, L. 2022, ApJ, 929, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Gräfener, G., Koesterke, L., & Hamann, W. R. 2002, A&A, 387, 244 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gräfener, G., Owocki, S. P., & Vink, J. S. 2012, A&A, 538, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grassitelli, L., Chené, A. N., Sanyal, D., et al. 2016, A&A, 590, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grassitelli, L., Langer, N., Grin, N. J., et al. 2018, A&A, 614, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gray, D. F. 2008, The Observation and Analysis of Stellar Photospheres (Cambridge University Press) [Google Scholar]
 Grevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, eds. N. Prantzos, E. VangioniFlam, & M. Casse (Cambridge University Press), 15 [Google Scholar]
 Hamann, W. R., & Gräfener, G. 2004, A&A, 427, 697 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Haser, S. M., Lennon, D. J., Kudritzki, R. P., et al. 1995, A&A, 295, 136 [NASA ADS] [Google Scholar]
 Hawcroft, C., Sana, H., Mahy, L., et al. 2021, A&A, 655, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hawcroft, C., Sana, H., Mahy, L., et al. 2024, A&A, in press, https://doi.org/10.1051/00046361/202245588 [Google Scholar]
 Hearn, A. G. 1972, A&A, 19, 417 [NASA ADS] [Google Scholar]
 Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2020, A&A, 633, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Herrero, A., Kudritzki, R. P., Vilchez, J. M., et al. 1992, A&A, 261, 209 [NASA ADS] [Google Scholar]
 Hillier, D. J. 1991, A&A, 247, 455 [NASA ADS] [Google Scholar]
 Hillier, D. J., & Lanz, T. 2001, in ASP Conf. Ser., 247, Spectroscopic Challenges of Photoionized Plasmas, eds. G. Ferland, & D. W. Savin, 343 [NASA ADS] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Howarth, I. D., Siebert, K. W., Hussain, G. A. J., & Prinja, R. K. 1997, MNRAS, 284, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, Y.F. 2023, Galaxies, 11, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, Y.F., Cantiello, M., Bildsten, L., Quataert, E., & Blaes, O. 2015, ApJ, 813, 74 [Google Scholar]
 Jiang, Y.F., Cantiello, M., Bildsten, L., et al. 2018, Nature, 561, 498 [NASA ADS] [CrossRef] [Google Scholar]
 Karp, A. H., Lasher, G., Chan, K. L., & Salpeter, E. E. 1977, ApJ, 214, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Kee, N. D., Owocki, S., & Sundqvist, J. O. 2016, MNRAS, 458, 2323 [Google Scholar]
 Keppens, R., Teunissen, J., Xia, C., & Porth, O. 2021, Comput. Math. Appl., 81, 316 [Google Scholar]
 Keppens, R., Popescu Braileanu, B., Zhou, Y., et al. 2023, A&A, 673, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kippenhahn, R., Weigert, A., & Weiss, A. 2013, Stellar Structure and Evolution (Springer) [Google Scholar]
 Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krticka, J., & Kubát, J. 2017, A&A, 606, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B. 1971, ApJ, 163, 95 [Google Scholar]
 Lucy, L. B., & Abbott, D. C. 1993, ApJ, 405, 738 [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (Oxford University Press) [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]
 Oskinova, L. M., Hamann, W. R., & Feldmeier, A. 2007, A&A, 476, 1331 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owocki, S. P. 2015, in Astrophys. Space Sci. Lib., 412, Very Massive Stars in the Local Universe, ed. J. S. Vink, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & Puls, J. 1996, ApJ, 462, 894 [Google Scholar]
 Owocki, S. P., & Rybicki, G. B. 1984, ApJ, 284, 337 [Google Scholar]
 Owocki, S. P., & Sundqvist, J. O. 2018, MNRAS, 475, 814 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., Gayley, K. G., & Shaviv, N. J. 2004, ApJ, 616, 525 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., Townsend, R. H. D., & Quataert, E. 2017, MNRAS, 472, 3749 [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 ASP Conf. Ser., 131, Properties of Hot Luminous Stars, ed. I. Howarth, 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, in ASP Conf. Ser., 214, IAU Colloq. 175: The Be Phenomenon in EarlyType Stars, eds. M. A. Smith, H. F. Henrichs, & J. Fabregat, 626 [NASA ADS] [Google Scholar]
 Petrovic, J., Pols, O., & Langer, N. 2006, A&A, 450, 219 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Poniatowski, L. G., Kee, N. D., Sundqvist, J. O., et al. 2022, A&A, 667, A113 [NASA ADS] [CrossRef] [EDP Sciences] [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., Markova, N., Scuderi, S., et al. 2006, A&A, 454, 625 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Najarro, F., Sundqvist, J. O., & Sen, K. 2020, A&A, 642, A172 [EDP Sciences] [Google Scholar]
 Rogers, F. J., & Iglesias, C. A. 1992, ApJS, 79, 507 [CrossRef] [Google Scholar]
 Rybicki, G. B., Owocki, S. P., & Castor, J. I. 1990, ApJ, 349, 274 [NASA ADS] [CrossRef] [Google Scholar]
 Sander, A., Hamann, W.R., & Todt, H. 2012, A&A, 540, A144 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 SantolayaRey, A. E., Puls, J., & Herrero, A. 1997, A&A, 323, 488 [NASA ADS] [Google Scholar]
 Sanyal, D., Grassitelli, L., Langer, N., & Bestenlehner, J. M. 2015, A&A, 580, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schultz, W. C., Bildsten, L., & Jiang, Y.F. 2020, ApJ, 902, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Schultz, W. C., Bildsten, L., & Jiang, Y.F. 2022, ApJ, 924, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Schwarzschild, M. 1958, Structure and Evolution of the Stars (Princeton University Press) [Google Scholar]
 Shaviv, N. J. 1998, ApJ, 494, L193 [NASA ADS] [CrossRef] [Google Scholar]
 SimónDíaz, S., Godart, M., Castro, N., et al. 2017, A&A, 597, A22 [CrossRef] [EDP Sciences] [Google Scholar]
 Sobolev, V. V. 1960, Moving Envelopes of Stars (Harvard University Press) [Google Scholar]
 Stothers, R. B., & Chin, C.W. 1993, ApJ, 412, 294 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., & Puls, J. 2018, A&A, 619, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Puls, J., & Feldmeier, A. 2010, A&A, 510, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., SimónDíaz, S., Puls, J., & Markova, N. 2013, A&A, 559, A10 [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., Kubát, J., Oskinova, L. M., & Feldmeier, A. 2012, A&A, 541, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Teunissen, J., & Keppens, R. 2019, Comput. Phys. Commun., 245, 106866 [Google Scholar]
 Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [NASA ADS] [CrossRef] [Google Scholar]
 ud Doula, A., & Owocki, S. P. 2002, ApJ, 576, 413 [NASA ADS] [CrossRef] [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Xia, C., Teunissen, J., Mellah, I. E., Chané, E., & Keppens, R. 2018, ApJS, 234, 30 [NASA ADS] [CrossRef] [Google Scholar]
To make suitable comparisons, we have here adjusted the standard input specifications of more recent versions of FASTWIND (Puls et al. 2005; Sundqvist & Puls 2018) by opting for the ‘Lucy temperature structure’ (Lucy 1971) without further corrections, neglecting any wind clumping effects, and choosing a solar abundance scale reflecting the calculations of this paper.
All Tables
All Figures
Fig. 1 Line force multiplier M(τ) as a function of the characteristic Sobolev optical depth τ. The blackdashed line is plotted by calculating M(τ) at temperature 30 kK and density 10^{−10} g cm^{−3}, whereas the red line is calculated by fitting M(τ) using M_{G}(τ) in Eq. (13). The best fit for the line force parameters for this particular temperature and density are $\overline{Q}=2278$, Q_{0} = 1418, and α = 0.59. 

In the text 
Fig. 2 Averaged opacity 〈κ〉 as a function of scaled radius x = 1 − R_{★}/r, with R_{★} as the 2D averaged optical photosphere (reddashed line). In the figure, the brown curve is the Rosseland mean opacity and the purple curve is the finitedisk corrected lineopacity. The black dashdotted curve is the total opacity using the hybrid opacity scheme. 

In the text 
Fig. 3 Initial 1D gas density (in g cm^{−3}) structure as a function of modified radius x_{0} = 1 − R_{0}/r, with R_{0} as the lower boundary radius. The cyan curve is corresponding FASTWIND density structure as described in Sect. 2.3. The reddashed vertical line is the optical photosphere R_{★}. The shadedpurple region is the instability region where we provide the initial perturbations (see zoomedin part). 

In the text 
Fig. 4 Number of total pressure height scale ${H}_{\text{p}}^{\text{tot}}$ resolved per cell in our simulation. Level 4 corresponds to the highest level of resolution in our models. Levels 3–1 are not shown since it is zoomed in to indicate the deep layers and the photosphere. See text for details. 

In the text 
Fig. 5 Mass loss rates (top) and effective temperature (bottom) for the O4 model as varying with time (shown here in 10^{3} seconds). The black curve displays the laterally averaged quantity varying across time written as 〈2D〉_{L}, and the blackdashed line is the lateraltemporal averaged quantity, here 〈2D〉. 

In the text 
Fig. 6 Colour maps of relative radiation temperature (top panel), relative density (second panel), and the radial velocity (third panel) and Γ (bottom panel), for several snapshots at different times (shown on top in 10^{3} seconds) at the beginning of our simulation for the O4 model. The figures are zoomed in to show the inner regions of our simulations. In the lateral direction, the figure extends to 0.2 R_{0} and in the radial direction to 1.5 R_{0}, with R_{0} as the lower boundary radius of our simulation. 

In the text 
Fig. 7 Similar to Fig. 6, but now at times after the models have dynamically relaxed from their initial conditions. The blackdashed line is the averaged 2D optical photosphere 〈R_{★}〉. 

In the text 
Fig. 8 Running averaged total luminosity 〈L_{tot}〉 as a function of scaled radius x = 1 − R_{★}/r. The thin lines represent lateral averages taken over a few time snapshots (with lighter colours representing earlier snapshots), whereas the black line represents a long timeaverage of the total luminosity. 

In the text 
Fig. 9 Running averages of the radial velocities for the O4^{*} model (top) and O4 model (bottom), taken over 10 snapshots starting from the beginning of the simulation, with lighter colours representing earlier times. The red dashdotted lines are the initial velocity profiles. The green dashdotted line is the final average for the O4^{*} model and the orange dashed line is for the O4 model. See text. 

In the text 
Fig. 10 Probability density maps for different quantities, from top to bottom radial velocity, gas density, Eddington Γ, and radiation temperature for different models O8 (left), O4 (middle), and O2 (right). 40 snapshots have been used to create these probability density maps. The red curves represent the lateraltemporal averages. The vertical orangedashed line is 〈R_{★}〉, whereas the green dashdotted vertical lines in the top panel are 〈R_{sonic}〉. The cyan horizontal lines in the Γ panels are where its value is unity. The colours signify the probability of finding a quantity at a certain cell, with yellow having a higher probability and blue having a lower one. 

In the text 
Fig. 11 2D average radial velocity (dashed lines) and its dispersion (solid lines) as a function of the scaled radius x = 1 − R_{★}/r. The colours correspond to the different models. 

In the text 
Fig. 12 Average relative density for each radial velocity at different radii for the O4 model. 〈ρ_{r}〉 is the average density at every radial cell. Colours here represent relative density (ρ/〈ρ_{r}〉) with yellow displaying underdense regions and blue representing overdense clumps. The red dashdotted line indicates the local escape velocity, the black dashed vertical line is the average photosphere 〈R_{★}〉 and the black dashdotted vertical line is the average sonic point 〈R_{sonic}〉. The zoomedin part highlights the regions below the optical photosphere. 

In the text 
Fig. 13 Initial 1D and resulting 〈2D〉 gas density [in g cm^{−3}] (top) and temperature [in K] (bottom) structure for the O4 model as a function of modified radius x = 1 − R_{0}/r, with R_{0} as the lower boundary radius. The ‘red dashdotted’ line is the 〈2D〉 photosphere and the ‘blue dashdotted’ is the 1D photosphere. The cyan dashdotted curve is produced using FASTWIND as described with specifications in Sect. 2.3. 

In the text 
Fig. 14 〈2D〉 gas density (in g cm^{−3}; top panel) and temperature (in K; bottom panel) as a function of the scaled radius x = 1 − R_{★}/r for the O4 model. The 1D structure is calibrated to the 〈T_{eff}〉 and 〈R_{★}〉 obtained from the 2D model. The orange line displays the calibrated 1D model with an α_{MLT} = 1 and a turbulent velocity of 90 km s^{−1} at the photosphere and reduced gradually in the deeper regions close to the hydrostatic envelopes. The cyan dashdotted curve is produced using FASTWIND, using the specifications outlined in Table 1 and Sect. 2.3. 

In the text 
Fig. 15 Convective flux to the total flux as a function of temperature for the O4 model. The black curve is the 〈2D〉, and the orange curve is the calibrated 1D model to match the 2D with α_{MLT} = 1. 

In the text 
Fig. 16 2D average gas pressure, radiation pressure, and turbulent pressure (top panel) and turbulent velocity, dispersion of radial velocity, and radial velocity (bottom panel) for the O4 model as a function of the scaled radius x = 1 − R_{★}/r with the upper axis as radiation temperature in 10^{3} K. 

In the text 
Fig. 17 Fluxporosity effect as a function of scaled radius for our O4 model. 

In the text 
Fig. 18 Mass loss rates in M_{⊙} yr^{−1} as a function of solar luminosity. The ★markers are the mass loss rates calculated from the models considered in this work, whereas the dotdashed line using Krtička & Kubát (2017), and the dashed line is calculated using Björklund et al. (2021). The triangle markers are plotted from the empirical mass loss rate study by Hawcroft et al. (2021; their optically thick GA rates). 

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.