Open Access
Issue
A&A
Volume 684, April 2024
Article Number A177
Number of page(s) 21
Section Stellar atmospheres
DOI https://doi.org/10.1051/0004-6361/202348206
Published online 23 April 2024

© The Authors 2024

Licence Creative CommonsOpen 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 line-driven outflow (Castor et al. 1975), atmospheric models of massive O-type stars (O stars) are required to treat the sub-surface 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, Santolaya-Rey 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ón-Dí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 time-dependent, radiation-hydrodynamical (RHD) simulations of the so-called ‘iron-opacity 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 line-opacity in a medium with significant velocity gradient (the ‘line-driving’ 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 opacity-formalism able to treat the sub-surface layers and the line-driven wind in a unified way and then applied this in multi-dimensional simulations of the atmospheres and outflows of classical Wolf-Rayet stars (Moens et al. 2022a).

In this work, we present the first results from our subsequent attempts to build a unified multi-dimensional 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 iron-opacity peak region, and then perturbed them and used them as initial conditions to solve the time-dependent RHD equations. We use the MPI-AMRVAC 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 MPI-AMRVAC (Moens et al. 2022b), and used the ‘Munich Atomic Database’ (Pauldrach et al. 1998; Puls et al. 2000) to calculate self-consistent line-opacities 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 sub-surface and atmospheric structures of standard 1D models. In follow-up 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 opacity-formalism 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 Radiation-hydrodynamics

The previous studies by Moens et al. (2022b,a) developed radiation-hydrodynamic (RHD) modules of the general hydrodynamics code MPI-AMRVAC (Xia et al. 2018) and used these to study multi-dimensional Wolf-Rayet stellar atmospheres and winds. Following their model and study, we also solve numerically the time-dependent RHD equations on a finite volume mesh with a ‘box-in-a-star’ 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: tρ+(ρv)=0,${\partial _t}\rho + \nabla \cdot (\rho {\bf{v}}) = 0,$(1) t(vρ)+(vρv+pgI)=fgrav +frad ,${\partial _t}({\bf{v}}\rho ) + \nabla \cdot \left( {{\bf{v}}\rho {\bf{v}} + {p_{\rm{g}}}} \right) = - {{\bf{f}}_{{\rm{grav }}}} + {f_{{\rm{rad }}}},$(2) teg+(egv+pgv)=fgravv+fradv+q˙.${\partial _t}{e_{\rm{g}}} + \nabla \cdot \left( {{e_{\rm{g}}}{\bf{v}} + {p_{\rm{g}}}{\bf{v}}} \right) = - {f_{{\rm{grav}}}} \cdot {\bf{v}} + {f_{{\rm{rad}}}} \cdot {\bf{v}} + \dot q.$(3)

Here, v is the velocity of the gas, ρ is the gas density, pg is the gas pressure, 𝕀 is a unit tensor, and eg is the total gas energy which comprises internal and kinetic components given by, eg=pgγg1+ρv22,${e_{\rm{g}}} = {{{p_{\rm{g}}}} \over {{\gamma _{\rm{g}}} - 1}} + {{\rho {v^2}} \over 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 right-hand 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 q˙$\dot q$. Since q˙$\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 field1 as: tErad+(Eradv+Fdiff)=q˙v:Prad.${\partial _t}{E_{{\rm{rad}}}} + \nabla \cdot \left( {{E_{{\rm{rad}}}}{\bf{v}} + {{\bf{F}}_{{\rm{diff}}}}} \right) = - \dot q - \nabla {\bf{v}}:{{\bf{P}}_{{\rm{rad}}}}$(5)

Here, E(4σ/c)Trad4=arTrad4$E \equiv (4\sigma /c)T_{{\rm{rad}}}^4 = {a_r}T_{{\rm{rad}}}^4$ represents the frequency-integrated radiation energy density for a radiation temperature Trad. Also, Prad is the frequency-integrated radiation pressure tensor and ar 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 ƒ = Prad/E. Additionally, v : Prad on the right-hand side of Eq. (5) is the dyadic product between the gradient of the velocity vector and Prad. 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 Fdiff is derived from Fdiff =cλκρE=DE,${F_{{\rm{diff }}}} = - {{c\lambda } \over {\kappa \rho }}\nabla E = - D\nabla E,$(6)

where κ is the flux-weighted opacity, c the speed of light, and λ a flux-limiter obtained from an approximation that bridges the analytic and opposite optically thick and free-streaming 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 frad=ρκFdiff c.${f_{{\rm{rad}}}} = \rho {{\kappa {{\bf{F}}_{{\rm{diff }}}}} \over c}.$(7)

As in Moens et al. (2022b,a), throughout this paper, we assume the flux-, Planck- and energy-weighted mean opacities to be the same, given by Eq. (10). Thus q˙=cκρE4πκρB,$\dot q = c\kappa \rho E - 4\pi \kappa \rho B,$(8)

where BσTg4/π$B \equiv \sigma T_{\rm{g}}^4/\pi $ is the frequency-integrated Planck function for gas temperature Tg = pgµmH/(kbρ). We assume a constant mean molecular weight µ = 0.61 and our method allows for different radiation and gas temperatures. kb is the Boltzmann constant and mH is the mass of a proton. Here, the gravitational acceleration is described by: fgrav =ρGMr2r^=ggr^${f_{{\rm{grav }}}} = \rho {{G{M_ \star }} \over {{r^2}}}\widehat {\bf{r}} = {g_{\rm{g}}}\widehat {\bf{r}}{\rm{, }}$(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 ‘line-driving’ (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: κ=κRoss +κline .$\kappa = {\kappa ^{{\rm{Ross }}}} + {\kappa ^{{\rm{line }}}}.$(10)

Here, κ is the total opacity of the medium, κline is the opacity due to lines in a supersonic medium (line-driving 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 × 106 lines, Pauldrach et al. 1998, 2001), such that κline =f*κ0iwv,iqi(1eqiτqiτ)=f*κ0M(τ),${\kappa ^{{\rm{line }}}} = {f^*}{\kappa _0}\sum\limits_{\rm{i}} {{w_{v,{\rm{i}}}}} {q_{\rm{i}}}\left( {{{1 - {e^{ - {q_{\rm{i}}}\tau }}} \over {{q_{\rm{i}}}\tau }}} \right) = {f^*}{\kappa _0}M(\tau ),$(11)

where qi is the line-strength, 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: τ=κ0cρ| dvr dr |1$\tau = {\kappa _0}c\rho {\left| {{{{\rm{d}}{v_r}} \over {{\rm{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 so-called 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, τ (∝ ρ/|dvr/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 < qiτ ≪ 1) → iwv,iqiQ¯$\sum\limits_i {{w_{v,i}}} {q_i} \equiv \bar Q$, where we have introduced Gayley (1995)’s Q¯$\bar Q$ (also see below). Following Poniatowski et al. (2022) this becomes Q¯=iBv,iαl,i/(Bρ)$\bar Q = \sum\limits_i {{B_{v,i}}} {\alpha _{l,i}}/(B\rho )$, where αl is the frequency-integrated 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 non-overlapping) lines and a Thomson scattering continuum. In an ‘intermediate’ region, however, lines could in principle be ‘double-counted’ 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 mass-loss rates of the simulations presented here and those computed for O stars by means of full frequency-dependent comoving-frame radiative transfer (in 1D and assuming a steady-state), 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 Sobolev-based ‘expansion opacity’ formulae suggested by Castor (2004) based on calculations of Fe III lines.

The results from the complete line-ensemble summation are fitted using (Gayley 1995): κline =f*κ0Q¯1α((1Q0τ)1α1)Q0τ=f*κ0MG(τ)${\kappa ^{{\rm{line }}}} = {f^*}{\kappa _0}{{\bar Q} \over {1 - \alpha }}{{\left( {{{\left( {1 - {Q_0}\tau } \right)}^{1 - \alpha }} - 1} \right)} \over {{Q_0}\tau }} = {f^*}{\kappa _0}{M_{\rm{G}}}(\tau ){\rm{. }}$(13)

Here, we choose the normalisation κ0 = 0.34 cm2 g−1, which reflects the Thomson scattering opacity in a fully ionised plasma with a chemical composition similar to that of the Sun. Also, Q¯$\bar Q$, Q0, and α are the line force fit parameters; Q¯$\bar Q$ sets the maximum value of the line opacity in the limit where all lines are optically thin (see above), whereas Q0 and α give an effective maximum line strength and power-law 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 MG(τ) to the integrated M(τ) calculated directly from our atomic database, namely, we obtain and tabulate α(ρ, T), Q¯(ρ,T),$\bar Q(\rho ,T),$, and Q0(ρ, T); an example is given in Fig. 1. These fitted and tabulated parameters are then used within the time-dependent 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 fit-parameters (such as Abbott’s ‘δ', see discussion in Poniatowski et al. 2022). The technique is similar to how OPAL opacity tables are constructed from pre-calculations 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 mass-fraction X = 0.70, helium mass-fraction Y = 1 − Z = 0.28 for metallicity Z = 0.02. We note that this differs slightly from the often-used 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 line-opacity component in the wind outflow above the stellar surface. To a large extent, the diminishing contribution of line-driving 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 finite-disk 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 non-radial 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 ‘over-loaded’ line-driven wind solutions in the radial streaming case. Indeed, strictly speaking, ƒ* as implemented here is valid only for spherically symmetric systems; a full multi-dimensional 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 finite-disk correction factor has been quite successfully applied in a range of other non-spherical line-driven systems like magnetically channeled O star winds (ud Doula & Owocki 2002). Nevertheless, this approximate 1D treatment should be improved upon in future work.

thumbnail Fig. 1

Line force multiplier M(τ) as a function of the characteristic Sobolev optical depth τ. The black-dashed 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 MG(τ) in Eq. (13). The best fit for the line force parameters for this particular temperature and density are Q¯=2278$\bar Q = 2278$, Q0 = 1418, and α = 0.59.

thumbnail Fig. 2

Averaged opacity 〈κ〉 as a function of scaled radius x = 1 − R/r, with R as the 2D averaged optical photosphere (red-dashed line). In the figure, the brown curve is the Rosseland mean opacity and the purple curve is the finite-disk corrected line-opacity. The black dash-dotted curve is the total opacity using the hybrid opacity scheme.

2.3 Initial conditions

To compute the initial conditions for our time-dependent simulations, we used a procedure similar to that outlined in Sect. 2 of Santolaya-Rey 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 vr(r) = v(1 − brT/r)β and mass density ρ = Ṁ/(4πr2vr) 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 = ƒciso,gas(Teff), where ciso,gas =pgρ${c_{{\rm{iso,gas }}}} = \sqrt {{{{p_{\rm{g}}}} \over \rho }} $ is the isothermal gas sound speed, the parameter ƒ ≈ 0.1–1, and b is derived from vT; 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 Rr(τR=2/3)${R_ \star } \equiv r\left( {{\tau _{\rm{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, pg, we solve numerically the following coupled differential equations in the layers below the intermediate boundary given by rT: dpgdm=g(1Γ)(Rr)2,${{{\rm{d}}{p_{\rm{g}}}} \over {{\rm{d}}m}} = {g_ \star }(1 - \Gamma ){\left( {{{{R_ \star }} \over r}} \right)^2}$(15) dT dm=3κTeff416T3(Rr)2,${{{\rm{d}}T} \over {{\rm{d}}m}} = {{3\kappa T_{{\rm{eff}}}^4} \over {16{T^3}}}{\left( {{{{R_ \star }} \over r}} \right)^2},$(16) dτRdm=κ(Rr)2,${{{\rm{d}}{\tau _{\rm{R}}}} \over {{\rm{d}}m}} = \kappa {\left( {{{{R_ \star }} \over r}} \right)^2},$(17) dr dm=1ρ,${{{\rm{d}}r} \over {{\rm{d}}m}} = - {1 \over \rho }$(18)

along with the definitions and relations σTeff 4FL4πR2,$\sigma T_{{\rm{eff }}}^4 \equiv {F_ \star } \equiv {{{L_ \star }} \over {4\pi R_ \star ^2}}$(19) gGMR2,${g_ \star } \equiv {{G{M_ \star }} \over {R_ \star ^2}},$(20) Γgrgg=Lκ4πGMc,$\Gamma \equiv {{{g_r}} \over {{g_g}}} = {{{L_ \star }\kappa } \over {4\pi G{M_ \star }c}}$(21) pg=ρkbTμmH,${p_{\rm{g}}} = {{\rho {k_{\rm{b}}}T} \over {\mu {m_{\rm{H}}}}},$(22)

where gr 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 dPrad/dm = gr/ρ = κF/c = κL/(4πr2c), Eq. (15) above may equivalently be formulated as a hydrostatic equation for the total pressure Ptot = Prad + pg of the radiating gas. To classify models, it is further useful to introduce the classical Eddington parameter, Γe, and luminosity, Ledd: ΓeLκe4πGMcLLedd ,${\Gamma _{\rm{e}}} \equiv {{{L_ \star }{\kappa _{\rm{e}}}} \over {4\pi G{M_ \star }c}} \equiv {{{L_ \star }} \over {{L_{{\rm{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 Santolaya-Rey et al. (1997) are that we (i) do not account for any variations in ‘Hopf-parameters’, but instead assume the Eddington approximation 3Prad = Erad = arT4 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 so-called ‘iron-opacity peak’ at T ~ 2 × 105 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 × 105 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 above2. The input parameters Teff, 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 above-mentioned iron-opacity peak, there is a clear density-inversion. 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 envelope-models 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 drop-off in density.

thumbnail Fig. 3

Initial 1D gas density (in g cm−3) structure as a function of modified radius x0 = 1 − R0/r, with R0 as the lower boundary radius. The cyan curve is corresponding FASTWIND density structure as described in Sect. 2.3. The red-dashed vertical line is the optical photosphere R. The shaded-purple region is the instability region where we provide the initial perturbations (see zoomed-in part).

2.4 Numerical grid

The lower boundary of our simulation is situated inside the stellar envelope, at a fixed R0 given by the initial condition calculation. The outer boundary extends well into the supersonic regions of the outflow, we set it here to r = 3R0. 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 MPI-AMRVAC is capable of adaptive mesh refinement (AMR), in our simulations, we use a fixed-radius refinement since otherwise, we run into methodological problems with small-scale ‘clumped’ regions in particularly the wind outflow. At our lowest resolution (Level 1), the numerical domain is 32 cells in the x-direction and 256 cells in the y-direction covering 0.2R0 and 3R0, respectively. The y-direction is here identified with the radial coordinate r, applying the spherical correction-terms discussed in Moens et al. (2022b,a), and the x-direction 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 R0 in length. This highest resolution is set from the lower boundary R0 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: Hptot=ciso,gas 2gβp,$H_{\rm{p}}^{{\rm{tot}}} = {{c_{{\rm{iso,gas }}}^2} \over {g{\beta _{\rm{p}}}}}$(24)

where β = pg/Ptot = pg/ (pg + Prad). Figure 4 confirms a posteriori that we do indeed resolve several total pressure scale heights in layers deeper than and around R.

Table 1

Fundamental parameters for the 〈2D〉 O star models studied in this paper.

thumbnail Fig. 4

Number of total pressure height scale Hptot $H_{\rm{p}}^{{\rm{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 sub-surface convec-tive region. This region may be identified following the linear analysis in Blaes & Socrates (2003), where instabilities in an optically thick radiation-pressure dominated gas were studied. The instability criterion for convection then becomes: γgpg4(γg1)Erad Ng2+13Nr2<0,${{{\gamma _{\rm{g}}}{p_{\rm{g}}}} \over {4\left( {{\gamma _{\rm{g}}} - 1} \right){E_{{\rm{rad }}}}}}N_{\rm{g}}^2 + {1 \over 3}N_{\rm{r}}^2 < 0,$(25)

where Ng2:=fgrav (1ρcs,gas2pglnρ),Nr2:=fgrav (13ρcs,rad2Eradlnρ).$\eqalign{ & N_{\rm{g}}^2: = - {f_{{\rm{grav }}}} \cdot \left( {{1 \over {\rho c_{{\rm{s}},{\rm{gas}}}^2}}\nabla {p_{\rm{g}}} - \nabla \ln \rho } \right), \cr & N_{\rm{r}}^2: = - {{\bf{f}}_{{\rm{grav }}}} \cdot \left( {{1 \over {3\rho c_{{\rm{s}},{\rm{rad}}}^2}}\nabla {E_{{\rm{rad}}}} - \nabla \ln \rho } \right). \cr} $

These are the Brunt-Vaisala 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. cs,rad and cs,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 purple-shaded parts.

In this convective region, we add small perturbations (max 20%) to the initial density profile (see the zoomed-in part in Fig. 3). Such perturbations have a plane wave dependence on r, namely, δρ~Aexp{ikr},$\delta \rho \~A\exp \{ {\bf{ik}} \cdot {\bf{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, kx[ 2πHptot ,2πLcell,x ],ky[ 2πHptot ,2πLcell,y ].${k_x} \in \left[ {{{2\pi } \over {H_{\rm{p}}^{{\rm{tot }}}}},{{2\pi } \over {{L_{{\rm{cell}},x}}}}} \right],{k_y} \in \left[ {{{2\pi } \over {H_{\rm{p}}^{{\rm{tot }}}}},{{2\pi } \over {{L_{{\rm{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 (kx,c, ky,c), (kx,c, − ky,c), (−kx,c, ky,c) and (kx,c, ky,c), respectively, with kx,c=(kx,min+kx,max)/2,ky,c=(ky,min+ky,max)/2.${k_{x,c}} = \left( {{k_{x,\min }} + {k_{x,\max }}} \right)/2,{k_{y,c}} = \left( {{k_{y,\min }} + {k_{y,\max }}} \right)/2.$(28)

Finally, perturbations are calculated by Fourier transforming N(k), so that δρ=0.2ρexp{ikr}N(k)dk$\delta \rho = 0.2\rho \int {\exp } \{ i{\bf{k}} \cdot {\bf{r}}\} N(k){\rm{d}}{\bf{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 Erad by a fixed input radiative luminosity at the bottom boundary, using a finite difference formulation of Eq. (6) to extract Erad from its gradient.

The gas energy eg is then set by assuming Trad = Tg 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 Erad, we analytically solve Eq. (6) from the outer boundary radius ro to r → ∞ by assuming constant opacity and flux-limiter, and a diffusive flux and mass density that varies according to r−2 outside ro, yielding Erad,o=Fdiff,o ro3Do${E_{{\rm{rad}},{\rm{o}}}} = {{{F_{{\rm{diff,o }}}}{r_{\rm{o}}}} \over {3{D_{\rm{o}}}}}$ where subscripts o in the parameters signify that values are taken at the radial outer boundary.

thumbnail Fig. 5

Mass loss rates (top) and effective temperature (bottom) for the O4 model as varying with time (shown here in 103 seconds). The black curve displays the laterally averaged quantity varying across time written as 〈2DL, and the black-dashed line is the lateral-temporal 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, td, as: td,a=Hptotcs,${t_{{\rm{d}},{\rm{a}}}} = {{H_{\rm{p}}^{{\rm{tot}}}} \over {{c_{\rm{s}}}}},$`(30)

where cs is the total adiabatic sound speed in the radiating gas as such, cs2=cs, rad 2+cs,gas2,$c_{\rm{s}}^2 = c_{{\rm{s}},{\rm{ rad }}}^2 + c_{{\rm{s}},{\rm{gas}}}^2,$(31)

for adiabatic radiation and gas sound speeds of cs,gas2=5pg3ρ,$c_{{\rm{s}},{\rm{gas}}}^2 = {{5{p_{\rm{g}}}} \over {3\rho }},$(32) cs,rad2=4Prad3ρ.$c_{s,{\rm{rad}}}^2 = {{4{P_{{\rm{rad}}}}} \over {3\rho }}.$(33)

On the other hand, the dynamics of the outflowing wind, starting around an ‘average’ sonic point Rsonicr [vr > ciso,gas], which in our simulations vary in space and time depending on the local conditions, is typically better characterised by, td,w=Rv,${t_{{\rm{d}},{\rm{w}}}} = {{{R_ \star }} \over 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, Hptot~0.03R$H_{\rm{p}}^{{\rm{tot}}}\~0.03{R_ \star }$, and cs ~ 150 km s−1, yielding td,a ~ 1000 s and td,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 radiation-dominated O star atmospheres by tth ~Erad Hptot /Fdiff ~103.4s${t_{{\rm{th }}}}\~{E_{{\rm{rad }}}}H_{\rm{p}}^{{\rm{tot }}}/{F_{{\rm{diff }}}}\~{10^{3. \ldots 4}}{\rm{s}}$. Alternatively, a thermal timescale may be estimated for the entire atmosphere following Grassitelli et al. (2016); Moens et al. (2022b), tth = GM Matm/ (RL), where Matm represents the mass of the stellar atmosphere. Integrating over the average radial density profile of our simulations from the iron-bump to the outer boundary, this gives typical numbers ≈ 40 000 s ≈ 4td,w. This estimate yields higher characteristic numbers than td,w, and in practice, we thus also examine whether our simulations have reached an approximate energetic steady-state 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, Teff, 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 〈Teff〉 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 simulation-times time-fluctuations in 〈〉 and 〈Teff〉 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.2R0 of the star in one of the two lateral directions, namely, only a fraction of the overall stellar and wind volume. As such, the time-variance displayed in this figure will notably not be representative of the predicted time-variability for observations of these global stellar and wind parameters. Specifically, 〈Teff〉 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 〈Teff〉, 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.

thumbnail 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 103 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 R0 and in the radial direction to 1.5 R0, with R0 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 iron-opacity peak region, which also moves spatially inward in the simulation because of the readjusting temperature and density scales. Simultaneously, a self-consistent line-driven wind outflow is initiated from the regions slightly above the optical photosphere at around 〈Rsonic〉. 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 td.w (or, ~90 ks) the simulation starts to reach a more dynamically relaxed state3.

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 time-averages of the total model luminosity 〈LtotL along the radius, where Ltot 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 time-averages of 〈Ltot〉 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 time-average of 4 snapshots within ~td,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 ~100td,w. As can be seen in the plot, the long time-average 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 time-snapshots from more advanced and dynamically relaxed stages of the O4 simulation. This shows how the Γ < 1 sub-Eddington layers close to the lower boundary have become quasi-stable, 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 super-Eddington regions (below the photosphere) with Γ > 1. In accordance with previous simulations of this opacity-driven unstable zone (Jiang et al. 2015), the net result is a highly structured and turbulent sub-surface 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 line-driven layers around and above the photosphere. This complex interplay creates a situation wherein the line-driven 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 line-opacity in the wind launching region. This leads to further structure formation and prevents the development of a stationary line-driven outflow as seen in 1D O star simulations with a similar Sobolev-based κ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 re-run the O4 simulation but starting now instead from an essentially hydrostatic envelope (achieved in practice by assuming a negligible mass-loss 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 line-driven 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 tw,d starting from the beginning of the simulation. Additionally, the figure also displays the initial conditions and a longer time-average taken after the wind has dynamically relaxed. The top panel shows the new ‘O4*’ model starting from the effectively no-wind condition (note the near-zero 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 line-driven 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, line-driving 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, line-driving 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 line-driven winds in our O star simulations are quite robust against specific details in the initial wind conditions.

thumbnail Fig. 7

Similar to Fig. 6, but now at times after the models have dynamically relaxed from their initial conditions. The black-dashed line is the averaged 2D optical photosphere 〈R〉.

thumbnail Fig. 8

Running averaged total luminosity 〈Ltot〉 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 time-average of the total luminosity.

3.3 Model O4

We go on to examine in greater detail our prototypical model O4 for an early O-supergiant in the Galaxy. As discussed above, after the initial relaxation, the O4 simulation instigates a turbulent deeper atmosphere driven by the iron-opacity peak, with line-driving 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 sub-surface iron-opacity 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 WR-stars 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 (vr, ρ, 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 over-plotted red curve is then the laterally and time ‘averaged’ quantity. As evident from Fig. 10, at the iron-opacity peak 〈Γ〉 indeed grazes unity before again decreasing as we move outwards in the atmosphere toward R. Since there is no net-outflow from these sub-surface 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 plateau-like feature. As 〈Γ〉 then decreases,〈ρ〉 starts to decrease more rapidly again until line-driving 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 line-driven supersonic outflow. This general behaviour is then also reflected in the 〈vr〉 profile. Namely, around the iron-opacity peak 〈vr〉 first increases slightly to around ~20 km s−1. The following decrease in 〈Γ〉 then produces a region with a net-infall of material before line-driving takes over above the average photosphere. As such, the net effect is that at R we interestingly enough observe a negative 〈vr〉 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 〈vdısp〉 indeed is close to zero in the deepest layers of our simulations but then increases rapidly once instabilities connected to the iron-opacity peak become effective. The moving gas particles in this region then shoot into the cooler upper atmosphere, interacting with the line-driven outflow launched from the turbulent surface layers.

For the O4-model, 〈vdı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 〈vdısp〉 values in our simulations are above the ‘gas’ sound speed, cs,gas, they are still below the ‘total’ sound speed, cs, 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 high-density clumps and yellow represents the under-dense structures. Although the under-dense 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 under-dense structures move at much higher velocities than their dense counterparts. As discussed above, line-driving is the main driver of the wind outflow here. However, since line-opacity and gas density essentially are inversely proportional to each other (see above), under-dense structures experience much stronger line-driving than dense clumps, resulting in significantly higher velocities for the former. The upshot is a strong anti-correlation 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 WR-stars). By contrast, below R* the gas is quite dense, line-driving 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.

thumbnail 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 dash-dotted lines are the initial velocity profiles. The green dash-dotted line is the final average for the O4* model and the orange dashed line is for the O4 model. See text.

thumbnail 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 lateral-temporal averages. The vertical orange-dashed line is 〈R〉, whereas the green dash-dotted vertical lines in the top panel are 〈Rsonic〉. 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.

thumbnail 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〉/Ledd, resulting in a higher 〈Teff〉, and one (O8) with a lower mass and luminosity, thus resulting in a lower Γe and 〈Teff〉. 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 iron-opacity peak, and line-driven 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 iron-opacity 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 sub-surface layers. But since the atmosphere there still is too dense for efficient line-driving, 〈Γ〉 drops again once the iron-opacity peak is overcome so that a deep sub-surface net-outflow is not launched, mimicking the O4 model. We note, however, that for this very luminous model, it is somewhat ‘easier’ for line-driving to become efficient, meaning a lower line-opacity force-multiplier boost is needed to overcome the effective gravity and launch an outflow. As such, in the O2 simulation, the average sonic point 〈Rsonic〉 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 line-driven winds (O stars, hot sub-dwarfs) to higher-luminosity objects with optically thick winds (‘slash’ or WNh-stars, classical WR-stars).

The variability in the sub-surface 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 〈vdı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.

thumbnail 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 under-dense regions and blue representing over-dense clumps. The red dash-dotted line indicates the local escape velocity, the black dashed vertical line is the average photosphere 〈R〉 and the black dash-dotted vertical line is the average sonic point 〈Rsonic〉. The zoomed-in part highlights the regions below the optical photosphere.

4 Comparison to 1D models

The previous section demonstrated the development of a very turbulent, time-dependent 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 R0. 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 re-ran the 1D model with the same stellar parameters (L, M, and R0) 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 Teff. Specifically, for the 〈2D〉 and 1D O4 models we obtain R = 16.7R and 18.7R; Teff = 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.

thumbnail 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 − R0/r, with R0 as the lower boundary radius. The ‘red dash-dotted’ line is the 〈2D〉 photosphere and the ‘blue dash-dotted’ is the 1D photosphere. The cyan dash-dotted 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 74. This treats energy transport within the standard mixing-length theory (MLT) framework for non-adiabatic convection, accounting for the effects of radiation pressure and cooling. The free input parameter, αMLT = /Hp, is the standard MLT parameter for mixing length, , and ‘total’ pressure scale height, Hp = |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 Pturb =ρvturb 2${P_{{\rm{turb }}}} = \rho v_{{\rm{turb }}}^2$ in the hydrostatic Eq. (15), where vturb is an input turbulent velocity parameter, such that now dP/dm = GM*/r2 for total pressure P = pg + Prad + Pturb. The connection of the sub-surface 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 Teff 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 iron-opacity peak, and a steeper density-profile around the photosphere.

First adding convective energy transport by setting αMLT = 1.0, but still keeping vturb = 0 km s−1, gets rid of the density-inversion also in the 1D structure. However, around the photosphere (above the convectively unstable region) the scale-height is still too small, producing the same steep density-profile.

In the next step, we first add a constant vturb ≃ 90 km s−1 throughout the atmosphere in our 1D model. This then indeed produces a shallower density-slope 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 vturb such that at the photosphere we have a maximum vturb ≃ 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 sub-surface layers of the 2D simulations may, in principle, be captured by these modifications of present-day 1D models. More quantitative comparisons and benchmarks will, however, require more in-depth 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 FASTWIND-like 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ρturb Hρrad (1+vturb 2/ag, iso 2).$H_\rho ^{{\rm{turb }}} \approx H_\rho ^{{\rm{rad }}}\left( {1 + v_{{\rm{turb }}}^2/a_{{\rm{g}},{\rm{ iso }}}^2} \right).$ For vturb ~ 90 km s−1 and αg,iso ~ 30 km s−1 (appropriate for the pho-tospheric layers) then, we get thus Hρturb 10Hρrad$H_\rho ^{{\rm{turb }}} \approx 10H_\rho ^{{\rm{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 vturb 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 radiation-hydrodynamic 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) r(ρvr2+pg)=ρGM*r2+ρκFdiff c,${\partial \over {\partial r}}\left( {\rho v_r^2 + {p_{\rm{g}}}} \right) = - \rho {{G{M_*}} \over {{r^2}}} + \rho {{\kappa {F_{{\rm{diff }}}}} \over c},$(35) vr(pg+Prad +eg+Erad +ρvr22ρGM*r)+Fdiff =Ftot =Ltot 4πr2,${v_r}\left( {{p_{\rm{g}}} + {P_{{\rm{rad }}}} + {e_{\rm{g}}} + {E_{{\rm{rad }}}} + {{\rho v_r^2} \over 2} - \rho {{G{M_*}} \over r}} \right) + {F_{{\rm{diff }}}} = {F_{{\rm{tot }}}} = {{{L_{{\rm{tot }}}}} \over {4\pi {r^2}}},$(36)

where Ltot = 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 3Prad = Erad, which is a good assumption in the deep optically thick layers, the stationary energy equation above may be equivalently written as: M˙(hg+hr+v22GM*r)+Ldiff =Ltot ,$\dot M\left( {{h_{\rm{g}}} + {h_{\rm{r}}} + {{{v^2}} \over 2} - {{G{M_*}} \over r}} \right) + {L_{{\rm{diff }}}} = {L_{{\rm{tot }}}},$(37)

for gas and radiation enthalpies hg = (5/2)pg/ρ and hr = 4Prad/ρ, respectively, and diffusive luminosity Ldiff = 4πr2 Fdiff (e.g. Owocki et al. 2017). These stationary relations invite us to identify the following suitable lateral- and time-averages as corresponding to ‘turbulent pressure’ and ‘convective flux’ in our simulations: Pturb ¯= ρvr2 ,$\overline {{P_{{\rm{turb }}}}} = \left\langle {\rho v_r^2} \right\rangle ,$(38) Fconv ¯= vr(pg+Prad +eg+Erad ) ,$\overline {{F_{{\rm{conv }}}}} = \left\langle {{v_r}\left( {{p_{\rm{g}}} + {P_{{\rm{rad }}}} + {e_{\rm{g}}} + {E_{{\rm{rad }}}}} \right)} \right\rangle ,$(39)

with associated turbulent and convective velocities vturb ¯=Pturb ¯ρ,$\overline {{v_{{\rm{turb }}}}} = \sqrt {{{\overline {{P_{{\rm{turb }}}}} } \over {\langle \rho \rangle }}} ,$(40) vconv ¯=Fconv ¯ pg+Prad +eg+Erad  .$\overline {{v_{{\rm{conv }}}}} = {{\overline {{F_{{\rm{conv }}}}} } \over {\left\langle {{p_{\rm{g}}} + {P_{{\rm{rad }}}} + {e_{\rm{g}}} + {E_{{\rm{rad }}}}} \right\rangle }}.$(41)

That is, the turbulent velocity is here simply identified as the density-weighted average root mean square (rms) velocity (see also Goldberg et al. 2022). As mentioned above, this expression for Pturb ¯$\overline {{P_{{\rm{turb }}}}} $ essentially assumes that the turbulent velocity field is isotropic in the relevant layers (specifically for our case that ρvr2 ρvt2 )$\left. {\left\langle {\rho v_r^2} \right\rangle \approx \left\langle {\rho v_t^2} \right\rangle } \right)$, 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 sub-surface 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 vturb ¯vconv ¯$\overline {{v_{{\rm{turb }}}}} \ne \overline {{v_{{\rm{conv }}}}} {\rm{. }}$. For the O star simulations here the additional kinetic and gravitational energy fluxes are small, and Fconv ¯ $\overline {{F_{{\rm{conv }}}}} {\rm{ }}$ 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 Fconv ¯ $\overline {{F_{{\rm{conv }}}}} {\rm{ }}$ 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. Fconv ¯ $\overline {{F_{{\rm{conv }}}}} {\rm{ }}$ is further centered around the iron-opacity 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 Fconv 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 Fconv ¯ $\overline {{F_{{\rm{conv }}}}} {\rm{ }}$, but with maximum levels of Fconv ¯ $\overline {{F_{{\rm{conv }}}}} {\rm{ }}$ 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 radiation-dominated atmosphere around the iron-opacity peak. As the opacity then decreases again, the turbulent pressure declines rapidly and the gas pressure, pg, 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 vturb ¯$\overline {{v_{{\rm{turb }}}}} $ for the O4 model. The qualitative behaviour of vturb ¯$\overline {{v_{{\rm{turb }}}}} $is quite similar to the straight average vdı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 density-weighted vturb ¯$\overline {{v_{{\rm{turb }}}}} $ is somewhat lower than vdısp. This is most likely related to the influence of the line-driven wind and the density-velocity anti-correlation 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, vconv ¯$\overline {{v_{{\rm{conv }}}}} $ only reaches a peak value (in the iron-opacity 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 iron-opacity peak region. Moreover, below the optical photosphere, although the gas sound speed cSıgas is around ~30 km s−1, the total sound speed cs is above 100 km s−1. As such, both vturb ¯ and vconv ¯$\overline {{v_{{\rm{turb }}}}} {\rm{ and }}\overline {{v_{{\rm{conv }}}}} $ stay well below the local total sound speed in our simulations.

Broadly, we also observe the same trends for vturb ¯$\overline {{v_{{\rm{turb }}}}} $ in the O2 and O8 models as for the O4 simulation, with vturb ¯$\overline {{v_{{\rm{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, pg and Prad 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 pg and Prad are similar to the O4 model, although Pturb ¯$\overline {{P_{{\rm{turb }}}}} $ approaches Prad much quicker than in the O4 model.

thumbnail 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 〈Teff〉 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 dash-dotted curve is produced using FASTWIND, using the specifications outlined in Table 1 and Sect. 2.3.

thumbnail 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.

thumbnail 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 103 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 multi-dimensional simulations as compared to 1D models based on their individual averages. That is, generally 〈ρκFdiff〉≠ 〈ρ〉〈κ〉〈Fdiff〉 (Shaviv 1998, see also discussion in Jiang 2023). Figure 17 displays this difference for our O4 simulation, illustrating that while the quasi-steady 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 quasi-static 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).

thumbnail Fig. 17

Flux-porosity 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 metallicity-scaling based on Björklund et al. (2021). Although this re-calibration 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 O-supergiants 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 line-driving 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 input-parameter 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 self-consistent emergent properties of the models.

thumbnail 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 dot-dashed 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 sub-surface iron opacity peak region of our model O stars. This means our simulations are generally located in the region of so-called ‘inefficient’ convection, which we may interpret through some simple scaling relations. For the radiation-dominated atmospheres studied here, an upper limit for convective flux would scale as (e.g. Gräfener et al. 2012) FConv ~ CsErad, so that FConv/Fdiff ~ (cs/c)(T/Teff,0)4 ~ (cs /c)τ, where Teff,04=L/(σ4πR02)$T_{{\rm{eff}},0}^4 = {L_ \star }/\left( {\sigma 4\pi R_0^2} \right)$ is an effective (flux) temperature evaluated below the iron-opacity peak at ~R0 and where the generic scaling τ ~ (T/Teff,0)4 has been applied further. Alternatively, we may view this scaling through the characteristic dynamical and radiative diffusion timescales for the region, td,a~Hptot/cs${t_{{\rm{d}},{\rm{a}}}}\~H_{\rm{p}}^{{\rm{tot}}}/{c_{\rm{s}}}$ and tdiff ~(τ/c)Hptot${t_{{\rm{diff }}}}\~(\tau /c)H_{\rm{p}}^{{\rm{tot}}}$, such that tdiff/td,a ~ (cs/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 quasi-stable radiative zone below. Setting tdiff = td,a then also gives us the ‘critical optical depth’ τc = c/cs as defined and used in Jiang et al. (2015) to distinguish between efficient and non-efficient 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 Fconv/Fdiff ~ tdiff/td,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 WR-stars, 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 LBV-like massive stars, the ratio of convective to diffusive fluxes may be significantly higher for the iron-opacity peak zone (Jiang et al. 2015, 2018).

Since the majority of the radiative flux thus cannot be carried by convective transport through the iron-opacity 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 line-broadening’ to match the line profiles computed by means of 1D atmospheric models to high-resolution 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ón-Díaz et al. 2017).

Additionally, 1D models often apply a moderate amount of (also ad-hoc) photospheric ‘microturbulence’ (~ 10–20 km s−1 in standard FASTWIND-based 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 macro-turbulence is usually motivated by means of optical depth arguments, as essentially borrowed from spectral analysis of cool, low-mass stars (see detailed overview in Gray 2008); in practice, this just means that microturbulence is explicitly included in the line-profile functions for computing number densities and synthetic spectra whereas macroturbu-lence is only applied afterward as a convolution that further broadens the spectral lines but preserves their equivalent widths (thereby also causing degeneracy-problems between rotational and macroturbulent line-broadening, see Sundqvist et al. 2013).

Overall, our simulations thus provide a natural explanation for the need for this ad hoc ‘extra line-broadening’ mechanism necessary to include in the 1D model atmosphere and spectrum synthesis codes (which in a sense mirrors results found already for low-mass sun-like 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 (so-called ‘P Cygni’ lines), which typically require turbulent velocities (usually included in the line-profile 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 line-deshadowing 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 Sobolev-based (Moens et al. 2022a) and LDI-based (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 ‘two-component 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 small-scale 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, small-scale clumps. Indeed, inspection of the probability colour maps in this paper really does not seem to find much basic support for a two-component medium ansatz; rather the most likely wind-density is close to the mean, with a large dispersion around it, and further with a strong anti-correlation 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 present-day 1D model atmosphere codes are in need of a careful re-investigation. Owocki & Sundqvist (2018) recently made a first such attempt to relax the standard two-component 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 spectro-scopic derivation of fundamental stellar parameters (Teff, 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 so-called ‘mass-discrepancy 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 Pturb ¯$\overline {{P_{{\rm{turb }}}}} $ in the hydrostatic equation, and, moreover, a moderate Fconv ¯ $\overline {{F_{{\rm{conv }}}}} {\rm{ }}$ 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 multi-dimensional 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 energy-weighted mean opacities are the same. In particular, this may affect q˙$\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 q˙$\dot q$ in the wind-parts 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 q˙$\dot q$ in such O star winds. As such, we have very recently attempted to extend our method for computing κline to (approximate) non-LTE 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 double-count line opacities in an ‘intermediate’ region wherein the Sobolev-based κ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 mass-loss rates overall agree well with those predicted by detailed 1D models based on full non-LTE and multi-frequency co-moving 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 mass-loss 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 frequency-range (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 Sobolev-based line-driving 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 line-driven 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) non-local integrations to obtain the corresponding optical depths. Typically, all LDI simulations carried out thus far have been based on some variant of the escape-integral 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 line-driving stemming from thousands of spectral lines within such an approach remains a challenging question for multi-dimensional simulations though.

Finally, in this study, we have only performed 2D time-dependent 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 multi-dimensional, time-dependent unified atmosphere with wind models for O stars. To this end, we use a ‘box-in-a-star’ approach (including spherical correction terms), while solving the RHD equations using an analytical flux-limiting closure approximation in the radiation energy equation, and accounting properly for line-driving 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 sub-surface iron-opacity 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 iron-opacity 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 vturb 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 line-driving becomes efficient; this is due to the inverse relation between density and line-driving opacity, the outflowing wind parts also display a clear anti-correlation 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 line-driving into our simulations.

A key follow-up 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 present-day 1D codes might have to be calibrated to account for basic multi-dimensional 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 energy-weighted mean opacities, using an approximate non-LTE 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: B-supergiants and LBVs.

Acknowledgements

The computational resources used for this work were provided by Vlaams Supercomputer Centrum (VSC) funded by the Research Foundation-Flanders (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 easy-to-follow (albeit in hand-written German, of course) koch rezept for MLT-based 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 KUL-EQUATION 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 mean-free-path between interactions with different lines is mfp = 1/(κeffρ), where κeff is the effective opacity. Here lmfp = Δl/P for Δl the physical mean distance between two lines and P=1eτs$P = 1 - {e^{ - {\tau _s}}}$ the line transition probability. By means of the Doppler shift and Sobolev line optical depth, we obtain: τs=αlcv0| dvldl |1,${\tau _s} = {{{\alpha _l}c} \over {{v_0}}}{\left| {{{d{v_l}} \over {dl}}} \right|^{ - 1}},$(A.1)

and the effective line opacity is thus: κeff=1ραlΔv1eτsτs,${\kappa _{{\rm{eff}}}} = {1 \over \rho }{{{\alpha _l}} \over {\Delta v}}{{1 - {{\rm{e}}^{ - {\tau _s}}}} \over {{\tau _s}}},$(A.2)

where αl is the frequency-integrated 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 κeff=1ρilines in ΔvαliΔv(1eτsiτsi)+κTH${\kappa _{{\rm{eff}}}} = {1 \over \rho }\sum\limits_i^{{\rm{lines in }}\Delta v} {{{\alpha _l^i} \over {\Delta v}}} \left( {{{1 - {{\rm{e}}^{ - \tau _s^i}}} \over {\tau _s^i}}} \right) + {\kappa _{{\rm{TH}}}}$(A.3) =dvldl1cρilines in ΔvviΔv(1eτsi)+κTH.$ = {{d{v_l}} \over {dl}}{1 \over {c\rho }}\sum\limits_i^{{\rm{lines in }}\Delta v} {{{{v^i}} \over {\Delta v}}} \left( {1 - {{\rm{e}}^{ - \tau _s^i}}} \right) + {\kappa _{{\rm{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 /qi (τ = 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, κeff=1ρilines in ΔvαliΔv(1eτsiτsi)+κRoss ${\kappa _{{\rm{eff}}}} = {1 \over \rho }\sum\limits_i^{{\rm{lines in }}\Delta v} {{{\alpha _l^i} \over {\Delta v}}} \left( {{{1 - {{\rm{e}}^{ - \tau _s^i}}} \over {\tau _s^i}}} \right) + {\kappa ^{{\rm{Ross }}}}$(A.5)

Connection to hybrid opacity model. The flux-weighted opacity of the above expansion opacity is: κeffFvΔvF${{{\kappa _{{\rm{eff}}}}{F_v}\Delta v} \over F}$(A.6)

when taken over the complete frequency range, this becomes: κ=κ0iall lines qiwi(1eτsiτsi)+κRoss $\kappa = {\kappa _0}\sum\limits_i^{{\rm{all lines }}} {{q_i}} {w_i}\left( {{{1 - {{\rm{e}}^{ - \tau _s^i}}} \over {\tau _s^i}}} \right) + {\kappa ^{{\rm{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 τ (∝ ρ/dvl/dl).

References

  1. Abbott, D. C., & Lucy, L. B. 1985, ApJ, 288, 679 [NASA ADS] [CrossRef] [Google Scholar]
  2. Asplund, M., Nordlund, A., Trampedach, R., Allende Prieto, C., & Stein, R. F. 2000, A&A, 359, 729 [NASA ADS] [Google Scholar]
  3. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  4. Björklund, R., Sundqvist, J. O., Puls, J., & Najarro, F. 2021, A&A, 648, A36 [EDP Sciences] [Google Scholar]
  5. Blaes, O., & Socrates, A. 2003, ApJ, 596, 509 [CrossRef] [Google Scholar]
  6. Brands, S. A., de Koter, A., Bestenlehner, J. M., et al. 2022, A&A, 663, A36 [Google Scholar]
  7. Cantiello, M., Langer, N., Brott, I., et al. 2009, A&A, 499, 279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Castor, J. I. 2004, Radiation Hydrodynamics (Cambridge University Press) [CrossRef] [Google Scholar]
  9. Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [Google Scholar]
  10. Conti, P. S., & Ebbets, D. 1977, ApJ, 213, 438 [CrossRef] [Google Scholar]
  11. Driessen, F., Sundqvist, J. O., & Kee, N. D. 2022a, PhD Thesis, KU Leuven, Belgium [Google Scholar]
  12. Driessen, F. A., Sundqvist, J. O., & Dagore, A. 2022b, A&A, 663, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Eversberg, T., Lépine, S., & Moffat, A. F. J. 1998, ApJ, 494, 799 [NASA ADS] [CrossRef] [Google Scholar]
  14. Freytag, B., Steffen, M., Ludwig, H. G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
  15. Friend, D. B., & Abbott, D. C. 1986, ApJ, 311, 701 [Google Scholar]
  16. Friend, D. B., & Castor, J. I. 1983, ApJ, 272, 259 [NASA ADS] [CrossRef] [Google Scholar]
  17. Gabler, R., Gabler, A., Kudritzki, R. P., Puls, J., & Pauldrach, A. 1989, A&A, 226, 162 [NASA ADS] [Google Scholar]
  18. Gayley, K. G. 1995, ApJ, 454, 410 [Google Scholar]
  19. Gayley, K. G., & Owocki, S. P. 1995, ApJ, 446, 801 [NASA ADS] [CrossRef] [Google Scholar]
  20. Glatzel, W. 1994, MNRAS, 271, 66 [NASA ADS] [CrossRef] [Google Scholar]
  21. Goldberg, J. A., Jiang, Y.-F., & Bildsten, L. 2022, ApJ, 929, 156 [NASA ADS] [CrossRef] [Google Scholar]
  22. Gräfener, G., Koesterke, L., & Hamann, W. R. 2002, A&A, 387, 244 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Gräfener, G., Owocki, S. P., & Vink, J. S. 2012, A&A, 538, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Grassitelli, L., Chené, A. N., Sanyal, D., et al. 2016, A&A, 590, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Grassitelli, L., Langer, N., Grin, N. J., et al. 2018, A&A, 614, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Gray, D. F. 2008, The Observation and Analysis of Stellar Photospheres (Cambridge University Press) [Google Scholar]
  27. Grevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, eds. N. Prantzos, E. Vangioni-Flam, & M. Casse (Cambridge University Press), 15 [Google Scholar]
  28. Hamann, W. R., & Gräfener, G. 2004, A&A, 427, 697 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
  30. Haser, S. M., Lennon, D. J., Kudritzki, R. P., et al. 1995, A&A, 295, 136 [NASA ADS] [Google Scholar]
  31. Hawcroft, C., Sana, H., Mahy, L., et al. 2021, A&A, 655, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Hawcroft, C., Sana, H., Mahy, L., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202245588 [Google Scholar]
  33. Hearn, A. G. 1972, A&A, 19, 417 [NASA ADS] [Google Scholar]
  34. Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2020, A&A, 633, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Herrero, A., Kudritzki, R. P., Vilchez, J. M., et al. 1992, A&A, 261, 209 [NASA ADS] [Google Scholar]
  36. Hillier, D. J. 1991, A&A, 247, 455 [NASA ADS] [Google Scholar]
  37. 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]
  38. Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
  39. Howarth, I. D., Siebert, K. W., Hussain, G. A. J., & Prinja, R. K. 1997, MNRAS, 284, 265 [NASA ADS] [CrossRef] [Google Scholar]
  40. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  41. Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
  42. Jiang, Y.-F. 2023, Galaxies, 11, 105 [NASA ADS] [CrossRef] [Google Scholar]
  43. Jiang, Y.-F., Cantiello, M., Bildsten, L., Quataert, E., & Blaes, O. 2015, ApJ, 813, 74 [Google Scholar]
  44. Jiang, Y.-F., Cantiello, M., Bildsten, L., et al. 2018, Nature, 561, 498 [NASA ADS] [CrossRef] [Google Scholar]
  45. Karp, A. H., Lasher, G., Chan, K. L., & Salpeter, E. E. 1977, ApJ, 214, 161 [NASA ADS] [CrossRef] [Google Scholar]
  46. Kee, N. D., Owocki, S., & Sundqvist, J. O. 2016, MNRAS, 458, 2323 [Google Scholar]
  47. Keppens, R., Teunissen, J., Xia, C., & Porth, O. 2021, Comput. Math. Appl., 81, 316 [Google Scholar]
  48. Keppens, R., Popescu Braileanu, B., Zhou, Y., et al. 2023, A&A, 673, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Kippenhahn, R., Weigert, A., & Weiss, A. 2013, Stellar Structure and Evolution (Springer) [Google Scholar]
  50. Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Krticka, J., & Kubát, J. 2017, A&A, 606, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Lucy, L. B. 1971, ApJ, 163, 95 [Google Scholar]
  53. Lucy, L. B., & Abbott, D. C. 1993, ApJ, 405, 738 [Google Scholar]
  54. Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (Oxford University Press) [Google Scholar]
  55. Moens, N., Poniatowski, L. G., Hennicker, L., et al. 2022a, A&A, 665, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Moens, N., Sundqvist, J. O., El Mellah, I., et al. 2022b, A&A, 657, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Oskinova, L. M., Hamann, W. R., & Feldmeier, A. 2007, A&A, 476, 1331 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. 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]
  59. Owocki, S. P., & Puls, J. 1996, ApJ, 462, 894 [Google Scholar]
  60. Owocki, S. P., & Rybicki, G. B. 1984, ApJ, 284, 337 [Google Scholar]
  61. Owocki, S. P., & Sundqvist, J. O. 2018, MNRAS, 475, 814 [NASA ADS] [CrossRef] [Google Scholar]
  62. Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
  63. Owocki, S. P., Gayley, K. G., & Shaviv, N. J. 2004, ApJ, 616, 525 [NASA ADS] [CrossRef] [Google Scholar]
  64. Owocki, S. P., Townsend, R. H. D., & Quataert, E. 2017, MNRAS, 472, 3749 [Google Scholar]
  65. Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] [Google Scholar]
  66. 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]
  67. Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Petrenz, P., & Puls, J. 2000, in ASP Conf. Ser., 214, IAU Colloq. 175: The Be Phenomenon in Early-Type Stars, eds. M. A. Smith, H. F. Henrichs, & J. Fabregat, 626 [NASA ADS] [Google Scholar]
  69. Petrovic, J., Pols, O., & Langer, N. 2006, A&A, 450, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Poniatowski, L. G., Sundqvist, J. O., Kee, N. D., et al. 2021, A&A, 647, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Poniatowski, L. G., Kee, N. D., Sundqvist, J. O., et al. 2022, A&A, 667, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Puls, J., Markova, N., Scuderi, S., et al. 2006, A&A, 454, 625 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Puls, J., Najarro, F., Sundqvist, J. O., & Sen, K. 2020, A&A, 642, A172 [EDP Sciences] [Google Scholar]
  76. Rogers, F. J., & Iglesias, C. A. 1992, ApJS, 79, 507 [CrossRef] [Google Scholar]
  77. Rybicki, G. B., Owocki, S. P., & Castor, J. I. 1990, ApJ, 349, 274 [NASA ADS] [CrossRef] [Google Scholar]
  78. Sander, A., Hamann, W.-R., & Todt, H. 2012, A&A, 540, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. 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]
  80. Santolaya-Rey, A. E., Puls, J., & Herrero, A. 1997, A&A, 323, 488 [NASA ADS] [Google Scholar]
  81. Sanyal, D., Grassitelli, L., Langer, N., & Bestenlehner, J. M. 2015, A&A, 580, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Schultz, W. C., Bildsten, L., & Jiang, Y.-F. 2020, ApJ, 902, 67 [NASA ADS] [CrossRef] [Google Scholar]
  83. Schultz, W. C., Bildsten, L., & Jiang, Y.-F. 2022, ApJ, 924, L11 [NASA ADS] [CrossRef] [Google Scholar]
  84. Schwarzschild, M. 1958, Structure and Evolution of the Stars (Princeton University Press) [Google Scholar]
  85. Shaviv, N. J. 1998, ApJ, 494, L193 [NASA ADS] [CrossRef] [Google Scholar]
  86. Simón-Díaz, S., Godart, M., Castro, N., et al. 2017, A&A, 597, A22 [CrossRef] [EDP Sciences] [Google Scholar]
  87. Sobolev, V. V. 1960, Moving Envelopes of Stars (Harvard University Press) [Google Scholar]
  88. Stothers, R. B., & Chin, C.-W. 1993, ApJ, 412, 294 [NASA ADS] [CrossRef] [Google Scholar]
  89. Sundqvist, J. O., & Puls, J. 2018, A&A, 619, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Sundqvist, J. O., Puls, J., & Feldmeier, A. 2010, A&A, 510, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Sundqvist, J. O., Simón-Díaz, S., Puls, J., & Markova, N. 2013, A&A, 559, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Sundqvist, J. O., Owocki, S. P., & Puls, J. 2018, A&A, 611, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Sundqvist, J. O., Björklund, R., Puls, J., & Najarro, F. 2019, A&A, 632, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Š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]
  95. Teunissen, J., & Keppens, R. 2019, Comput. Phys. Commun., 245, 106866 [Google Scholar]
  96. Turner, N. J., & Stone, J. M. 2001, ApJS, 135, 95 [NASA ADS] [CrossRef] [Google Scholar]
  97. ud Doula, A., & Owocki, S. P. 2002, ApJ, 576, 413 [NASA ADS] [CrossRef] [Google Scholar]
  98. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  99. Xia, C., Teunissen, J., Mellah, I. E., Chané, E., & Keppens, R. 2018, ApJS, 234, 30 [NASA ADS] [CrossRef] [Google Scholar]

1

Here, this is the frequency-integrated 0th angular moment of the time-dependent radiative transfer equation in the comoving frame.

2

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.

3

By which we mean the simulations have relaxed in a statistical and global sense; the atmospheres remain time-dependent and dynamically very active throughout the simulations.

4

In practice, our implementation of this formalism follows a ‘koch rezept’ stemming from a set of notes developed by Jo Puls.

All Tables

Table 1

Fundamental parameters for the 〈2D〉 O star models studied in this paper.

All Figures

thumbnail Fig. 1

Line force multiplier M(τ) as a function of the characteristic Sobolev optical depth τ. The black-dashed 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 MG(τ) in Eq. (13). The best fit for the line force parameters for this particular temperature and density are Q¯=2278$\bar Q = 2278$, Q0 = 1418, and α = 0.59.

In the text
thumbnail Fig. 2

Averaged opacity 〈κ〉 as a function of scaled radius x = 1 − R/r, with R as the 2D averaged optical photosphere (red-dashed line). In the figure, the brown curve is the Rosseland mean opacity and the purple curve is the finite-disk corrected line-opacity. The black dash-dotted curve is the total opacity using the hybrid opacity scheme.

In the text
thumbnail Fig. 3

Initial 1D gas density (in g cm−3) structure as a function of modified radius x0 = 1 − R0/r, with R0 as the lower boundary radius. The cyan curve is corresponding FASTWIND density structure as described in Sect. 2.3. The red-dashed vertical line is the optical photosphere R. The shaded-purple region is the instability region where we provide the initial perturbations (see zoomed-in part).

In the text
thumbnail Fig. 4

Number of total pressure height scale Hptot $H_{\rm{p}}^{{\rm{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
thumbnail Fig. 5

Mass loss rates (top) and effective temperature (bottom) for the O4 model as varying with time (shown here in 103 seconds). The black curve displays the laterally averaged quantity varying across time written as 〈2DL, and the black-dashed line is the lateral-temporal averaged quantity, here 〈2D〉.

In the text
thumbnail 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 103 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 R0 and in the radial direction to 1.5 R0, with R0 as the lower boundary radius of our simulation.

In the text
thumbnail Fig. 7

Similar to Fig. 6, but now at times after the models have dynamically relaxed from their initial conditions. The black-dashed line is the averaged 2D optical photosphere 〈R〉.

In the text
thumbnail Fig. 8

Running averaged total luminosity 〈Ltot〉 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 time-average of the total luminosity.

In the text
thumbnail 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 dash-dotted lines are the initial velocity profiles. The green dash-dotted line is the final average for the O4* model and the orange dashed line is for the O4 model. See text.

In the text
thumbnail 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 lateral-temporal averages. The vertical orange-dashed line is 〈R〉, whereas the green dash-dotted vertical lines in the top panel are 〈Rsonic〉. 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
thumbnail 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
thumbnail 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 under-dense regions and blue representing over-dense clumps. The red dash-dotted line indicates the local escape velocity, the black dashed vertical line is the average photosphere 〈R〉 and the black dash-dotted vertical line is the average sonic point 〈Rsonic〉. The zoomed-in part highlights the regions below the optical photosphere.

In the text
thumbnail 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 − R0/r, with R0 as the lower boundary radius. The ‘red dash-dotted’ line is the 〈2D〉 photosphere and the ‘blue dash-dotted’ is the 1D photosphere. The cyan dash-dotted curve is produced using FASTWIND as described with specifications in Sect. 2.3.

In the text
thumbnail 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 〈Teff〉 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 dash-dotted curve is produced using FASTWIND, using the specifications outlined in Table 1 and Sect. 2.3.

In the text
thumbnail 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
thumbnail 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 103 K.

In the text
thumbnail Fig. 17

Flux-porosity effect as a function of scaled radius for our O4 model.

In the text
thumbnail 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 dot-dashed 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.