© The Authors 2025

1 Introduction

Traditionally, the assumptions of a spherically symmetric, stationary one-dimensional (1D) atmosphere and wind have been made in spectroscopic studies of O stars (for example Hillier & Miller 1998; Kudritzki & Puls 2000; Puls et al. 2005; Sander et al. 2012). However, theoretical studies have long indicated that the coupled envelopes, atmospheres, and line-driven winds (Castor et al. 1975) of O stars are highly structured and variable (due to convective and radiative envelope instabilities, as well as wind instabilities, Hearn 1972; Owocki & Rybicki 1984; Blaes & Socrates 2003; Cantiello et al. 2009; Jiang et al. 2015; Van der Sijpt et al. 2025; Key et al. 2025). Other strong indications can be found in various observed phenomena; for example optical absorption lines that indicate very large ‘micro’ and/or ‘macroturbulent’ velocities (the latter even going up to >100 km/s, see for example Simón-Díaz et al. 2017), optical emission lines that show line profile variability, and ‘wind clumping’ (Conti & Ebbets 1977; Eversberg et al. 1998; Puls et al. 2006; Sundqvist et al. 2010; Simón-Díaz et al. 2017).

Recently, we performed time-dependent, two-dimensional (2D) unified simulations of O-star atmospheres with winds using a flux-limiting radiation-hydrodynamical (RHD) finite volume modelling technique, accounting properly for the effects of line driving (Debnath et al. 2024). In these simulations, opacities were computed using a hybrid approach (Poniatowski et al. 2021) that combined tabulated Rosseland mean opacities with calculations of the enhanced line opacities expected for supersonic flows. The latter were here based on the Sobolev (1960) approximation. In these simulations then, structure formation can already be found just below the iron-opacity peak (located at approximately 200 kK, Iglesias & Rogers 1996). By means of radiative acceleration exceeding gravity, local pockets of gas shoot up from these deep layers into the atmosphere above. Once in the upper atmosphere, these pockets then interact with the overlying line-driven wind outflow, creating a highly structured, turbulent atmosphere. This interaction also gives rise to large turbulent velocities in the photospheric layers of such unified atmosphere and wind simulations for O stars. The order of the turbulent velocities is ~30–100 km/s, with higher values for models with higher luminosity-to-mass ratios. This is in general good agreement with observations of photospheric macroturbulence in O stars (Simón-Díaz et al. 2017), and also overall agrees well with results from recent independent three-dimensional (3D) massive-star modelling by Schultz et al. (2023), albeit without accounting for the effects of line driving and wind outflows. In that work, the authors used the Monte-Carlo radiative transport framework SEDONA (Kasen et al. 2006) to produce absorption lines from 3D RHD simulations (which were run using ATHENA++, Jiang et al. 2015; Schultz et al. 2022).

In this paper, we present the first results from our 3D unified atmosphere and wind simulations of O stars. We focus here on the spectroscopic signatures of such simulations using 3D radiative transfer techniques to create synthetic images and spectra (building from the general code package presented in Hennicker et al. 2020, 2022).

2 3D unified atmosphere and wind simulations of O-type stars

2.1 RHD simulation method

Previously, Debnath et al. (2024) studied unified atmosphere and wind simulations of O-type stars using the RHD module from Moens et al. (2022b) of MPI-AMRVAC (Xia et al. 2018; Keppens et al. 2023) in 2D. We took their prototypical ‘O4’ model, and used the same numerical set-up to calculate new simulations by solving the time-dependent RHD equations on a finite volume in three spatial dimensions, including correction terms for spherical divergence (see Moens et al. 2022b). This means that we solved the equations of mass, momentum, and energy conservation, including the effects of the radiation field in both the total momentum (radiation force) and energy (radiative heating and cooling). Gravity was included for a fixed point stellar mass, M_⋆. The radiation field was computed from the frequency-integrated time-dependent radiation energy equation in the co-moving frame. The necessary closure relations between radiation energy density, pressure, and flux were obtained analytically with a flux-limiting procedure that preserves the optically thick and thin limits. Mean opacities include Rosseland mean opacities and line driving effects in the same way as was first outlined in Poniatowski et al. (2022), and as in previous work we assumed that the flux, energy, and Planck mean opacities were the same. For further details, including exact equations as well as specifications of initial and boundary conditions, see Moens et al. (2022a,b); Debnath et al. (2024). The radial direction extends (approximately) thrice the lower boundary radius, R₀. This R₀ is defined as the radius in which the temperature in the initial conditions reaches 500 kK, here ≃ 13.54 R_⊙. Both the tangential directions cover 0.2 R₀. Radially, the box has 128 grid points for the lowest resolution, and both tangential directions have 16 grid points. The highest level of refinement is level four, which has a resolution that is a factor of 2³ higher than the lowest resolution (first level).

2.2 Basic characteristics of early O-star 3D model

