© The Authors 2025

1 Introduction

Due to their extreme properties and role in galactic evolution, massive stars are of fundamental importance to our understanding of the Universe. They are the primary source of chemical enrichment of the interstellar medium, which is released during stellar outflows and supernova explosions. The individual evolution of a massive star is mainly controlled by its stellar wind combined with stellar rotation. Winds from massive stars are called line-driven winds; they are driven by intense radiation pressure where the radiation field transfers momentum to the stellar plasma by scattering in metal lines (Lucy & Solomon 1970). Currently, the hydrodynamic theory that describes these winds (m-CAK theory) is based on the pioneering work of Castor (1974), Castor et al. (1975), and Abbott (1982) and the improvements made by Friend & Abbott (1986) and Pauldrach et al. (1986); see a review from Curé & Araya (2023). From the standard m-CAK theory, the line-force parameters (α, k, and δ) provide scaling laws for the mass-loss rate $(\dot{M})$ and terminal velocity (v_∞) of the wind (Kudritzki et al. 1989). The parameters α and k describe the statistical dependence of the number of lines on frequency position and line strength (Puls et al. 2000), while the δ parameter is related to changes in the ionisation throughout the wind. Curé et al. (2011) found, in the frame of this m-CAK theory, a new spherical hydrodynamic solution (called the δ-slow solution), which is different from the standard m-CAK solution (called the fast solution). This solution is obtained when high values of the δ parameter (≳0.3) are considered. Although this value might seem high, Puls et al. (2000) showed that for winds with neutral hydrogen as a trace element, the δ value is ∼0.3. Moreover, for extremely low-metallicity winds, Kudritzki (2002) showed that δ might reach much higher values than 0.3.

The main characteristic of the δ-slow solution is its lower terminal velocities than the fast solution. This δ-slow solution might explain the observed terminal velocities from the winds of A-type supergiants (Curé et al. 2011) and B-hypergiants (Venero et al. 2024).

From an observational point of view, the latest generation of large telescopes has opened a wide range of possibilities for studying massive stars. Continuous superior instrumentation development and improving the quality of the spectroscopic data are essential to understand the nature of these objects better and their chemical and dynamic evolution. Moreover, supergiant stars present themselves as very promising distance indicators through the application of the wind-momentum luminosity relationship (WLR, Kudritzki et al. 1995) and the flux-weighted gravity-luminosity relationship (FGLR, Kudritzki et al. 2003). Hence, to fully understand the nature of massive stars, it is necessary to perform accurate analyses of quantitative spectroscopy on large samples of these objects to obtain their physical properties (Simón-Díaz 2020).

The only way to derive information from massive stars and their winds is by decoding their emitted radiation. Quantitative spectroscopy can extract numerous properties of massive stars and their winds through detailed non-local thermodynamic equilibrium (NLTE) modelling of their atmospheres and winds, generating a synthetic stellar spectrum for comparison with observational line profiles. Therefore, it is essential to use a state-of-the-art NLTE code capable of predicting a synthetic spectrum for a star with specific stellar and wind parameters and inhomogeneities. Whilst current multi-object spectrographs are undoubtedly capable of producing extensive collections of spectra (see, e.g. Evans et al. 2005; Simón-Díaz et al. 2011a), accurately modelling their atmospheres and winds is an inherently complex task due to the presence of NLTE processes. This complexity demands the use of advanced stellar modelling codes such as FASTWIND (Santolaya-Rey et al. 1997; Puls et al. 2005; Rivero González et al. 2012), CMFGEN (Hillier 1987; Hillier & Miller 1998; Hillier & Lanz 2001; Hillier 2012), and POWR (Hamann & Schmutz 1987; Todt et al. 2015).

Regarding analysis techniques, several alternatives have been proposed by different authors to minimise the subjective component present in the widely used ’by-eye’ techniques by introducing objective, automatic, and fast methods. These methods include using genetic algorithms employed by various studies such as Mokiem et al. (2005), Tramper et al. (2014), Hawcroft et al. (2021), and Brands et al. (2022). Another approach is a grid-based method, used by Lefever et al. (2010) and Simón-Díaz et al. (2011b), or the principal components analysis (PCA) algorithm designed by Urbaneja et al. (2008). Given the multi-dimensionality of the parameter space considered during the modelling process, the solution is not unique¹, and the existence of similarly acceptable solutions with different combinations of stellar parameters and/or abundances can be obtained. Indeed, correctly identifying the range of acceptable parameters and abundances is as essential as determining the best fitting (Simón-Díaz 2020).

In the pioneering work of Pauldrach et al. (1986), the authors fitted a set of hydrodynamical models with the law now known as the β-law, with three parameters: stellar radius, terminal velocity, and the β parameter. Nowadays, the β parameter is empirically used as input in these codes, with typical values of β = 0.7 to 1.5 (or even larger) found in O- and B-type stars (Puls et al. 1996; Kudritzki & Puls 2000). However, in the case of B- and A-type supergiant stars, the β value shows a tendency towards higher values (Lefever et al. 2007; Markova & Puls 2008; Haucke et al. 2018; Rivet et al. 2020; de Almeida et al. 2022, among others), and in some cases, values as high as 3–4 (Stahl et al. 1991; Crowther et al. 2006). These higher values, however, lack a physical and hydrodynamic justification, and their use is justified posteriorly by the quality of the line profile fits achieved (Kudritzki & Puls 2000).

In this work, we present the following elements:

• i)

  A new grid called the grId of Stellar atmOSphere and hydrodynamiC modELs for massivE Stars (ISOSCELES). The novelty of this new grid is that it uses the hydrodynamic solutions instead of the β approximation for the velocity field. The hydrodynamic calculations were obtained by the stationary code HYDWIND (Curé 2004) within the framework of the m-CAK theory. We used these hydrodynamical solutions as input in the FASTWIND code. In addition, these hydrodynamic solutions consider both the fast and δ-slow solutions.

