1. Introduction
Very massive stars (VMS) are born with masses above 
$100 M_\odot$ and are driving engines in the evolution of their host galaxies. Their strong winds and extreme luminosities mean that – despite their rarity – they are able to contribute significantly to their environment through the injection of kinetic energy, momentum, and ionizing flux. Understanding their evolution is therefore key to studying their impact throughout cosmic time, especially in populations in which individual stars cannot be resolved (e.g., Crowther et al. Reference Crowther, Caballero-Nieves, Bostroem, Maz Apellániz, Schneider, Walborn, Angus, Brott, Bonanos, de Koter, de Mink, Evans, Gräfener, Herrero, Howarth, Langer, Lennon, Puls, Sana and Vink2016; Higgins et al. Reference Higgins, Vink, Hirschi, Laird and Sabhahit2023; Upadhyaya et al. Reference Upadhyaya, Marques-Chaves, Schaerer, Martins, Pérez-Fournon, Palacios and Stanway2024).
Among massive stars, VMS are unique in foregoing a cool supergiant phase of evolution, instead transitioning directly from O to WNh-type Wolf-Rayet stars while still undergoing core hydrogen burning due to their proximity to the Eddington limit and their resulting high mass-loss rates (e.g., Ekström et al. Reference Ekström, Georgy, Eggenberger, Meynet, Mowlavi, Wyttenbach, Granada, Decressin, Hirschi, Frischknecht, Charbonnel and Maeder2012; Köhler et al. Reference Köhler, Langer, de Koter, de Mink, Crowther, Evans, Gräfener, Sana, Sanyal, Schneider and Vink2015; Gormaz-Matamala et al. Reference Gormaz-Matamala, Romagnolo and Belczynski2025).
The strong winds of VMS make them challenging to model using stellar evolution codes, since these tend to assume a plane-parallel gray atmosphere. This deviation from the real physics of VMS introduces inaccuracies in radius and effective temperature right from the beginning of their evolution, which ultimately propagate back into the entire evolution itself, for instance through mass-loss prescriptions. For this reason, we aim to combine current evolution tracks with detailed atmosphere models of VMS in order to gain deeper insights into the underlying physics and resulting uncertainties that must be considered in this regime. The work presented here is a first step in this effort and constitutes a summary of the results of Josiek et al. (Reference Josiek, Sander, Bernini-Peron, Ekström, González-Torà, Lefever, Moens, Ramachandran and Schösser2025). We refer to this paper for a more in-depth discussion of the presented phenomena.
2. Evolution models of VMS
The evolution of VMS is highly dependent on the mass-loss rate, which sets how fast the helium-enriched layers are exposed on the surface. Recently, Sabhahit et al. (Reference Sabhahit, Vink, Higgins and Sander2022) proposed a new mass-loss recipe taylored to VMS, which consists of using the standard O-star mass-loss rate by Vink et al. (Reference Vink, de Koter and Lamers2001) with an additional boost (“kink”) as the star approaches the Eddington limit. We have implemented this mass-loss recipe in the stellar evolution code GENEC (Ekström et al. Reference Ekström, Georgy, Eggenberger, Meynet, Mowlavi, Wyttenbach, Granada, Decressin, Hirschi, Frischknecht, Charbonnel and Maeder2012; Eggenberger et al. Reference Eggenberger, Meynet, Maeder, Hirschi, Charbonnel, Talon and Ekström2008) and computed a set of VMS tracks along the main sequence, which we show in Fig. 1.
Hertzsprung-Russell diagram of VMS on their main sequence computed with GENEC. The open circles over the 
$150\,M_\odot$ track represent the snapshots used for the atmosphere models discussed in Sec. 3.
A striking feature in the models (especially at the lower-mass end of the considered range) is the inflation that they undergo during their main-sequence evolution. Closer inspection of the structure models reveals that this is purely a surface effect (as opposed to core-driven expansion) caused by the high opacity induced by iron ionization switched known as the “hot iron bump” (Sabhahit and Vink Reference Sabhahit and Vink2025). This high opacity translates to a high outward radiative force, which ultimately gets injected into the gas above this layer causing its expansion. Since this significantly affects the path of the star in the HRD, we conclude that our predictions of VMS evolution strongly depend on the treatment of the iron-bump region in the evolution codes. These predictions are then fed into calculations of stellar population and spectral synthesis, which ultimately results in the question of whether an LTE plane-parallel atmosphere is still an appropriate assumption for the iron-bump layers. Furthermore, we should consider the scenario where inflation is bypassed due to effects not included in the evolution code (for example the driving of turbulent motion).
3. Connecting model atmospheres
In order to address the potential shortcomings of evolution codes regarding VMS physics, we computed detailed stellar atmosphere models using PoWR (Gräfener et al. Reference Gräfener, Koesterke and Hamann2002; Sander et al. Reference Sander, Hamann, Todt, Hainich and Shenar2017), using our 
$150\,M_\odot$ evolution track as an input. When connecting atmosphere models to an underlying structure model, one has to deal with the fundamental decision of defining the inner boundary of the atmosphere model. There is a-priori no physically motivated reason to choose where to place this cut, and the effect of this arbitrary choice must be considered as part of the modeling uncertainties. In order to test this, we first divided the main-sequence evolution into 16 snapshots spaced roughly evenly in current stellar mass, and computed two atmosphere models at each snapshot connected at a different depth. Shallow connection: The standard way of defining an atmosphere model is by choosing some Rosseland mean optical depth 
$\tau$ that contains the entire line formation region of the atmosphere, motivated observationally by the use of atmosphere models for spectral fitting. Here we choose 
$\tau=20$. Deep connection: Instead of an optical depth requirement, we can alternatively choose the radius of the lower boundary of the atmosphere and assign it a corresponding density. For the models to cover the hot iron bump region, we chose a connection density of 
$\rho = 10^{-7}\,\mathrm{g}\,\mathrm{cm}^{-3}$ and define the radius to match the underlying evolution track. This is a new approach which allows us to ensure a continuous connection in the density profile between the interior structure and the atmosphere. If the atmosphere is in hydrostatic equilibrium
It is important to note that in the standard version of the PoWR code which we use here, the atmosphere density profile is not fully computed in a hydrodynamically consistent way. Instead, only the inner region of the atmosphere is integrated hydrostatically, and connected to an exterior wind regime with a prescribed velocity law. The velocity field sets the density profile in the wind via the continuity equation. Our sets of atmosphere models therefore do not show the true effect of the radiation field on the atmosphere structure itself, but represent two limit cases of a shallow versus a deep wind onset.
The evolution of the density profiles is shown in Fig. 2. One can clearly recognize the formation of a density plateau in the evolution models, which is avoided by the deep models. Consequently, the deep models represent a more compact and hot scenario for the main-sequence evolution of VMS.
Comparison of the density profile in shallow and deep atmospheres across evolution, resulting from a predefined velocity field.
4. Effects on spectra
We show in Fig. 3 the evolution of the emergent spectra for our 
$150\,M_\odot$ star in the two scenarios considered with our atmosphere modeling. The results show that in the absence of subsurface inflation (deep models), the spectra do not fully transition to pure WNh type morphology, as the H
$\beta$ line never reaches full emission. We conclude that subsurface inflation directly affects the abundance of expected WNh-type VMS, yielding a direct observational constraint on the nature of VMS atmospheres, even in layers so deep that they do not directly contribute to spectral line formation. The indirect effect of the iron bump on the atmosphere’s structure could well have consequences on the spectral signature of VMS stars.
Spectral evolution of a 
$150\,M_\odot$ star using shallow and deep atmosphere models. The bottom spectra represent the zero-age main sequence, and the evolution occurs upwards. The shallow model sequence is interrupted due to convergence issues for the maximally inflated models at the end of the main sequence.
Furthermore, the results show a striking difference in spectral predictions from the same atmosphere model using the same underlying physical stellar parameters, varying only the numerical boundary condition. This introduces an additional uncertainty that affects spectroscopic diagnostics of VMS as long as the structure of their atmosphere is not well-understood. Some of our test calculations resulted in a discrepancy of up to 0.4 dex in mass-loss rate determinations when using our deep vs. shallow models to fit observed spectra.
5. Conclusions and outlook
The atmosphere structure of VMS and their wind-launching mechanism remains an poorly understood problem of stellar astrophysics. Its treatment in stellar evolution codes is currently severely lacking and requires insights from other approaches in order to obtain a more accurate picture of VMS. In this work, we performed a simple modeling experiment to understand the effect of a deep vs. shallow wind in VMS. One set of models follows the evolution track close to the surface and follows the predicted inflation throughout the main sequence evolution. The other set of models is a new approach, changing the boundary conditions of the atmosphere models to suppress the inflated region and test the alternative possibility of compact VMS.
This experiment revealed striking differences in the spectra predicted by both approaches. These must be resolved using techniques that can directly assess the impact of the hot iron bump on the atmosphere structure. A promising avenue for this is the use of hydrodynamical modeling, including multi-D models with turbulence (Moens et al. Reference Moens, Poniatowski, Hennicker, Sundqvist, El Mellah and Kee2022). Previous such models have already suggested that the opacity of the hot iron bump can drive strong macroscopic turbulent motion and remove enough energy from the radiation field to significantly suppress inflation Debnath et al. (Reference Debnath, Sundqvist, Moens, Van der Sijpt, Verhamme and Poniatowski2024). Implementing this kind of opacity-depended turbulence directly into a 1D code such as PoWR remains a significant challenge, but insights from the multi-D simulations represent a valuable first step in this effort (see González-Torà et al. Reference González-Torà, Sander, Moens, Sundqvist, Debnath, Delbroek, Josiek, Lefever, Van der Sijpt, Verhamme and Bernini-Peron2025a; González-Torà et al. Reference González-Torà, Sander, Sundqvist, Debnath, Delbroek, Josiek, Lefever, Moens, Van der Sijpt and Verhamme2025b).
150M⊙
150M⊙
Open access