As in the 2D simulations by Debnath et al. (2024), structure starts appearing in our 3D model just below the iron-opacity peak, where the gas and radiation temperature are ~200 kK. As energy transport by convection in these atmospheres is inefficient (Debnath et al. 2024), the very sensitive dependence of the Rosseland mean opacity on temperature around the iron opacity peak yields local pockets of gas that are radiatively accelerated exceeding gravity (g_rad > g_grav). As a result, these pockets shoot up with velocities higher than the gas sound speed (however, less than the gas+radiation sound speed, see Debnath et al. 2024) into the upper and cooler atmosphere. There, Rosseland mean opacities are lower than in the deeper layers, and the gas is decelerated. While some gas parcels fall back into the deep atmosphere, others become subject to line driving and re-accelerate to launch a supersonic wind outflow. This complex interplay creates a very turbulent atmosphere characterised by large density and temperature fluctuations and high velocity dispersions. A more complete analysis of this process has been explored in Debnath et al. (2024) for a set of 2D O-star models, and Moens et al. (2022a) for a 3D stripped, hot dwarf model.

Fig. 1 shows a volume rendering of relative density in a selected part of the 3D atmosphere, taken at a snapshot well after the simulation has adjusted from the initial conditions. Relative density here means that the displayed colours are log₁₀ρ/⟨ρ⟩ where angle brackets denote lateral averaging at each radial position. The 3D rendering shows the part of the atmosphere from just below the ‘iron bump’ where inhomogeneities first occur, across the turbulent and variable photosphere, and into the out-flowing line-driven wind (up to 1.5 above the average stellar surface). The figure clearly illustrates the turbulent, filamented, and variable nature of the O-star envelope and atmosphere. This can also been seen in Fig. 2, in which we see one snapshot sliced at different laterally averaged electron scattering optical depths (and hence different heights above the inner boundary of the snapshot). The upper three panels correspond to a wind part of the snapshot, the middle three panels are around the optical photosphere and the bottom three deep down in the atmosphere. Inspecting the panels in Fig. 2, it can be seen that these three different regions have different structures, velocities, and temperature fluctuations. As noted above, in the deep atmosphere around the iron bump, opacity has a very strong temperature dependence. On the other hand, Rosseland mean opacities are also positively correlated with density, meaning that we sometimes observe blobs of high-density gas around the photosphere that are radiatively accelerated (see also Jiang et al. 2018) and thus have rather high positive velocities (middle panels). By contrast, radiative acceleration due to line driving is negatively correlated with density, implying that high density wind regions (upper panels) do not experience efficient line driving but are rather dragged along (with relatively low velocities) by the surrounding rapidly accelerating low-density wind (see also discussions and figures in Moens et al. 2022a; Debnath et al. 2024). Moreover, radiation temperature fluctuations are highly damped in the wind parts of the simulation as compared to the large fluctuations around the photosphere and in the deeper layers.

Table 1 lists fundamental average parameters of the 3D RHD model. Lateral averages are here taken over 65 snapshots each separated in time by ~500 seconds, taken well after the simulation has relaxed from initial conditions. Photospheric parameters are then computed at the average vertical position in the atmosphere corresponding to continuum ⟨τ_c⟩ = 2/3, whereas the average mass-loss rate, ⟨$\dot{M}$⟩, and maximum wind speed, ⟨v_max⟩, are taken close to the outer boundary in the outflowing parts of the simulation¹. Finally the Eddington luminosity is defined as usual, L_Edd ≡ 4πGM_⋆c/κ_e, where κ_e is the electron scattering opacity, which we also take to be the sole continuum opacity in the computation of the optical depth. Similar to Debnath et al. (2024), we computed that the characteristic dynamical timescale in the deeper atmosphere is of the order of t_d,a ~ 1000 s, while the dynamical timescale in the wind is of the order of t_d,w ~ 10 000 s.

Fig. 1 Volume rendering of relative density (displayed colours are log10ρ/⟨ρ⟩) for the 3D model atmosphere and wind. Angle brackets denote lateral averaging at each radial position.

Fig. 2 Maps of logarithm of density, radial velocity, and radiation temperature, in three selected vertical planes with laterally averaged continuum optical depth according to the legends. The lateral X and Y axes are labelled in units of R0.

3 Methodology of radiative transfer calculations

Three-dimensional radiative transfer techniques (building from Hennicker et al. 2020) are now applied to the 3D unified atmosphere and wind simulations of O stars described above, more specifically on the whole domain of the spherical models, constructed from the 3D unified atmosphere and wind simulations. The construction of these spherical models is outlined in Section 3.3.

3.1 Equations of radiative transfer

Specific intensities were obtained by solving the time-independent equation of radiative transfer:

$$\mathbf{n}\mathrm{\nabla }{I}_{\nu }={\eta }_{\nu }-{\chi }_{\nu }{I}_{\nu }={\chi }_{\nu }(S_{\nu }-I_{\nu }),$$(1)