• ii)

  A new automatic grid-based procedure was implemented over ISOSCELES, specifically designed to study OBA-type stars in the optical and infrared range. With this tool, we expect to derive stellar (T_eff, log g), wind (v_∞, M), and line-force (α, k, δ) parameters.

Our work is organised as follows: Section 2 describes our grid of models (hydrodynamic and stellar atmosphere), and Section 3 describes the methods for the spectral line fitting. The ongoing results are shown in Section 4, and, finally, the discussion and conclusions of our work are given in Section 5.

2 ISOSCELES grid

ISOSCELES is the first grid involving the m-CAK hydrodynamics (instead of the generally used β-law²) and NLTE radiative transport to generate synthetic line profiles for massive stars. ISOSCELES has been set up to cover almost the complete parameter space of OBA-type stars with solar metallicity, neglecting the effects of the stellar rotation³.

The winds of massive stars are often modelled using a clumping factor, a parameter that accounts for small-scale density inhomogeneities (or ‘clumps’) within the wind, rather than assuming a smooth, homogeneous outflow. In the current implementation of ISOSCELES, we adopted a clumping factor of unity (i.e. unclumped winds) across the entire grid. This choice allowed for a significant reduction in computational complexity when constructing millions of models; however, it also introduces notable limitations in the determination of wind parameters. In particular, density-squared diagnostics such as Hα may lead to overestimated mass-loss rates, since a clumped wind with the same emission profile would require a lower average mass-loss rate. Moreover, the exclusion of clumping neglects porosity and optical depth effects that can affect both optical and UV line profiles. As a result, the current grid may not fully capture the impact of wind inhomogeneities in stars with stronger outflows. Future versions of ISOSCELES will incorporate a parame-terised treatment of optically thin clumping to better constrain wind parameters and reduce degeneracies.

To produce this grid of synthetic line profiles with the code FASTWIND (version 10.1.7), we first computed a grid of hydrodynamic wind solutions with the stationary code HYDWIND for both m-CAK spherical solutions (fast and δ-slow). FASTWIND modified the lower part of the velocity profile from HYDWIND, replacing the corresponding quasi-hydrostatic stratification to account for a better description of the velocity gradient (instead of using the Sobolev approximation). The surface gravities range from log g = 4.3 down to approximately 90% of the Eddington limit, with increments of 0.15 dex. We considered 58 effective temperature values, from 9000 to 45 000 K, using steps of 500 K below 30 000 K and steps of 1000 K above it. Note that roughly below 13 000 K, the density structure predicted by FASTWIND becomes inaccurate because of missing Fe II lines in the optical, leading to underestimated radiative acceleration and consequently underestimated log g values. However, we have implemented ISOSCELES from 9000 K to avoid observed stars too close to the grid’s borders. For calculating the non-LTE stellar atmosphere models, we used micro-turbulent velocities of 8 km s⁻¹ (T_eff < 15 000 K), 10 km s⁻¹ (15 000 K ≤ T_eff < 20 000 K), and 15 km s⁻¹ (T_eff ≥ 20 000 K). For all synthetic line profiles (as calculated in the final formal integral of the FAST-WIND code), micro-turbulent velocities of 1, 5, 10, 15, 20, and 25 km s⁻¹ were employed.

In the current version of the grid, silicon (Si) is the only chemical element whose abundance is varied. We adopted five values of Si abundance, log ϵ_Si = 7.21, 7.36, 7.51 (solar), 7.66, and 7.81, based on the solar reference from Asplund et al. (2009). The choice to vary only Si at this stage is motivated by its strong diagnostic lines in the optical and near-infrared for B-type stars, such as Si III 4552 Å and Si IV 4212 Å, which are key for determining T_eff in supergiants and early main-sequence stars. The abundances of all other species are set to the solar values, and the He abundance is fixed at He/H = 0.10. Future updates of ISOSCELES will be made available to us, including variations in additional elements (e.g. C, N, O) to enhance its diagnostic power and applicability to different metallicities and evolutionary phases. Thus, the grid points were selected to cover the region of the T_eff–log g plane as shown in Fig. 1, where massive stars are located from the main sequence up to the blue super-giant phase. In this figure, each dot is further described by the other five parameters (see details in Section 2.2).

Venero et al. (2024) showed that the winds of some B super-giants might be described by the δ-slow solution. Therefore, we increased the number of models below 30 000 K to allow a more detailed analysis of the winds of BA supergiants.

Fig. 1 Location of ISOSCELES models in Teff–log g plane (red dots). Grey solid lines represent the evolutionary tracks from 2 M⊙ to 60 M⊙ without rotation (Ekström et al. 2012). Solid black lines correspond to both zero-age main sequence (ZAMS) and terminal-age main sequence (TAMS). Red (empty) circles indicate non-converged models in ISOSCELES.

2.1 Grid of hydrodynamic models

We computed a grid of hydrodynamic wind solutions with the code HYDWIND. These wind solutions are then input into the FASTWIND code. Each HYDWIND model is described by the stellar parameters (T_eff, log g, R ) and the line force parameters from the m-CAK theory (α, k∗, and δ). All these spherically symmetric models are calculated without rotation and satisfy the boundary condition τ = 2/3 at the stellar surface for the optical depth τ. For each given pair (T_eff, log g), the radius was calculated using the FGLR (Kudritzki et al. 2003, 2008), namely: $$M_{\text{bol}}^{\text{FGLR}}=3.41\times \left(\log {\text{gF}}_{}-1.5\right)-8.02,$$(1) where $M_{\text{bol}}^{\text{FGLR}}$ is the bolometric magnitude, and log g_F = log g − 4 log(T_eff/10⁴).

While FGLR was originally calibrated for supergiant stars, we applied this relationship more broadly to estimate stellar radii for all stars in the grid, including those in lower luminosity classes. This approximation is adopted to ensure a uniform and computationally consistent method across the entire grid. However, we acknowledge that the FGLR may not be strictly valid outside the supergiant regime, and the derived radii for dwarfs should be interpreted as internal scaling parameters within the grid rather than precise physical values. For quantitative spectroscopic analyses of individual stars, these radii can be later rescaled using observed magnitudes, distances, and extinction, as discussed in Section 2.2 (see below) and in Lefever et al. (2010). In future updates, we plan to incorporate evolutionary models with appropriate luminosity-dependent radius estimates for each class.