where I_ν is the observer’s frame specific intensity at frequency ν, η_ν the emissivity, χ_ν the opacity, S_ν the source function, and n the observer’s direction. In the main part of this work we assumed that opacities and source functions could be calculated locally, using an approximate non-local thermodynamic equilibrium (aNLTE) technique following Lucy & Abbott (1993); Springmann (1997); Puls et al. (2000). This method approximately corrects level population numbers computed in local thermodynamic equilibrium (LTE) by using potentially different radiation and gas temperatures, as well as a modified spherical dilution factor. This modified spherical dilution factor is computed from the mean photospheric radius and optical depth, as is explained in, for example, Springmann (1997). Since the aim of the code is to work with non-monotonic velocity fields (which are fully accounted for, see Hennicker et al. 2020), the equation of radiative transfer is solved in the observer’s frame. This is done because a frequency-dependent comoving-frame formulation is very complicated to implement when dealing with highly non-monotonic velocity fields. Local radiation and gas temperatures were directly taken from the RHD simulation described in the previous section, and the dilution factor was computed from the radially averaged optical depth scale of the model. For continuum opacities and source functions, we assumed a constant Thomson scattering opacity and a Planck function using the local gas temperature. In Sect. 6, however, we also provide first test results from computation of scattering-dominated source functions. Boundary conditions are as follows (Hennicker et al. 2018, 2020): if rays originated from the stellar core the radiation at the lower boundary was set to $I_{\nu }^{+}=B_v(T_{rad},v)$, with the other quantities obtained from trilinear interpolation from the 3D model atmosphere. Inside the core, we used $I_{\nu }^{+}=S_L=S_C$ (with S_L and S_C the line and continuum source functions, respectively). At the outer boundary, the intensities that come from outside this boundary were set to zero and the intensities coming from inside were calculated with our RT scheme. For rays that do not intersect the core, these boundary conditions were applied at both sides of the hemisphere. Finally, in this work we only treated single lines and we further assumed pure Doppler profiles with the widths set by the local thermal speed of the considered ion.

3.2 Emergent intensities and fluxes

The specific intensity found from Equation (1) could be used to determine surface intensities of the 3D models. This emergent intensity was then integrated over the projected surface to obtain emergent fluxes. For planar models this is straightforward, but some care must be taken when accounting for sphericity effects. We used a cylindrical co-ordinate system (for example, Sundqvist et al. 2012; Hennicker et al. 2022):

$$F_{\nu }=\frac{1}{d^2}{\int }_0^{2\pi }{\int }_{-R_{max}}^{R_{max}}{I}_{\nu ,n}(p,\zeta ,z=R_{max})p\mathrm{d}p\mathrm{d}\zeta ,$$(2)

with d being the distance from the observer, R_max the maximum size of the emitting object, (p, ζ) the cylindrical co-ordinates describing the projected disc of the emitting object perpendicular to the observer’s direction, n, and I_ν,n(p, ζ, z = R_max) the emergent specific intensity into that direction evaluated at a distance z = R_max. The z-axis is aligned with the direction to the observer. The geometry is illustrated in Fig. 1 of Hennicker et al. (2022). Radiation is presumed to be free-streaming outside of the simulation domain, and as such the intensities at the outer radius do not change until the projected surface at z = R_max. If the rays go through the core, our limits are not −R_max to R_max, but the core radius up to R_max.

3.3 Building 3D spherical models

Before computing emergent flux line profiles, we reconstructed a global 3D spherical model from our local box-in-star RHD simulations. Since the width of the RHD simulations is much smaller than the radius of the star (as can be seen in Figures 2 and 3), we need many (in this specific case 100) simulation boxes (from here on called snapshots) to cover the full stellar and wind surface. Figure 4a gives a figurative illustration of how typical snapshots look (shape-wise). Once the snapshots were chosen from which the 3D spherical model would be constructed, the selected snapshots were stacked together (see Fig. 4b for an illustration of this process) and transformed into a sphere using cube-sphere transformations. The cubed-sphere co-ordinates split the sphere into six separate sectors (see Fig. 4c), where each of the sectors has its own local co-ordinate system (r, ξ, ζ). Here, r is the radial co-ordinate and ξ and ζ are two angular co-ordinates. The transformation from our pseudo-planar co-ordinates from the RHD simulations to the cubed-sphere co-ordinates, happens by interpreting the radial co-ordinate from the simulations as the radial co-ordinate on the cubed sphere, and interpreting the two lateral co-ordinates as the two angular ones. A final transformation turned this into a spherical star as illustrated in Fig. 4d. See for example Ronchi et al. (1996); Dimitrijevic et al. (2016) for more details on the cube-sphere transformation process.

A justified concern that should be noted is that the cubesphere transformation (going from Figs. 4b to 4d) stretches and distorts the atmospheric structures. This distortion grows linearly with the radius as we interpret the lateral co-ordinate in the RHD simulations as an angular co-ordinate in the construction of the sphere. Most importantly, however, the radial distance from the stellar core for each point stays the same. Sphericity effects have been accounted for as much as possible already in the RHD simulations (see for example Moens et al. 2022b, Appendix A) and as such the radial structure is consistent.

In this work, we used a set of 65 snapshots, each spaced approximately 0.5 t_d,w in time. We built two kinds of 3D spherical models. First, we created ‘snapshot-spheres’. For these models, one snapshot was copied 100 times to construct our stacked snapshot array, and in total 600 times to create our full 3D spherical model (since the stacked snapshot array is copied six times to form our sphere). Second, we created ‘mixed-spheres’. For these models, we mixed our 65 snapshots to build our spherical model. We created 20 mixed-sphere models (in which we made sure that the snapshot placement was different from the other mixed-sphere models). Note that since we needed 100 snapshots to assemble our stacked snapshot array, we allowed that snapshots be used (at most) twice during this construction process. An important detail to keep in mind here, is that both types of models (mixed- and snapshot-spheres) used the exact same set of snapshots during their formation process; that is, the snapshot-spheres and mixed-spheres used the same time series. During the spherical model construction and the spectral synthesis process, we degraded the 3D information, and also transformed it onto other grids. The final cylindrical grid on which the formal integral was solved, has 397, 157, and 800 points in (p, ζ, z), respectively. We tested that doubling this resolution has a negligible impact on the final resulting flux profiles.