The range of the line-force parameters is shown in Table 1, where the values of δ determine fast (δ ≲ 0.24) or δ-slow solutions (δ ≳ 0.28). Due to numerical problems, not all combinations of these parameters converge to a physical stationary hydrodynamic solution (see, e.g. Venero et al. 2016). As δ-slow solutions converge less frequently than fast solutions, we used a denser grid of δ values from δ ≳ 0.26.

Topologically, fast and δ-slow solutions differ significantly in the singular (or critical) point location (Curé & Araya 2023). The fast solution has this point near the stellar surface, while the δ-slow solution has a singular point very far from this surface. This property is illustrated in Fig. 2, which presents a histogram of the locations of the singular points (r_crit) derived from the m-CAK equation of motion for the converged (fast and δ-slow) solutions.

2.2 ISOSCELES: Grid of stellar atmosphere models

A crucial aspect of analysing a massive star’s spectra involves creating a grid of synthetic spectral lines and comparing it to the observed spectra by reproducing their different features. Because of its high computational efficiency, we used the stellar atmosphere code FASTWIND. This code considers non-LTE effects in spherical symmetry, with an explicit treatment of the stellar wind effects such as line blanketing/blocking and clumping. We used a version that enables the hydrodynamic wind solution from HYD-WIND, ensuring a smooth transition between the pseudo-static photosphere and the lower wind layers. The transition point is located at v = 0.1v_sound, where v_sound is the sound speed.

The main advantage of FASTWIND concerning other similar codes is the possibility of generating realistic models in a short time interval, a crucial point for building a large dataset of synthetic line profiles. FASTWIND provides, in the optical and infrared wavelength ranges, H and He lines by default and lines from other elements, in dependence on purpose and available atomic data, in our case from Si. Using some of these lines, it is possible to determine the star’s photospheric properties and the characteristics of its wind. In the ISOSCELES grid, each FASTWIND model is described by the following eight main parameters, T_eff, log g, R_∗, α, k, δ, log ϵ_Si, and v_mic, corresponding to effective temperature, the logarithm of the surface gravity, stellar radius, line-force parameters (α, k and δ), the logarithm of the Si abundance, and the micro-turbulence velocity, respectively.

Regarding the radius, we must remember that a rescaling to the real value is required once concrete objects are analysed. The actual radius should be determined from the visual magnitude, the star’s distance, and reddening (Lefever et al. 2010).

In the end, 3 440 598 synthetic spectra models (considering the six micro-turbulent velocities) were obtained, each containing 56 line profiles; these lines are tabulated in Appendix A. Figure 1 shows models with solutions (red dots) in the (T_eff, log g) space. The grid was transformed to a binary format (FITS files) to reduce size. Currently, it has a size of 289 GB for the entire grid.

It is important to note that our method is not fully self-consistent since the radiation-wind interaction goes both ways. In principle, the wind structure and the spectrum must be obtained self-consistently through an iterative process. In addition, the m-CAK theory is considered a good approximation in the supersonic parts of the wind, but a questionable assumption in the sonic and sub-sonic portions due to the usage of the Sobolev approximation (e.g. Vink 2022 and references therein). For some self-consistent works avoiding this approximation, see, for example, Sander et al. (2017), Krtička & Kubát (2017), Gormaz-Matamala et al. (2019), and Björklund et al. (2021). However, despite these drawbacks, the results obtained from the ISOSCELES grid should have a better justification than the currently and widely calculated models with a β-law in quantitative spectroscopic analyses. Also, the ISOSCELES grid is helping us to clarify if the δ-slow solution is a good candidate to solve some discrepancies, for example between predicted and observed mass-loss rates or the high values of β obtained for some velocity laws.

Another assumption considered in our models is a smooth wind, i.e., a wind without clumps (regions of overdensity). These regions could become optically thick (especially in strong resonance lines), leading to porosity effects that may affect the line profiles and the mass-loss rate measurements. Similarly, optically thin clumps affect the mass-loss rate measurement from lines depending on ρ². Then, wind diagnostic lines that ISOSCELES cannot reproduce could be potential targets of interest for detailed investigations considering clumpy winds.

While the use of HYDWIND velocity fields instead of a prescribed β-law might mitigate certain degeneracies in stellar and wind parameter determination, the exclusion of clumping as a free parameter introduces others. In particular, line diagnostics sensitive to ρ² (e.g. Hα) can be reproduced with different combinations of mass-loss rate and clumping factor, making it difficult to uniquely constrain $\dot{M}$ in the absence of additional observables or assumptions.

Fig. 2 Histogram showing rcrit values from converged solutions. Blue bars represent fast solutions, and red bars correspond to δ-slow solutions. Both axes are on a logarithmic scale.

3 Quantitative analysis procedure

To identify the optimal synthetic model for each star under analysis, we developed a semi-automatic code that analyses a set of different observational line profiles for hydrogen, helium, and silicon. Specific lines of these elements can be selected to extract key stellar parameters and information about the wind parameters.

Each selected line contributes differently to the determination of stellar and wind parameters. In general, the Balmer lines (Hγ and Hδ) are primarily used to constrain the surface gravity. The Hα profile is highly sensitive to the mass-loss rate $(\dot{M})$ due to its strong dependence on the wind density.

Different helium ionisation states are good indicators to estimate the effective temperature values: He I 4471 and 6678 Å increase in strength towards cooler temperatures, while He II 4541 and 4686 Å are stronger in hotter stars and can also serve as wind indicators in the presence of emission components. Similarly, the degree of ionisation of the silicon lines such as Si II 4130, Si III 4552, and Si IV 4212 Å helps constrain T_eff. This analysis is sensitive to silicon abundance variations, accounted for in the ISOSCELES grid. The combined fitting of these lines allows for a robust determination of the stellar parameters across a wide range of spectral types. The code additionally requires the rotational parameter, v sin i, and the macroturbulence velocity, v_macro, as inputs, which were previously obtained with the IACOB-BROAD⁴ tool (Simón-Díaz & Herrero 2014).

Optionally, the user can provide estimated ranges for the effective temperature, surface gravity, and solution type (fast or δ-slow) to optimise the search procedure. If these optional parameters are not provided, the code uses effective temperatures from the Simbad astronomical database⁵ to constrain the search process or uses the entire grid if there is no information. On the other hand, log g is exhaustively searched throughout the whole ISOSCELES grid.

The code utilises multi-processing to enhance efficiency, leveraging all available processors to execute: (i) convolution routines on synthetic line profiles; (ii) interpolate the convolved line profiles within a user-defined spectral window at the wavelength of the observed spectra; and (iii) perform a χ² test, namely $${\chi }^2=\frac{1}{n_{\text{lines\hspace{0.17em}}}}{\displaystyle \sum }_j^{n_{\text{lines\hspace{0.17em}}}}\frac{1}{n_v}{\displaystyle \sum }_i^{n_v}\frac{{\left(y_{\text{i},\text{j}}-y_{\text{obs}}\right)}^2}{{\sigma }_{\text{L}}^2},$$(2) where n_lines represents the number of spectral lines to be fitted, y_obs is the normalised observed flux, y_i,j is the flux of the convolved synthetic model, σ_L is the inverse signal-to-noise ratio (S/N) of the observed spectrum, and n_ν corresponds to the number of wavelengths from the observed spectral line. This χ² test was adapted from Holgado et al. (2018).The code compares the observed spectral lines against ISOSCELES’s convolved and interpolated synthetic line profiles.

Finally, the code compiles the χ² results, sorts them in ascending order, and identifies the best-fitting model, defined by the lowest χ² value, averaged with equal weighting across all the used line profiles⁶. The best model is then exported to a PDF file that includes a graphical comparison of the synthetic spectral lines with the observational ones.

The uncertainties in the derived stellar and wind parameters are primarily determined by the resolution of the ISOSCELES, which is defined by the step size of each parameter within the grid. Here, we did not account for the implementation of determined systematic errors of Holgado et al. (2018); this study will be the scope of future work. The database is constructed using discrete intervals for key parameters such as effective temperature, surface gravity, and the line-force parameters (α, k, and δ). As a result, any derived parameter’s uncertainty corresponds to the interval’s size between adjacent grid points⁷. For example, the step size for effective temperature (T_eff) and surface gravity (log g) directly defines the range within which the best-fit values may vary. Although this approach provides a systematic and well-defined estimate of uncertainties, future improvements could involve interpolating between grid points or employing Bayesian methods to further refine the estimates and reduce parameter uncertainties.

To move beyond assigning uncertainties based solely on the grid resolution, we conducted an additional analysis, presented in Appendix B, to evaluate the sensitivity of the error estimates to the number of top-ranked models considered. Following a quantile-based method similar to that of Turis-Gallo et al. (2025), we tested different values of N and found that the resulting uncertainties, particularly for wind parameters, depend strongly on this choice. This highlights the limitations of arbitrary model selection thresholds and motivates the future implementation of a more rigorous statistical framework.

As an additional robustness check, we also performed a separate temperature fitting using only helium and, alternatively, only silicon lines, while keeping the wind parameters fixed. Both diagnostics returned effective temperature values consistent within the estimated uncertainties. The details of this test are described in Appendix B.

4 Quantitative spectroscopy

To test the capabilities of the ISOSCELES database in accurately deriving stellar and wind parameters across a range of physical conditions, we selected six stars, a subset of three dwarf stars and another with three supergiant stars; each subset has a star that lies on the hot side of the bi-stability jump (see, e.g. Vink 2022, hereafter BS-jump). The selection of stars on both sides of the BS-jump allows us to test whether the ISOSCELES framework can account for the different wind behaviours associated with this transition. Information about the used observations (dates, SNR, and instruments) is summarised in Appendix C.

Considering that the mass-loss rate, $\dot{M}$, strongly determines the shape of the Hα line profile, the selected stars show different shapes of this line, allowing us to test the performance of our ISOSCELES procedure for a wide range of $\dot{M}$. We analysed a set of different spectral lines, including Hα, and contrasted the results with previous studies that use CMFGEN/FASTWIND codes with the β-law velocity profiles.

4.1 Analysis of dwarf stars

In this section, we test the ability of the ISOSCELES grid to model the spectra of dwarf stars. Since the wind of these stars is too weak, we only focus on the stellar parameters, T_eff and log g, testing the photospheric structure.

As the analysis of line profiles is highly sensitive to the region where the transition between the quasi-hydrostatic layers and the stellar wind occurs, which is very close to the stellar surface, it is crucial to assess whether the inclusion of m-CAK wind solutions affects the outer layers of the photosphere. The selected dwarf stars representing different stellar types are HD 14633, HD 35299, and HD 35912.

4.1.1 HD 14633

This is an ON8.5 V star, and the spectrum is from the IACOB⁸ database (Simón-Díaz et al. 2011a). Figure D.1 (see Appendix D) shows nine spectral lines and their best fast and δ-slow fits, where the fast solution is the better one. Some fitted line profiles are commonly used as a criterion to determine stellar parameters, summarised in Table 2, where we compared our results with the ones from Aschenbrenner et al. (2023).

The main differences between the observation and the synthetic line profiles are found in the core of the lines formed in the photosphere and the wind. The fast wind solution fits most of the line profiles better. Note that while He II 4686 Å is well described, the core of Hα is too shallow for the fast solution.

4.1.2 HD 35299

This second tested dwarf star is a B1.5 V from the IACOB database. Figure D.2 shows a set of nine spectral lines that are best suited for this spectral type (see, e.g. Evans et al. 2015). We compare our results (see Table 3) with the ones from Nieva & Przybilla (2012). Both wind regimes predict almost similar line profiles, implying that the mass-loss rate is very low.

4.1.3 HD 35912