Fig. 3 Emergent optical continuum intensities of local patches on the 3D O star model, for six different snapshots. Each panel here corresponds to one μ=1 (local surface normal and line of sight are aligned) 3D patch simulation, and is showing the complete patch. Colour bars display local emergent radiation temperature defined through Iν ≡ Bν(Trad). The spatial extent of each X and Y axis is 0.2 R0.

Fig. 4 Figurative illustrations of (a) the snapshots (3D unified atmosphere and wind simulations of O-type stars), (b) the stacked snapshots (stacked snapshot array), (c) one of the six sectors that is constructed from our stacked snapshot array and (d) the sphere that is constructed. For this illustration only one snapshot was used; hence a repetitive pattern is visible. The colour represents an arbitrary field in the simulation.

4 The resolved O-star surface

Before investigating the O-star surface of our spherical constructed models we first examine the individual snapshots. These are not subject to the deformations and stretching discussed in Section 3.3.

4.1 Continuum variation of local patch on O-star surface

Figure 3 displays emergent intensities for six snapshots of our 3D simulation, computed for a typical continuum wavelength in the optical assuming constant electron scattering opacity. The snapshots are separated by ~2000 seconds (which is in between the t_d,a ~ 1000 s and t_d,w ~ 10 000 s) and taken well after the model has relaxed and adjusted itself from the initial conditions. The figure shows how the O-star surface is characterised by very large spatial and time fluctuations; labelled with the intensities expressed as a radiation temperature of I_ν ≡ B_ν(T_rad), typical surface variations are as high as of the order of ~10⁴ K. Since the total spatial dimension of these local patches is 0.2 R₀, the sharp lanes as well as the larger bubble-like features visible in the figure are small in comparison to the star itself. Overall, this reveals a dynamically very active surface in stark contrast with the homogeneous surface assumed by present-day standard 1D models used for spectroscopic studies of O stars.

Fig. 5 Surface brightness plots of the O-star simulations described above, at the wavelength of the O III 5594 line: (a) for a snapshot-sphere model, (b) for a mixed-sphere model.

4.2 Surface brightness maps from spherical models

Using the 3D spherical surface volume reconstruction method described in Sect. 3.3, Figs. 5a and 5b next show two examples of simulated surface brightness plots for the complete star. The first panel in each figure corresponds to the continuum intensity, the second panel to the continuum + line total emergent intensity. We note the visual difference between Figs. 5a and 5b; the former displays a repetitive pattern, as the same snapshot was copied over the entire sphere (the ‘snapshot-sphere’ models). The latter, on the other hand, filled in the sphere using 65 distinct snapshots and as such avoids the regular patterns. Because of this, we focus primarily on these ‘mixed-sphere’ models in the rest of this work. For example, Fig. 6 displays the emergent intensity at line-centre across the surface for the three key diagnostic lines O III 5594 Å, Hα, and C IV 1548 Å. These lines were chosen for Fig. 6 because they represent three different kinds of important diagnostics in O stars (photospheric absorption line, wind-influenced recombination line, UV resonance wind line). As we in practice cannot resolve the O star surface in these wavelengths, the images in Fig. 6 are primarily for illustration purposes. Moreover, since these lines were calculated using an aNLTE approach, particularly the C IV 1548 Å line has to be treated with care as this line most definitely has a very scattering-dominated source function (see Sect. 6). Nonetheless, what already can be clearly seen from these images is that this (mixed-sphere) star is very non-uniform, which is in contrast with the typical assumption that massive stars are pretty smooth and stable at the photosphere.

5 Line profile results and trends

5.1 Line profile variability

Spectroscopic line profile variability has been observed in O-type stars (for example Fullerton et al. 1996; Kholtygin et al. 2003; Aerts et al. 2018). In our simulations, variability occurs naturally. Hence, let us first find approximate upper and lower limits of variability predicted by our RHD simulations. This variation in our line profiles can be looked at quite naturally, since it comes directly from our models (something that cannot be done in stationary 1D models). We focus the analysis in this section on three photospheric optical lines of carbon, nitrogen, and oxygen that are often used in spectroscopic analysis of O stars.

A lower limit can be found from the mixed-sphere models (see Sect. 3.3), since these average out the variability of individual snapshots (different assemblies of the same set of snapshots will give similar spectroscopic features). The snapshot-sphere models do not apply such averaging, and thus provide approximate upper limits to the variability of emergent flux profiles (different assemblies of different snapshots will give different spectroscopic features). Since the snapshots used to create the snapshot-sphere models (and also mixed-sphere models) are separated in time, one could say that the line profile variability is caused by time variability of the different snapshots.

Fig. 8a displays results from 65 snapshot-sphere models along with the mean flux profile, illustrating a large (and very likely overestimated) level of variability in three characteristic optical absorption lines stemming from this method. Figure 7 compares the typical level of variability from the two spherical reconstruction methods. As expected, the variability (that can also be seen in Fig. 8b) in the mixed-sphere models is much lower than in the snapshot-sphere models. This can be quantified by calculating the standard deviation of the variability. This standard deviation was calculated for each frequency-point or wavelength-point of the calculated profile. For clarity, we only provide the standard deviation (in units of normalised flux: F_tot/F_cont) at line-centre for each line. For the snapshot-sphere models, they are approximately 0.017 for the O III 5594 Å line, 0.018 for the C III 5697 Å line, and approximately 0.009 for the N IV 4059 Å line. For the mixed-sphere models, these values are approximately 0.0019 for both the O III 5594 Å and the C III 5697 Å line, and 0.0009 for the N IV 4059 Å line. This indeed shows that the variability of the mixed-sphere models is much lower than the variability of the snapshot-sphere models. Importantly, however, the average emergent profiles agree very well between the two methods (as can be seen in Fig. 8c), demonstrating that for the purpose of evaluating gross 3D effects the specific averaging choice is less important.

5.2 Averaged line profiles

We next analysed averaged line profiles, using measurements of average full width at half maximum (FWHM) and the associated standard deviation σ, as well as the equivalent width (EW). We note that for the calculation of σ, we assume that we are working with an isotropic Gaussian: $\sigma =\mathrm{F}\mathrm{W}\mathrm{H}\mathrm{M}/2\sqrt{2\ln (2)}$. The line averaging procedure is a simple averaging, where simulated data at each velocity point are added together and divided by the total number. The results of this averaging for the snapshotspheres can be found in Fig. 8a, where the red line corresponds to the averaged line, and Fig. 8c, where the blue line corresponds to the averaged line. For the mixed-spheres, the averaged line can be seen in red in Fig. 8b and in red in Fig. 8c. Although there is a clearly visible difference in the variability of these models (as can also be seen in Fig. 7) the averages of the line profiles are almost identical. This shows that even though different model-construction-methodologies were used the average line profiles that emerge are very similar. The resulting FWHM, σ, and EW for the mixed-spheres can be found in Table 2 (see Appendix A for corresponding values for snapshot-spheres, including Tables A1 and A2 for further details).

The results in Table 2 show that for the photospheric lines O III 5594 Å, C III 5697 Å and N IV 4059 Å stemming from our 3D mixed-sphere models, FWHM of the averaged lines lie in between 120 km/s and 130 km/s, with a corresponding σ in the range 52–55 km/s. The small differences in FWHMs and σ’s between the lines could be due to them having slightly different formation regions (and hence also slightly different turbulent velocities in the formation regions). Note that these characteristic widths were obtained here directly from the 3D model without any inclusion of ad hoc micro- or macroturbulence broadening; and indeed, they are on the same order as typical observed photospheric O-star line profiles (see the introduction).

Fig. 6 Surface brightness illustrations showing the emergent specific intensity as seen by an observer in a certain direction for the O III 5594 Å line, the H alpha line, and the C IV 1548 Å line at the line centre.

Fig. 7 Coloured-in regions showing variability for the O III 5594 Å line. The blue corresponds to the variability of the snapshot-sphere models (65 models in total), while the red corresponds to the variability of the mixed-sphere models (20 models in total).

5.3 Comparison with 1D, static atmosphere

We next compare the results of our synthetic photospheric line production to corresponding results of 1D codes by re-computing the lines from 1D radial averages of our snapshotsphere and mixed-sphere models, and setting all velocities to zero. This approach thus mimics what would come out of an equivalent 1D, static model photosphere analysis. For the mixed-sphere models, results are shown for the same three lines analysed above in Fig. 8d. The FWHMs, σ’s, and EWs of the mixed-sphere models are given in Table 2. As previously, corresponding values for the snapshot-sphere models can be found in Appendix A. Inspecting Fig. 8d and Table 2 one can see that, as expected, the line profiles for the 1D static model are much narrower (FWHM of the order of 10–20 km/s instead of 120–130 km/s), and also deeper. This is a direct illustration of the effect that the 3D velocity fields have upon photospheric absorption lines, offering a natural explanation for the need of ‘extra’ ad hoc broadening mechanisms in 1D model atmosphere codes. Additionally, we note that also the EWs of the lines are different, with the ones computed from 3D models clearly higher (see discussions below). This difference in EWs shows that this broadening is not just a macroscopic broadening effect, but also affects the overall strengths of the lines (most often mimicked in 1D codes by adding ad hoc ‘microturbulence’). Using our 1D static model lines (by artificially broadening them to fit the 3D lines), we can find what micro- and macroturbulence velocities are needed to reproduce the 3D results.

For the macroturbulence fit, we use two options: an isotropic Gaussian profile method and the radial-tangential (rad-tan) model by Gray (1975) with equal radial and tangential contributions. Fig. 9 shows that for isotropic macroturbulence we find best-fit values v_vmac ~ 70 km/s. This is consistent with the discussion above regarding velocity dispersion, since the characteristic velocity in this macroturbulence model is directly related to σ by $v_{\text{mac }}=\sigma \sqrt{2}$. For the rad-tan model we as expected find higher characteristic velocities (~110–120 km/s); it is a well-known result (e.g. Gray 1975; Sundqvist et al. 2013) that an isotropic macroturbulence model returns lower characteristic velocities than the rad-tan method. To reproduce the EWs of the 3D lines we need to apply an isotropic microturbulence ~10–20 km/s, where the different best-fit numbers for the different lines could reflect different formation heights or be a saturation effect caused by our neglect of scattering in the calculation of the source function (see discussion below). We note further that it is currently unclear to what extent 3D effects on EWs can be fully captured by such microturbulence, or whether it rather is indicative of 3D chemical abundance effects (see discussion in Section 7). Another important note is that the lines in Fig. 8d are slightly shifted towards the red in comparison to the 1D averaged lines. The shift values were approximately −17 km/s for the O III 5594 Å line, −15 km/s for the C III 5697 Å line, and −20 km/s for the N IV 4059 Å line. This shift is due to the fact that the photospheric surface is a dynamic one, and as such, the emitting surface is shifting the line. This net-infall of material around the photosphere was noted also in Debnath et al. (2024) and can be seen in the middle panel of Fig. 2, which shows that more than half of the surface has infalling material.