Our last dwarf star has a spectral type B2V; we also obtained an observed spectrum from the IACOB database. Figure D.3 and Table 4 summarise the results of the fitting procedure and a comparison with the work from Nieva & Simón-Díaz (2011). The predicted line profiles from both wind regimes exhibit a high degree of concordance.

As expected, the solutions for absorption lines are very similar for both δ-slow and fast regimes, indicating that the wind has minimal impact on these lines. Exceptions to this include Hα and He II 4686 Å, where more pronounced differences are observed in the cores due to the ‘fill-in’ effect caused by the tenuous wind.

As a final remark of this section, the photospheric structure of the ISOSCELES models computed using HYDWIND and FASTWIND delivers values of the stellar parameters close to published ones. This gives us confidence in fitting more complex atmospheres that include strong winds, as in the case of supergiant stars.

The successful fits obtained for the three dwarf stars are largely expected, given that their winds are extremely weak and have negligible influence on the line formation in the optical regime. Consequently, the hydrodynamic input from HYDWIND plays a minor role, and the synthetic spectra are dominated by the NLTE atmospheric structure and radiative transfer handled by FASTWIND. In this sense, our comparison with previous analyses based on DETAIL/SURFACE or TLUSTY models (e.g. Nieva & Simón-Díaz 2011; Nieva & Przybilla 2012) serves more as a consistency check between independent codes rather than a full test of the hydrodynamic component of ISOSCELES. Nonetheless, this validation step is important to ensure that the grid can reliably recover stellar parameters for stars with negligible winds before moving on to more complex cases, such as supergiants, where the wind solution significantly affects the atmospheric structure and emergent spectrum.

Fig. 3 Hα spectral-line fits for O-type star HD 14947. The black line represents the observed spectrum, while the best fit from the ISOSCELES database with a fast solution and δ-slow solution are shown in blue and red, respectively. The fast regime provides a better fit to this line profile.

4.2 Analysis of supergiant stars

In this section, we analyse three early-type supergiant stars, all with P-Cygni or emission Hα line profiles, which are very sensitive to the mass-loss rate. Due to this feature, we added the wind parameters in addition to the stellar ones.

4.2.1 HD 14947

The star HD 14947, classified as spectral type O4.5 If (Sota et al. 2011), was studied using spectroscopic data from the IACOB database. Figure 3 presents the Hα line profile, while Fig. 4 displays eight additional spectral lines. The results show that the fast solution yields a better overall fit to the observed profiles than the δ-slow solution. However, it is important to note that the poor fit in the wings of the He II 4686 Å line suggests that the outer velocity field is not accurately reproduced. This discrepancy highlights the need for further refinement in modelling the wind structure using UV lines. The velocity profiles used as input for generating the synthetic spectra are displayed in Fig. 5.

Table 5 shows stellar and wind parameters for both hydrodynamical solutions. For comparison purposes, in the last column, we add the results from Bouret et al. (2012), which used a β-law description for the velocity field. It is worth noting that ISOSCELES is a grid of models, all with a clumping factor of one ( f_cl = 1), while the results from Bouret et al. (2012) consider an optically thin clumping factor with a value of f_cl = 33 ( f_∞ = 0.03).

4.2.2 HD 206165

The star HD 206165, classified as spectral type B2 Ib (Negueruela et al. 2024), was analysed using spectroscopic observations from the IACOB database. Figure 6 shows the Hα line profile, where the δ-slow solution fits this line profile much better than the fast one. Eight additional lines are shown in Fig. 7. The synthetic lines from both wind solutions are indistinguishable, suggesting that Hα is the only sensitive line in our set of optical line profiles to distinguish the different wind regimes. Fig. 8 shows the corresponding velocity profiles. Table 6 summarises the stellar and wind parameters from our fits and the results from the work of Markova & Puls (2008), which used a β-law to describe the velocity field.

Fig. 4 Best spectral-line fits for HD 14947. The black line depicts the observed spectrum. The best fit with a fast solution is shown in blue, while the best fit with a δ-slow solution is in red. Clearly, the fast solution better describes the wind of this early-type star.

Fig. 5 Velocity profiles as function of inverse radial coordinate u for the best-fit solutions of HD 14947. The velocity profile for the fast solution is shown in blue, while the δ-slow solution is depicted in red. Additionally, the green line represents the β-law velocity field, based on the parameters from Bouret et al. (2012).

Fig. 6 Same as Fig. 3, but for HD 206165. Here, the δ-slow solution is the one that best fits the observed line profile

4.2.3 HD 53138

The star HD 53138 is a blue supergiant of spectral type B3 Ia (Negueruela et al. 2024). Figure 9 shows the Hα line profile, and Figure 10 shows the other eight available lines from our spectrum used in the fitting procedure. The δ-slow solution fits the observed profile better than the fast solution, mainly because the red wing of the Hα line has a better agreement with the observation. Fig. 11 shows the corresponding velocity profiles.

In Table 7, some stellar and wind parameters for both hydrodynamical solutions are tabulated. For comparison purposes, in the last column, we added the results from Haucke et al. (2018), which used a β-law description for the velocity field.

The standard fitting procedure for massive star winds typically relies on the empirical β-law for the velocity field and allows for different values of the clumping factor, introducing degeneracies between the derived mass-loss rate and terminal velocity. In this work, we adopted a different approach: we employed self-consistent hydrodynamic m-CAK solutions (both fast and δ-slow) for the wind-velocity structure, but we fixed the clumping factor to f_cl = 1 (i.e. assuming unclumped winds). This strategy significantly reduces the degeneracy between $\dot{M}$ and v_∞ by eliminating the additional free parameters associated with the β-law and variable clumping, allowing for a more physically motivated determination of wind parameters. Nevertheless, it is important to note that using a single value of the clumping factor is a simplification; real stellar winds are expected to be inhomogeneous, and future works will need to explore a range of f_cl values to fully capture the impact of clumping on wind diagnostics. While optical lines such as Hα remain sensitive to both $\dot{M}$ and v_∞, UV resonance lines provide crucial constraints on v_∞ and the wind’s ionisation structure, helping to break these degeneracies. The comparison with recent UV-based studies (e.g. Bernini-Peron et al. 2023, for HD 206165 and HD 53138) highlights the added value of UV diagnostics, which often yield higher terminal velocities than those inferred from optical-only, δ-slow fits. This motivates our planned extension of ISOSCELES to include UV synthetic spectra and integrated optical+UV fitting routines. By combining hydrodynamic wind solutions with a more comprehensive treatment of clumping in future work, we aim to minimise degeneracies further and improve the reliability of derived stellar wind parameters.