Fig. 8 O III 5594 Å line, the C III 5697 Å line and the N IV 4059 Å line for the 3D O-star simulation, (a) using 65 snapshot-sphere models (the red lines show the result when all lines are averaged), (b) using mixed-sphere models (the red lines show the result when all lines are averaged), (c) showing averaged lines from both the snapshot-spheres (blue line) and mixed spheres (red line) and (d) showing averaged lines for the mixed-sphere models for both the ‘3D’ (unaltered) models and the 1D averaged, velocity-field-turned-off models.

Fig. 9 Micro- and macroturbulent velocities found for one mixed-sphere model, using the 1D static lines as input and the 3D lines as pseudo-observations.

6 First line profile results including scattering

This section discusses first test results when line profiles are computed using scattering-dominated source functions. Basic assumptions and methods of this are the following (see Hennicker et al. 2020). The continuum and line source functions S_C and S_L are, respectively:

$$S_C=(1-{ϵ}_C)J_{\nu }+{ϵ}_C{B}_{\nu },$$(3)

$$S_L=(1-{ϵ}_L)\overline{J}+{ϵ}_L{B}_{{\nu }_0},$$(4)

with ϵ_C and ϵ_L the thermalisation parameters (0 is full scattering, 1 is full thermal), J_ν the mean intensity, $\overline{J}$ the profile-weighted mean line intensity and B_ν0 the Planck function at line-centre ν₀. The determination of source functions assumes a two-level atom with a thermal+scattering continuum for the line formation. For these first tests, we simply took both ϵ values as fixed input parameters (work is underway to instead compute them from atomic line and continuum data). Source functions are derived from 3D short-characteristics (SC) solutions to the equations of radiative transfer, discretised on a non-uniform Cartesian grid. The solution is augmented by an ‘accelerated-lambda-iteration’ (ALI) scheme, using non-local operators to ensure convergence. The procedure accounts fully for Doppler shifts from arbitrary 3D non-relativistic velocity fields, as described in detail by Hennicker et al. (2020). The line opacities were still calculated using the approximate NLTE approach described above. The Cartesian grid used for these calculations consisted of 219 points along each axis (so 110 in each positive or negative direction of x, y, and z). Using this part of our code package further allows us to compare short-characteristics (SC) fluxes to corresponding fluxes predicted by the analytic FLD relation used in our RHD simulations. In this respect, however, it is important to keep in mind that the strong frequency-dependent nature of fluxes makes it computationally prohibitive to apply the SC technique directly in the RHD simulations, at least in cases where opacity from a multitude of spectral lines is important. As such, we here instead compared FLD and SC fluxes by applying a simplified gray opacity model, and present results from this in Appendix B (see also Section 5.2. in Moens et al. 2022a for more such comparisons).

Including scattering typically reduces source functions in the line forming layers. As a consequence, this could potentially desaturate photospheric lines that may be saturated in our previous results. This is important for future applications, for example comparing synthetic lines to observations, since it can affect the line depths (and hence EWs). We note, however, that in our tests the FWHMs and σ derived above were barely affected by the inclusion of scattering. The deepening effect on the line stemming from scattering can be seen in Fig. 10 for the O III 5594 Å line, which shows that the scattering line source function is significantly depressed as compared to the aNLTE source function.

The scattering line source function also allows us to model the strong UV resonance lines that typically appear as P-Cygni lines in O-star spectra. Figure 11 shows the carbon UV resonance line at 1548 Å. In this figure, it is clearly visible that we have a much better P-Cygni wind profile compared to our aNLTE method. Note, however, that the underlying 3D RHD O-star simulation did not go out far enough in radius to capture the complete C IV 1548 Å line formation region, which means that we may expect the blue edge of the line to reach somewhat higher velocities than in this figure. Moreover, since the investigated C IV line in reality is a doublet, this must be accounted for before a direct comparison to observations can be done. Nonetheless, the figure illustrates how our 3D models are able to naturally capture key elements of such P-Cygni lines, such as the softening of the blue absorption edge by the turbulent velocity field. In summary, Figs. 10 and 11 show first promising results from using scattering-dominated source functions, taking us one step closer to performing quantitative spectroscopy using 3D models.

Fig. 10 O III 5594 Å line for one mixed-sphere model, using two different calculation methods. First, the previously shown aNLTE approach (the black line). Second, the scattering-dominated source function approach (red and blue line). Note that the difference between the red and blue line is the used values of the thermalisation parameters ϵC and ϵL.

Fig. 11 C IV 1548 Å line for one mixed-sphere model, using two different calculation methods. First, the previously shown aNLTE approach (the black line). Second, the scattering-dominated source function approach (red line).

7 Discussion and conclusion

When inspecting our 3D RHD simulation of a typical early O star in the Galaxy, we found turbulent velocities of approximately 48 km/s, alternatively measured as velocity dispersions (${\sigma }_{vd}^2=<v_r^2>-<v_r{>}^2$) of approximately 57 km/s at the average optical photosphere. Computing averaged photospheric absorption line profiles (see Table 2), we found broadening with standard deviation of approximately 54 km/s. That is, we may directly relate the magnitudes of turbulent velocities measured in the RHD simulations to the typical broadening expected for corresponding absorption line profiles. Moreover, as our fits to corresponding 1D models applying microturbulence and macroturbulence show, the characteristic photospheric velocities in our 3D models are broadly on the same order as typically observed for O stars (for example Simón-Díaz et al. 2010; Simón-Díaz & Herrero 2014; Simón-Díaz et al. 2017; Serebriakova et al. 2023, 2024).

This large natural broadening of massive-star absorption lines has been found in other recent 3D modelling as well (Schultz et al. 2023, see for example their Figs. 9 and 10), using Monte-Carlo radiative transport to produce absorption lines from the 3D ATHENA++ RHD simulations by Jiang et al. (2015); Schultz et al. (2022). Since their chosen massive stars to model did not have the same fundamental parameters as the O star presented here, the broadening results can however not be directly quantitatively compared. Also, in addition to different radiative transfer techniques, the underlying 3D RHD simulation methods differ between the two groups, as their ATHENA++ framework on the one hand uses a more sophisticated radiation closure (based on an Eddington tensor method), but on the other hand does not account for line driving effects (thus their simulations lack a wind outflow). Nonetheless, the overall qualitative agreement in results stemming from these two independent developments is very reassuring for future 3D massive-star model atmosphere work.

In addition to line broadening, the average atmospheric properties also change in multi-D simulations. Perhaps most strikingly for O stars, the strong turbulent pressure results in much larger photospheric scale-heights, and thus in significantly shallower density profiles (see also Debnath et al. 2024). As recently investigated by González-Torà et al. (2025) this affects in particular surface gravity determinations from observed spectral line profiles, and it is suggested that this might present a solution to the so-called ‘mass-discrepancy problem’ that is seen between evolutionary and spectroscopy mass determinations (Herrero et al. 1992). In this paper we have demonstrated how the photospheric turbulent velocity is directly related to the measured dispersion of absorption lines. While awaiting full 3D-based spectroscopic studies, one can thus use observed velocity dispersions in photospheric lines (at least in the absence of dominant rotational broadening) to calibrate the turbulent velocity needed to correct the hydrostatic equations used in all present-day 1D atmospheric models.

Finally, we note that we have potential 3D effects on abundance determinations, since the EWs found for the lines are different than in corresponding 1D models. This may lead to differences in derived abundances for the same objects, although the relative role of such 3D abundance corrections and the ‘microturbulence’ used in 1D models is not clear at the present (also since this strongly depends on the treatment of line scattering which will be further studied in future work). In comparison to studies on low-mass stars, abundance determinations done using 3D (time-dependent hydrodynamical) models of the Sun’s atmosphere have even changed the solar chemical composition yardstick (Asplund et al. 2009). For such low-mass stars in general, studies over the past decade have shown that results derived from 1D average structures are not adequate as substitutes for full 3D spectroscopy (Lind & Amarsi 2024).

In this paper, we have shown first promising results towards spectroscopy using 3D models also for hot, massive stars. In future work, we shall extend our formalism to include doublets, damping wings, as well as full scattering and NLTE effects, enabling a direct comparison to observations by means of quantitative spectroscopy.

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. The authors gratefully acknowledge support from the European Research Council (ERC) Horizon Europe under grant agreement number 101044048, from the Belgian Research Foundation Flanders (FWO) Odysseus program under grant number G0H9218N, from FWO grant G077822N, and from KU Leuven C1 grant BRAVE C16/23/009. The authors would also like to thank previous members of the KUL-EQUATION group for their earlier contributions. Finally, the authors would like to thank all members of the KUL-EQUATION group for fruitful discussion, comments, and suggestions. LD would also like to thank Conny Aerts for giving feedback on the draft, Levin Hennicker for the help and discussions in order to understand the general code package better and Stan Owocki for the discussions and insights. We thank the referee for their useful comments on our manuscript. We made significant use of the following packages to analyse our data: NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), matplotlib (Hunter 2007), Python amrvac_reader (Keppens et al. 2021), PyVista (Sullivan & Kaszynski 2019).

References

Appendix A Additional tables

Appendix B Short-Characteristics method versus Flux Limited Diffusion method

Here, we compare the Short-Characterics (SC) method that our radiative transfer code package uses for calculation of the source function in section 6 with the Flux Limited Diffusion (FLD) method that our RHD model atmospheres with winds simulations use. This is done in the following way: we take one snapshot, and apply our SC method on this snapshot to find the frequency-integrated Eddington flux using a pure absorption gray Thomson opacity model and the same 3D Cartesian technique as the corresponding SC calculations presented in Moens et al. (2022a). From this snapshot data, we then also derive the Eddington flux using the analytic FLD approximation as outlined in Moens et al. (2022b).