As initially motivated, the selection of the six targets was designed to sample stars on both sides of the so-called bi-stability jump (BS-jump), a regime where significant changes in wind properties, particularly terminal velocity and mass-loss rate, are observed near T_eff ∼ 21 000–25 000 K (see, e.g. Vink 2022). Interestingly, our results show that the δ-slow hydrodynamic solution systematically provides better fits to the cooler stars below the BS-jump, such as HD 206165 and HD 53138, while hotter stars such as HD 14947 are better described by the fast solution. This pattern suggests that ISOSCELES may offer a physically motivated mechanism, within the m-CAK theory, for interpreting the BS-jump as a bifurcation between distinct wind regimes. While this result alone does not constitute a predictive proof of the BS-jump, it supports the idea that the δ-slow solution could provide a consistent hydrodynamic explanation for this observational phenomenon. Further investigation is warranted using a larger sample and, ideally, the inclusion of UV diagnostics to evaluate whether the δ-slow/fast dichotomy aligns consistently with the BS-jump across the HR diagram.

Fig. 7 Same as Fig. 4, but for HD 206165. Both solutions show the same behaviour, demonstrating that these lines are not as sensitive to the mass-loss rate and the velocity profile as Hα.

Fig. 8 Same as Fig. 5, but for HD 206165. The β velocity profile is from the analysis of Markova & Puls (2008).

Fig. 9 Same as Fig. 3, but for HD 53138. In this case, the δ-slow solution is the one that best fits the observed line profile.

Fig. 10 Best spectral-line fits for HD 53138. The black line depicts the observed spectrum. The best fit with a fast solution is shown in blue, while the best fit with a δ-slow solution is in red.

Fig. 11 Same as Fig. 5, but for HD 53138. The β velocity profile is from the analysis of Haucke et al. (2018).

5 Discussion and conclusions

The ISOSCELES database, introduced in this study, represents an advancement in the quantitative spectroscopic analysis of massive stars. We propose that by incorporating hydrodynamic wind solutions derived from the m-CAK theory, instead of relying solely on the conventional β-law approximation, we provide a more physically robust framework for modelling stellar winds. Including fast and δ-slow solutions within ISOSCELES potentially permits a more precise representation of wind dynamics. This is particularly relevant for addressing previously observed discrepancies, specifically concerning terminal velocities and mass-loss rates in B-supergiants, where the δ-slow solutions might offer a better fit. However, multi-wavelength analyses and a larger sample of stars are still necessary to draw a definitive conclusion. We plan to address this in future work to maximise the potential of our approach.

In applying ISOSCELES to a limited, yet carefully chosen, sample of dwarf stars (HD 14633, HD 35299, and HD 35912) and supergiant stars (HD 14947, HD 206165, and HD 53138), we have observed its ability to reproduce their optical spectral profiles with reasonable fidelity. For the O-type star HD 14947, the fast solution delivered a superior fit to the Hα line profile, aligning with expectations for stars possessing high terminal velocities. This is consistent with previous studies and provides a degree of validation for the fast solution within the ISOSCELES framework. Conversely, for the two B-type supergiant analysed, the δ-slow solution demonstrated a better fit, particularly in the Hα profile. This underscores the potential of the δ-slow solution in modelling stars with lower terminal velocities on the order of a few hundred km/s and showcases the advantages of ISOSCELES. These findings, while promising, are preliminary, and we agree on the need for further refinement and expansion of the ISOSCELES database.

The derived α values from our best-fit models fall within the expected theoretical range for line-force parameters. However, for B-supergiants, the retrieved values are slightly higher than those predicted by theoretical studies, such as those presented by Puls et al. (2000). This deviation may reflect the influence of additional physical processes not fully captured in current hydrodynamic models. This result underscores the importance of further investigation into the parametrisation of radiative acceleration. It suggests that improved empirical constraints alongside refined theoretical modelling will be crucial for advancing our understanding of radiation-driven winds in this regime.

Given the computational expense of incorporating clumping, the assumption of smooth, unclumped winds in our models is an (currently) inevitable/forced simplification. However, we acknowledge that clumping in actual stellar winds can influence the observed spectra and derived parameters, particularly affecting Hα and UV lines, leading to overestimates of mass-loss rates (see, e.g. Hawcroft et al. 2021).

Similarly, while computationally expedient, the dependence on the Sobolev approximation may not fully encapsulate the intricate physics of the sonic and subsonic wind regions. In these regions, the velocity field is either strongly curved or its gradient is slight, and the Sobolev approximation can break down. However, as previously noted, this limitation does not affect our analysis, as we employed an alternative approach for calculating the wind structure that does not depend on the Sobolev approximation in these regions.

Furthermore, although non-solar metallicity environments are crucial for understanding stellar evolution across different galaxies, they have not yet been explored in this study. Applying our database to other stellar populations, such as those in the Magellanic Clouds or other nearby galaxies, will require the development of additional grids with appropriately adjusted metallicities. We plan to develop a more comprehensive grid of models, including an optimised χ² fitting procedure and, where justified, a finer sampling of δ values. This approach will enable a more robust and quantitative exploration of the parameter space than is possible through visual comparison alone.

Incorporating grids for supergiants with varying metallicities is essential to fully leverage the X-Shooting ULLYSES project (Vink et al. 2023). This will enhance our ability to constrain stellar and wind parameters with greater precision, ultimately improving our understanding of massive stars in diverse galactic environments.

It is worth emphasising that ISOSCELES is still a developing project with the potential to evolve in several important directions. One key avenue for improvement involves refining the physics incorporated into the grid generation. This includes addressing current limitations, such as the treatment of clumping and the reliance on the Sobolev approximation in the wind calculations, both of which can impact the accuracy of wind parameter estimates. Additionally, a more precise description of the velocity field could be achieved by incorporating a detailed treatment of line overlaps, as proposed by Poniatowski et al. (2022). Their iterative procedure self-consistently accounts for the complex interactions between radiation and matter, providing a more accurate representation of the wind dynamics and radiation field coupling. Implementing such enhancements would further improve the reliability of ISOSCELES in modelling massive star winds and interpreting observational data.

Secondly, we could explore alternative techniques for fitting observed spectra, such as employing machine learning algorithms (see, e.g. Olivares et al. 2024) or Bayesian inference methods to enhance parameter estimation and provide a more robust quantification of uncertainties. Another promising approach is to deconvolve the observed line profile to that of a non-rotating star (see, e.g. Escárate et al. 2023). This would eliminate the need to convolve all synthetic profiles during the comparison process, significantly reducing computational time and resource requirements.

Crucially, we plan to extend the ISOSCELES analysis to include the UV spectral range, which is rich in diagnostic information about the winds of massive stars. The UV spectrum contains numerous resonance lines that are highly sensitive to the wind density and velocity structure. By simultaneously fitting spectral features in both the optical and UV ranges, we can achieve much tighter constraints on key wind parameters, including the terminal velocity, mass-loss rate, and the parameters governing the velocity law. This multi-wavelength approach is essential for breaking degeneracies that arise when relying solely on optical lines, leading to a more comprehensive and robust characterisation of the stellar wind. This enhancement will be integrated into the ISOSCELES database, improving the accuracy of derived physical and wind parameters.

One of our future objectives is to incorporate CNO elements into the ISOSCELES framework to further enhance its precision in parameter derivation. UV data are critical as they provide firm constraints on the properties of stellar winds. The UV range encompasses many important spectral lines, including resonance lines, intercombination lines, and subordinate lines. Resonance lines, in particular, are invaluable for constraining low mass-loss rates (of the order of ∼10⁻⁹M_⊙/yr), where the Hα line becomes insensitive (Hillier 2020).

The synergy between ISOSCELES and UV spectra will offer a significant potential for advancing our understanding of mass loss in B supergiants. Another key advantage of UV lines, especially those involving resonance transitions and large optical depths, is their ability to provide precise measurements of terminal wind velocities. As a result, incorporating the UV regime will be critical for obtaining accurate and reliable information on terminal velocities and other wind properties, thereby enhancing the overall effectiveness of our spectroscopic models. However, to fully exploit the diagnostic power of the UV, it is essential to account for optically thick clumping.

In conclusion, we present the ISOSCELES database as an evolving and valuable resource designed to benefit the scientific community, especially with its planned extension to the UV spectral range. While we acknowledge its limitations, we are optimistic that ISOSCELES will contribute significantly to our understanding of massive stars and their role in galactic evolution with ongoing development. Its ability to handle a wide range of stellar parameters, along with its potential for expansion to different metallicities and wavelength ranges, positions it as a promising tool for analysing large datasets from current and future observational facilities, such as the ELT and JWST.

We believe ISOSCELES will be instrumental in interpreting the vast amounts of data expected from these next-generation instruments, and we aim to enhance its accuracy and efficiency over the coming years. To foster collaboration and support further advancements in the field, we welcome community engagement and are open to collaborations upon reasonable request. We also plan to make the ISOSCELES database available in the future.

Acknowledgements

The authors thank Joachim Puls for his valuable comments and insights on this manuscript. We also thank the anonymous referee for the valuable comments and suggestions, which helped to improve the clarity and quality of this manuscript. The authors are grateful for support from ANID / FONDO 2023 ALMA / 31230039. M.C., I.A., and C.A. thank the support from ANID FONDECYT project 1230131. M.C. and C.A. acknowledge partial support from Centro de Astrofísica de Valparaíso. N.M. thanks the support from ANID BECAS / DOCTORADO NACIONAL 21221364. R.O.J.V. acknowledges financial support from CONICET (PIP 11220200101337CO) and the Univer-sidad Nacional de La Plata (Programa de Incentivos 11/G192 and 11/G193). This work has been possible thanks to the use of AWS-U.Chile-NLHPC credits. Powered@NLHPC: this research was partially supported by the NLHPC’s supercomputing infrastructure (ECM-02). This project has received funding from the European Union’s Framework Programme for Research and Innovation, Horizon 2020 (2014-2020), under the Marie Skłodowska-Curie Grant Agreement No. 823734, and is also co-funded by the European Union (project 101183150 -OCEANS). We also thank Graeme Candlish and Omar Cuevas for allowing us to use part of their computer facilities for the initial calculation of our grid.

References

Appendix A Lines available in ISOSCELES

The lines available in ISOSCELES are listed in Table A.1 and Table A.2 for the optical and infrared ranges, respectively.

Appendix B Uncertainties estimation sensitivity

In this appendix, we present an analysis of the sensitivity of uncertainty estimates based on the number of top-ranked models included in the calculation. Models are sorted by increasing χ² value, and uncertainties are derived from the distribution of a selected number N of best-fitting models, following a quantile-based approach similar to that of Turis-Gallo et al. (2025). Specifically, we use HD 206165 and HD 14633 as test cases and compute the 0.25 quantile (q_0.25) for N = 3,000, 5,000, 10,000, 15,000 and 20,000. The results are presented in Table B.1, additionally, the table includes the best-fit parameter values.