We also compute Eddington fluxes using a third method, namely SC solutions applied to the corresponding snapshotsphere model taking full sphericity effects into account. However, these Eddington fluxes cannot directly be compared to the RHD results, since this part of our radiative transfer code package (that deals with complete spherical models) is developed for post-processing purposes and so only contains frequency-dependent quantities. To remedy this, we set a T_rad,v using the Planck function from the SC Eddington fluxes of the corresponding snapshot-sphere model, calculate a frequency-integrated mean intensity $J=\frac{\sigma }{\pi }{T}_{rad}^4$ from this, and use this mean intensity to determine our Eddington flux using the analytic FLD approximation. This gives a rough indication of full sphericity effects upon the fluxes. The results of these comparisons can be seen in Fig. B1, where we have taken radial averages. Fig. B1 shows that in the inner layers of the atmosphere we find nice agreement between the different methods. However, for the outer optically thin wind layers there are clear differences. When inspecting the H_FLD curve we can see that the spherical correction we apply in the RHD simulations is, on average, ‘intermediate’ between the Cartesian SC result (H_SC) and the full spherical SC result (H_SC(spher), which clearly displays the expected 1/r² spherical dilution effect). Of course, in addition to sphericity one must not forget that the FLD closure relation itself will have a direct effect on the difference between SC and FLD fluxes; this will require further work in order to improve our closure relation in the RHD simulations. In Fig. B2 we show slices of the 3D RHD model at different optical depths for a given snapshot, focusing now on direct comparison to the 3D Cartesian SC fluxes. In deeper layers, there are only relatively small differences present, and in particular the main structures of the FLD and SC fluxes are very similar. In the outer optically thin layers, on the other hand, the SC flux becomes larger than the FLD flux; this again arises due to a combination of the approximate FLD closure and sphericity effects. Moreover, the SC fluxes in these outer regions are significantly smoother than the FLD fluxes; this is because the SC method integrates over large distances whereas the FLD flux is based on a local gradient approach.

Fig. B.1 1D radial averages of the Eddington flux, using an SC method, an FLD method, and a translated-from-spherical-SC method. The optical depth on the x-axis is the average continuum optical depth.

Fig. B.2 The Eddington flux calculated using the SC method, FLD method, and the ratio of them, for different vertical planes in the 3D RHD model with average continuum optical depths according to the labels. The lateral X and Y axes are labelled in units of R0.

All Tables

All Figures

	Fig. 1 Volume rendering of relative density (displayed colours are log10ρ/⟨ρ⟩) for the 3D model atmosphere and wind. Angle brackets denote lateral averaging at each radial position.
In the text

	Fig. 2 Maps of logarithm of density, radial velocity, and radiation temperature, in three selected vertical planes with laterally averaged continuum optical depth according to the legends. The lateral X and Y axes are labelled in units of R0.
In the text

Fig. 3 Emergent optical continuum intensities of local patches on the 3D O star model, for six different snapshots. Each panel here corresponds to one μ=1 (local surface normal and line of sight are aligned) 3D patch simulation, and is showing the complete patch. Colour bars display local emergent radiation temperature defined through Iν ≡ Bν(Trad). The spatial extent of each X and Y axis is 0.2 R0.

In the text

Fig. 4 Figurative illustrations of (a) the snapshots (3D unified atmosphere and wind simulations of O-type stars), (b) the stacked snapshots (stacked snapshot array), (c) one of the six sectors that is constructed from our stacked snapshot array and (d) the sphere that is constructed. For this illustration only one snapshot was used; hence a repetitive pattern is visible. The colour represents an arbitrary field in the simulation.

In the text

	Fig. 5 Surface brightness plots of the O-star simulations described above, at the wavelength of the O III 5594 line: (a) for a snapshot-sphere model, (b) for a mixed-sphere model.
In the text

	Fig. 6 Surface brightness illustrations showing the emergent specific intensity as seen by an observer in a certain direction for the O III 5594 Å line, the H alpha line, and the C IV 1548 Å line at the line centre.
In the text

	Fig. 7 Coloured-in regions showing variability for the O III 5594 Å line. The blue corresponds to the variability of the snapshot-sphere models (65 models in total), while the red corresponds to the variability of the mixed-sphere models (20 models in total).
In the text

Fig. 8 O III 5594 Å line, the C III 5697 Å line and the N IV 4059 Å line for the 3D O-star simulation, (a) using 65 snapshot-sphere models (the red lines show the result when all lines are averaged), (b) using mixed-sphere models (the red lines show the result when all lines are averaged), (c) showing averaged lines from both the snapshot-spheres (blue line) and mixed spheres (red line) and (d) showing averaged lines for the mixed-sphere models for both the ‘3D’ (unaltered) models and the 1D averaged, velocity-field-turned-off models.

In the text

	Fig. 9 Micro- and macroturbulent velocities found for one mixed-sphere model, using the 1D static lines as input and the 3D lines as pseudo-observations.
In the text

	Fig. 10 O III 5594 Å line for one mixed-sphere model, using two different calculation methods. First, the previously shown aNLTE approach (the black line). Second, the scattering-dominated source function approach (red and blue line). Note that the difference between the red and blue line is the used values of the thermalisation parameters ϵC and ϵL.
In the text

	Fig. 11 C IV 1548 Å line for one mixed-sphere model, using two different calculation methods. First, the previously shown aNLTE approach (the black line). Second, the scattering-dominated source function approach (red line).
In the text

	Fig. B.1 1D radial averages of the Eddington flux, using an SC method, an FLD method, and a translated-from-spherical-SC method. The optical depth on the x-axis is the average continuum optical depth.
In the text

	Fig. B.2 The Eddington flux calculated using the SC method, FLD method, and the ratio of them, for different vertical planes in the 3D RHD model with average continuum optical depths according to the labels. The lateral X and Y axes are labelled in units of R0.
In the text