As expected, the estimated uncertainties increase with N, especially for wind-related parameters such as $\dot{M}$ and v_∞, which are more degenerate. However, the growth rate diminishes for large N, suggesting a saturation effect. For instance, the uncertainty in $\dot{M}$ for HD 206165 increases significantly between N = 3,000 and N = 10,000, but shows a slower increment beyond this point. This behaviour implies that the quantile-based estimation (here q_0.25) stabilises as more of the relevant model space is sampled. Despite this stabilisation, we acknowledge that this method is not statistically rigorous and is sensitive to the grid’s local density. Therefore, this analysis should be interpreted as an indicative rather than definitive error estimate. In future implementations, we plan to incorporate more robust uncertainty estimation techniques, such as Bayesian inference or profile likelihood approaches, to properly account for both model degeneracies and observational uncertainties.

In addition to the numerical analysis, we include a figure (Fig. B.1) showing the variation of the total χ² value as a function of the number of selected models N. As expected, the χ² value tends to increase with larger N, since more models with poorer fits are progressively included in the sample. This behaviour reinforces the idea that, while selecting a higher number of models may provide broader estimates of parameter uncertainties, it also introduces a growing contribution from suboptimal fits, which may bias the interpretation. Therefore, there is a compromise between uncertainty estimation and fitting quality, which must be considered when choosing the value of N.

To further investigate the diagnostic capabilities of individual spectral lines, we carried out an additional test using helium and silicon lines for the stars HD 206165 and HD 14633. In this controlled simulation, the wind-related parameters were held fixed at their best-fit values, and the effective temperature parameter was independently derived by fitting either only the helium lines or only the silicon lines.

This test enabled us to assess the constraining power of each group of lines on key parameters such as the effective temperature. For HD 206165, the silicon-only fit reproduced the same effective temperature as the full-line analysis (20 000 K), while the helium-only fit yielded a slightly lower value, differing by about 500 K, which is within the estimated uncertainty. A similar outcome was obtained for HD 14633: the silicon-only fit recovered the global solution (34 000 K), and the helium-only fit returned a temperature of 33 000 K, i.e. a difference of 1 000 K that also lies within the uncertainty range. These results confirm the robustness of the fitting procedure across different luminosity classes, from supergiants to dwarfs.

Although these results confirm the consistency of our method, implementing a weighted χ² fitting approach that accounts for the varying sensitivity of spectral lines to specific parameters would require a broader set of diagnostics, including multiple ionisation stages. While such an extension lies beyond the scope of this study, it will be the focus of future work.

Fig. B.1 Total χ2 value as a function of the number of selected models N used to estimate uncertainties. As more models are included, the total χ2 increases steadily due to the incorporation of progressively less optimal fits.

Appendix C Information about the observed stars

The information about the spectrographs, signal-to-noise (SNR), and observation dates is summarised in Table C.1.

Appendix D Spectral line fittings for the dwarf stars

We present the comparisons between the observed spectra and the best-fitting ISOSCELES synthetic spectra for the dwarf stars in our sample. These figures illustrate the quality of the fits and complement the discussion in Section 4.1.

Fig. D.1 Observed line profiles of HD 14633 are shown in solid black. The best fast and δ-slow solutions from ISOSCELES are plotted in blue and red, respectively.

Fig. D.2 Same as Fig. D.1, but showing the observed and synthetic spectra for HD 35299.

Fig. D.3 Same as Fig. D.1, but showing the observed and synthetic spectra for HD 35912.

All Tables

All Figures

	Fig. 1 Location of ISOSCELES models in Teff–log g plane (red dots). Grey solid lines represent the evolutionary tracks from 2 M⊙ to 60 M⊙ without rotation (Ekström et al. 2012). Solid black lines correspond to both zero-age main sequence (ZAMS) and terminal-age main sequence (TAMS). Red (empty) circles indicate non-converged models in ISOSCELES.
In the text

	Fig. 2 Histogram showing rcrit values from converged solutions. Blue bars represent fast solutions, and red bars correspond to δ-slow solutions. Both axes are on a logarithmic scale.
In the text

	Fig. 3 Hα spectral-line fits for O-type star HD 14947. The black line represents the observed spectrum, while the best fit from the ISOSCELES database with a fast solution and δ-slow solution are shown in blue and red, respectively. The fast regime provides a better fit to this line profile.
In the text

	Fig. 4 Best spectral-line fits for HD 14947. The black line depicts the observed spectrum. The best fit with a fast solution is shown in blue, while the best fit with a δ-slow solution is in red. Clearly, the fast solution better describes the wind of this early-type star.
In the text

	Fig. 5 Velocity profiles as function of inverse radial coordinate u for the best-fit solutions of HD 14947. The velocity profile for the fast solution is shown in blue, while the δ-slow solution is depicted in red. Additionally, the green line represents the β-law velocity field, based on the parameters from Bouret et al. (2012).
In the text

	Fig. 6 Same as Fig. 3, but for HD 206165. Here, the δ-slow solution is the one that best fits the observed line profile
In the text

	Fig. 7 Same as Fig. 4, but for HD 206165. Both solutions show the same behaviour, demonstrating that these lines are not as sensitive to the mass-loss rate and the velocity profile as Hα.
In the text

	Fig. 8 Same as Fig. 5, but for HD 206165. The β velocity profile is from the analysis of Markova & Puls (2008).
In the text

	Fig. 9 Same as Fig. 3, but for HD 53138. In this case, the δ-slow solution is the one that best fits the observed line profile.
In the text

	Fig. 10 Best spectral-line fits for HD 53138. The black line depicts the observed spectrum. The best fit with a fast solution is shown in blue, while the best fit with a δ-slow solution is in red.
In the text

	Fig. 11 Same as Fig. 5, but for HD 53138. The β velocity profile is from the analysis of Haucke et al. (2018).
In the text

	Fig. B.1 Total χ2 value as a function of the number of selected models N used to estimate uncertainties. As more models are included, the total χ2 increases steadily due to the incorporation of progressively less optimal fits.
In the text

	Fig. D.1 Observed line profiles of HD 14633 are shown in solid black. The best fast and δ-slow solutions from ISOSCELES are plotted in blue and red, respectively.
In the text

	Fig. D.2 Same as Fig. D.1, but showing the observed and synthetic spectra for HD 35299.
In the text

	Fig. D.3 Same as Fig. D.1, but showing the observed and synthetic spectra for HD 35912.
In the text