© ESO 2021

1 Introduction

Massive stars (M_ini > 8 M_⊙) have a profound influence on their surroundings, whether chemically through their ejecta enriched with heavy elements or physically with ionising radiation, momentum, and kinetic energy (Chiosi & Maeder 1986; Matteucci 2008; Geen et al. 2020). Their strong outflows shape the interstellar medium and can play an integral role in triggering star formation (Mac Low et al. 2005; Chen et al. 2007; Elmegreen 2012; Kennicutt & Evans 2012). Massive stars are also proxies of the first generations of stars and progenitors of a range of exotic phenomena from supernovae to black holes (Haiman & Loeb 1997; Zaroubi et al. 2012). The breadth of the influence of massive stars in astronomy is difficult to overstate and an understanding of their evolution is imperative. In this context we focus on one of the most critical factors in massive star evolution, the mass loss through stellar winds (Langer 2012).

The source of mass loss in these stars is their strong line-driven winds (see e.g. Puls et al. 2008), which are inherently unstable due to the line deshadowing instability (LDI; MacGregor et al. 1979; Carlberg 1980; Owocki & Rybicki 1984). Full time-dependent, radiation-hydrodynamic wind models show the LDI leading to the development of a highly structured wind morphology consisting of out-flowing over-dense clumps, where wind opacity can increase, separated by gaps of tenuous material, which allows more photons to escape (see e.g. Owocki et al. 1988; Feldmeier 1995; Dessart & Owocki 2003; Sundqvist & Owocki 2013; Sundqvist et al. 2018; Driessen et al. 2019).

Most modern spectroscopic models of hot stars approximately include the effects of the wind instability by assuming that all the wind mass is concentrated in optically thin clumps, occupying a fraction of the total volume of the wind (e.g. CMFGEN, Hillier & Miller 1998 or PoWR, Gräfener et al. 2002). But if the clumps start to become optically thick there will be additional light leakage effects, due to the Doppler shift, through porous channels in velocity space. In particular, there can be significant light leakage for spectral lines formed in the wind (e.g. Oskinova et al. 2007; Sundqvist et al. 2010; Šurlan et al. 2013). This effect is unaccounted for in such a volume filling factor approach. In this paper we employ the unified, non-local-thermodynamic-equilibrium (NLTE), steady-state stellar atmosphere and wind modelling code FASTWIND (Santolaya-Rey et al. 1997, and follow-up papers as specified in Sect. 3), where the effect of clumps of arbitrary optical thickness has been fully implemented (Sundqvist & Puls 2018). Due to the flexibility and speed of FASTWIND, we are able to examine effects of velocity porosity (and a potentially non-void ‘interclump medium’) upon strategic diagnostics in a systematic way that has not been possible before.

To model the inherent physical processes of line transitions, stellar atmosphere models solve the stellar wind physics in a NLTE atmosphere computation and the production of the synthetic spectrum is handled by the formal solution of the radiative transfer equation. This is one of the key differences between the wind porosity included in FASTWIND and other tools used in studies such as Oskinova et al. (2007) and Šurlan et al. (2013). For example, these other methods look only at what is the impact of the wind opacity and porosity on spectral lines in the formal integral. We are able to more realistically account for these effects with FASTWIND, which includes wind opacity and porosity in the computation of the stellar and wind structure, for example, how the modified radiation field affects temperature and ionisation balances throughout the stellar atmosphere and wind.

To derive good empirical constraints on the various parameters describing the wind, spectra that include both recombination and resonance wind lines are required, each of which probe different aspects of stellar wind physics. This highlights the importance of both optical and ultraviolet coverage. In the optical, various lines produced by recombination processes are observed, while in the ultraviolet, lines are often produced by resonance scattering. Rates of recombination processes are proportionalto the density squared, implying that the strength of a spectral line resulting from a recombination process is shaped strongly by the assumed overdensities of the wind clumps in the inner wind region. On the other hand, resonance scattering is proportional to density, meaning these diagnostics generally react differently to a clumped and potentially porous wind. As a result, lines arising through each process are affected uniquely by the different wind parameters.

Currently, mass loss is one of the most significant sources of uncertainty in the evolution of massive stars, exemplified by the well-established discrepancy between theoretical predictions of mass loss and those derived from observations (see e.g. Crowther et al. 2002; Hillier et al. 2003; Bouret et al. 2005; Fullerton et al. 2006; Puls et al. 2006; Mokiem et al. 2007; Tramper et al. 2014; Bestenlehner et al. 2014). Another current problem is that no simultaneous fits to optical and ultraviolet spectra, including the P V λλ1118–1128 line, can be achieved without considering the wind structure (e.g. Oskinova et al. 2007; Sundqvist et al. 2010, 2011; Šurlan et al. 2013). These problems are strong motivations for an empirical study of the wind structure and mass loss of massive stars. Here, we focus on the end of the main sequence, with a sample of stars close to the terminal-age main sequence (TAMS). At this stage, the future evolution of the star is highly sensitive to the mass loss; it can help define the post-main-sequence evolution and the stars’ eventual end fate. This sample also represents typical rotation rates of O-type stars; coupled with the evolutionary status, this provides a good opportunity to study CNO enrichment on the main sequence, as detailed in Bouret et al. (2012). We may also obtain clues of inhomogeneities from other physical processes (rotation, stellar pulsations, etc.) by comparing theoretical and observational wind structure properties.

The layout of the paper is as follows. Section 2 describes the stellar sample and observations used. Section 3 covers our modelling software, procedure, and assumptions. Sections 4 and 5 give a generalised overview of stellar and wind diagnostics, respectively. In Sect. 6, we present the results of our analysis, including stellar and wind parameters. Section 7 is a discussion of our findings, especially wind parameters. Finally our conclusions are outlined in Sect. 8. Appendices A and B include an in-depth discussion of the analysis of each object and figures of final best fits obtained, respectively.

2 Sample and observations

The sample in the present paper is the full sample of eight O supergiants adopted from Bouret et al. (2012, hereafter BHL12). Almost all of these objects are spectral standards for their respective class (Sota et al. 2011) and their classification has not changed in follow-up studies. As this is a pilot study of wind porosity, the adopted sample presents the advantage of minimising the influence of binarity on the spectral analysis. While the majority of massive stars are found in binary systems (see e.g. Sana et al. 2013), if we are to include the effects of optically thick clumping for the first time, in the spectral fitting of a representative sample, it is best to ensure there are as few other processes as possible impacting the winds of these stars.

Table 1 shows an overview of the instruments utilised to collect the data analysed in this study while Table 2 shows an overview of the sample used in this study. All objects have IUE and FUSE spectra. BHL12 used observations in the far-ultraviolet (FUV) obtained by Copernicus for HD 66811 (ζ Puppis) andso these observations are also used here (further discussion on these observations can be found in Morton & Underhill 1977; Pauldrach et al. 1994). Northern objects have spectra from ELODIE, except λ Cep, which was observed with NARVAL. The southern objects with FEROS spectra are HD 163758 and HD 66811. For further details on the observations and data reduction we refer to BHL12. It is important to note that for this analysis we use the spectra as they are presented in BHL12, the normalisation is not adjusted.

3 Methods

Our analysis is performed with synthetic spectra produced by the stellar atmosphere and wind modelling code FASTWIND (Santolaya-Rey et al. 1997; Puls et al. 2005; Rivero González et al. 2011; Carneiro et al. 2016, specifically the version presented inSundqvist & Puls 2018). FASTWIND solves NLTE number density rate equations, assuming a steady-state spherically symmetric envelope of material encompassing the photosphere and wind of the star, where wind-inhomogeneities (clumping) are treated in a statistical way. FASTWIND accounts for thousands of spectral lines from various metals, using detailed model atoms (H, He, C, N, O, Si & P), with full co-moving frame radiation transfer for the elements mentioned here. A simplified approach is applied for other elements to account for the effects of line-blanketing while minimising runtime (see Puls et al. 2005).

As detailed in Santolaya-Rey et al. (1997), FASTWIND uses a unified approach to compute the complete atmosphere of the star by connecting a hydrostatic photosphere to an outflowing wind at a transition point defined at 10% of the sound speed. The wind follows a β-law velocity (see e.g. Eq. (2) in Sundqvist & Puls 2018), while the density structure is derived from the equation of continuity for a pre-defined mass-loss rate, Ṁ. The user inputs include a number of standard photospheric parameters: effective temperature T_eff (defined at Rosseland optical depth τ = 2∕3), surface gravity log g, stellar radius R (R_eff), and chemical abundances. If the outflows are assumed to be homogeneous, we define: Ṁ, β, and the terminal wind speed v_∞. For clumped winds a number of extra inputs are required, these are detailed in the following section (Sect. 3.1), where the number of parameters is larger for optically thick clumping than for optically thin clumping.

A microturbulent velocity (v_turb) is included with a depth-dependent prescription. During the NLTE computation, microturbulence is a fixed constant, in our case v_mic = 15 km s⁻¹. In the formal integral we select the maximum value between v_turb = 15 km s⁻¹ and v_turb = 0.1v, where v is the wind speed. This essentially acts as a fixed value of microturbulence at the photosphere, which then increases linearly with wind speed in the outer, wind-dominated atmosphere. BHL12 also implement a similar depth-dependent microturbulence, of the form v_turb = v_min + (v_max − v_min)v(r)∕v_∞, where v_min = 15 km s⁻¹ and v_max is of the order of 0.1v_∞. This inclusion of supersonic turbulence is physically motivated as a way to account for velocity dispersion in the wind. Such velocity dispersion is predicted by LDI simulations and observed as saturation over larger wavelength regions than predicted with subsonic turbulence in UV P-Cygni troughs (Hamann 1981; Lucy 1982; Puls et al. 1993; Sundqvist et al. 2011). In Sect. 3.3, we perform a robustness check to investigate what effect changing the minimum microturbulence has on the spectrum, but we do not change our microturbulence description from what is defined in this paragraph in any of our fits.

One significant difference between our study and BHL12 is that we do not include the effects of shock-generated X-ray emission. The impact of the inclusion of X-rays in our models is discussed further in Sect. 3.4. In any spectral fitting with optically thin clumping we therefore take the same assumptions as BHL12, with the exception of X-rays. Another key difference is the fitting method employed; BHL12 optimised models by eye while we use a genetic fitting algorithm, this is discussed in Sect. 3.2.

3.1 Wind clumping formalism

The standard clumping formalism is to include only the effects of optically thin clumps, parameterised by a generalised clump volume filling factor f_vol and onset velocity v_cl. The optically thin clumping formalisation is detailed in, for example, BHL12. The volume filling factor f_vol generally acts as the inverse of the clumping factor f_cl, discussed below, and becomes exactly 1/f_vol in the limit of optically thin clumping and a void interclump medium.

The adopted version of FASTWIND (FW v10.3, Sundqvist & Puls 2018) includes the leakage of light, where the clumping is parameterised with five main factors: (i) a clumping factor f_cl, (ii) a velocity filling factor f_vel, (iii) the interclump density factor f_ic, (iv) the clumping onset velocity v_cl and (v) the velocity at which maximum clumping is reached v_cl,max. For the assumption of optically thin clumping only f_cl, v_cl, and v_cl,max are required¹. The wind can be thought of as a two-component medium consisting of a mixture of dense clumps and tenuous interclump material. The material in the wind is compressed into clumps, described through the clumping factor, f_cl. The discreet clumps occupy a fraction of the total volume, this volume is the aforementioned clump volume filling factor f_vol (≈ 1∕f_cl in the case of a small f_ic). These factors are defined in relation to the mass density ρ of the mean wind as (1)

Here, the average mass density conserves mass relative to a smooth medium with the same mass-loss rate. The density of the material in the clumps is ρ_cl, and ρ_ic is the density of the tenuous interclump medium. The impact of f_cl on observed line profiles is discussed further in Sect. 5.1. The velocity filling factor can be thought of as an equivalent scaling in velocity space (2)

where δv is the velocity span in a clump and Δv is the velocity separation between clumps. The impact of this parameter on observed line profiles is discussed further in Sect. 5.3. The interclump density factor (3)

defines the density of the tenuous interclump material.

An effective opacity formalism is used to approximate this two-component medium. This effective opacity χ_eff is given as (4)

where τ_cl is the clump optical depth, which takes different forms for continuum and line absorption (see Eqs. (9) and (11) in Sundqvist & Puls 2018). The interclump density f_ic is defined in Eq. (3). The mean opacity ⟨χ⟩ is defined as ⟨χ⟩ = σn∕f_cl, where σ and n are the atomic cross-section and occupation number density (corrected for stimulated emission) of a fiducial clump, respectively. This mean opacity is computed first, using the NLTE occupation numbers of fiducial clumps with densities ρ_f = ⟨ρ⟩f_cl = ρ_smf_cl where the ‘smooth’ density ρ_sm = Ṁ∕4πv(r)r² comes from the input mass-loss rate, assuming a β velocity law (see Eq. (2) in Sundqvist & Puls 2018).

Using these general expressions one can relate back to the more commonly used case of optically thin clumps by considering appropriate limits: In the case τ_cl ≪ 1, and assuming a void interclump medium, f_ic = 0, one recovers f_cl = 1∕f_vol. Moreover, Eq. (4) for the effective opacity then yields χ_eff = ⟨χ⟩. In this limiting case, the description for the mean opacities is consistent with previous implementations of optically thin clumping (e.g. CMFGEN, Hillier & Miller 1998 or PoWR, Gräfener et al. 2002). In the case of opacities that linearly depend on density: (5)

and in the case of opacities that depend on the square of the density: (6)

The above parametrisation then aims to capture the overall effect of the time-dependent structure expected from the LDI while implementing a time-independent constant mass loss and a parameterised form for the clumping based on a set of input parameters. The time-dependent nature of the LDI is visually represented by Fig. 2 in Sundqvist & Owocki (2013) and Fig. 3 in Driessen et al. (2019).

3.2 Fitting procedure: genetic algorithm

To perform a global fit to each atmosphere of our sample of stars we utilised a genetic algorithm (GA) wrapper for FASTWIND, adapted from similar codes (Mokiem et al. 2005; Tramper et al. 2014) originally based on the pikaia framework (Charbonneau 1995). This algorithm is designed to effectively explore a large parameter space and avoid local minima, which is ideally suited to attempt to simultaneously fit for all the stellar parameters that we are able to determine spectroscopically using FASTWIND. The GA works analogously to genetic evolution and natural selection. We begin with an initial population randomly sampling the full parameter space, the ‘genes’ of which are encoded with the input parameters for FASTWIND. The models are then ranked by a fitness metric F defined as: (7)

where N is the number of spectral lines in the fit and is the reduced chi-square of each line (Mokiem et al. 2005). We compute the χ² for each individual line by comparing the modelled line profile to the observed spectrum. In this study we weigh all lines equally.

Within a given generation, models with a higher fitness are more likely to be selected to pair with another model to then produce two offspring models. These offspring models comprise the next generation and the process is repeated iteratively until convergence to an optimised fitness. To ensure the parameter space is thoroughly explored we introduce random mutations to digits of the offspring genes. These mutations are applied at a variable rate. The idea being that after a number of generations the parameter space exploration will be refined and narrowed around the models with high fitness; at this point the mutation rate is increased to re-explore the parameter space. Error margins associated with the GA exploration are described in Sect. 6. The GA has been used to fit a multitude of stellar types and datasets (see e.g. Mokiem et al. 2006, 2007; Tramper et al. 2011, 2014; Ramírez-Agudelo et al. 2017; Abdul-Masih et al. 2019). The version we use is presented in Abdul-Masih et al. (2021).

The updates that we make to the GA for this study include allowing clumps to become optically thick (Sect. 3.1), a modified list of diagnostic lines (Table 3) and a bolometric luminosity anchor (Sect. 4.1). We also tailor the population size and generations as necessary to converge on a best fit with suitable confidence intervals. For a population with 250 individuals we iterate for 200–300 generations, yielding a grand total of around 50 000–75 000 FASTWIND models per fit.

The key motivation of this work is the inclusion of optically thick clumping. To isolate the effects of a modified effective wind opacity and porosity in our analysis we split the fitting procedure into two stages. The initial fit strictly follows the prescription of optically thin clumps. In this case we aim to constrain the photospheric properties of these stars. The line list over which we define our global quality of fit can be seen in Table 3. We do not include the P V λλ1118–1128 lines when computing models with optically thin clumping due to known issues in replicating this line with optically thin clumps (BHL12), which we also reproduced when fitting models assuming optically thin clumps. These issues are discussed in Sect. 7.5. For the second round of fitting we extend our wind prescription in the atmosphere modelling to allow for optically thick clumping. For clarity we present only the analysis with optically thick clumping in the main text and all discussion of fits with optically thin clumping with FASTWIND are included in Appendix A.

The commonly used ‘by-eye’ method for massive-star spectroscopic fitting relies on optimising models to fit a number of specific diagnostics to determine individual parameters, for example iterating the surface gravity of the model to fit the wings of Balmer lines. While we are not using such a direct method to constrain individual parameters, as the exploration of the parameter space is automated based on the fitness ranking of the models in a population, we ensure such diagnostic lines are included. The general line list used for most GA runs is shown in Table 3. If, for an individual object, any line is omitted we discuss this on a case-by-case basis in Appendix A.

If used blindly, such automated fitting is not guaranteed to ensure each parameter is constrained to its individual true value. The spectra must be manually prepared and fitting regions carefully selected. As mentioned previously we do not adjust the normalisation of these spectra from BHL12 but we do remove any blends (that are not included in the model) within our selected fitting regions from the spectra as these would negatively influence the fitness metric. Figures of all line regions used for fitting are provided in Appendix B and indicate the diagnostic lines used as well as the omitted spectral ranges. We also monitor the fitness line-by-line for each parameter to investigate the sensitivity of each line to individual parameters. Therefore, with careful manual consideration of the input observations and physical shortcomings of input models, the GA provides an automated method to converge on unique, reproduceable solutions with a well defined statistical foundation for evaluation. The overall agreement between GA-based fitting and a by-eye approach has been investigated by Markova et al. (2020) and was found to be very good.

Fig. 1 Effect of the inclusion of X-rays on the resonance line N V λ1239–1243. GA best fit model without X-rays – red solid line, best fit with X-rays – blue dashed line. We note that N V λ1239–1243 is contaminated with Lyman-α in the blue part of the spectral line, which causes the discrepancy in width of the absorption trough between the spectrum and model.

3.3 Robustness: minimum microturbulence

When varying minimum microturbulence (v_turb as defined earlier in Sect. 3) in these stars the largest effect is on photospheric metal lines and weak helium lines: an increase in v_turb will increase the strength of these lines (McErlean et al. 1998). The effect on strong resonance lines in the UV is marginal. Therefore, in these models the microturbulence stands only to impact abundance diagnostics. As abundance determination is not a main focus of this study we leave more extensive analysis of the microturbulence prescription to future studies.

Fig. 2 Effect of changing effective temperature on diagnostic lines He I λ4471 and He II λ4541. GA best-fit model – red solid line, best fit with 1kK lower Teff – blue dashed line, best fit with 1 kK higher Teff – green dotted line.

3.4 Robustness: X-rays

The inclusion of X-rays, as currently parameterised in FASTWIND, does not have a significant effect on any line profiles investigated in this study (see Table 3 for the list of line profiles). Therefore, we do not include X-rays and all best-fit parameters were obtained using models without X-rays. For this robustness check we manually run an individual FASTWIND model, shown in Fig. 1. In this model we include X-ray parameters from Abdul-Masih et al. (2019), based on recommendations in Carneiro et al. (2016). These recommendations are tailored for typical O stars, meaning for this robustness check we have not examined potential effects of higher X-ray luminosities than those normally assumed for O stars. On the other hand, for lines from a high ionisation stage such as N V λ1239–1243, X-rays have a significant effect on the profile. Included in Fig. 1 is the high ionisation stage UV line N V λ1239–1243. For this line, and others of high ionisation stage, X-rays are required to reproduce the line profile, these lines are excluded from this analysis.

4 Stellar parameter diagnostics

4.1 Radius

The stellar radius is an important FASTWIND input. We adjust the model radius to balance with the varying temperatureso as to maintain a fixed luminosity throughout the GA run. For this analysis we opted to use a bolometric luminosity anchor, taking luminosities calculated in BHL12 and adjusting such that in any model the radius followed , where L is the bolometric luminosity, σ_SB is the Stefan-Boltzmann constant, and the effective temperature is a free parameter.

4.2 Effective temperature

Effective temperature is often constrained by the ionisation balance of available elements. For example, if multiple helium absorption lines are present in the spectrum from different energy levels and ions, comparing their relative strengths (equivalent widths) can be a robust temperature estimate. In our χ² based approach, we rely on the absolute strength of all diagnostics lines, including (ideally) He, O and N lines of different ionisation stages (Table 3) and the model sensitivity in the parameter regime. Figure 2 shows the effect of changing T_eff on diagnostic helium lines.

4.3 Surface gravity

The surface gravity is typically constrained using the hydrogen Balmer series, which is affected by Stark broadening. This shapes the wings of these profiles and is strongly dependent on the electron density at the surface of the star. Figure 3 shows the effect of changing logg on Hγ. Mostly, Balmer transitions from higher levels (Hδ and Hγ) are uncontaminated by wind infilling emission and are ideal surface gravity diagnostics. Hβ remains dominated by an absorption profile in most objects but the wind emission begins to fill in the line core, which can also affect the line wings, in objects with strong winds.

The effective temperature also affects the electron pressure so there is some degeneracy between temperature and surface gravity. To help break this degeneracy our list of diagnostic lines (Table 3) includes temperature diagnostics that are not as sensitive to log g, as discussed in Sect. 4.2. The present GA based approach allows us to take these correlations into account while converging towards the optimal set of atmospheric and wind parameters.

Fig. 3 Effect of changing surface gravity on Hγ. GA best-fit model – red solid line, best fit with 0.2 lower log g – blue dashed line, best fit with 0.2 higher log g – green dotted line.

4.4 Surface abundances

Determining surface abundances is not one of the main focuses of this study; therefore, we do not include an extensive line list. This means we do not rigorously examine the abundances, rather we include a subset of lines and compare with what is found in BHL12 where the surface abundances were a focal point. We include multiple lines from N, C and O in both the optical and UV (Table 3). We aim to include at least two different ionisation stages from each element, and for each ionisation stage a photospheric line profile relatively free from contamination from other processes like winds or blends. However, for Si we only include one ion so we will not discuss the Si abundance further in our analysis.

4.5 Rotation and turbulence

After the formal integral the synthetic spectra are broadened with an instrumental, rotational and macroturbulent profile. The rotational and macroturbulent velocities are free parameters in the GA, and metal absorption lines are good diagnostics for rotational and macroturbulent broadening. Figures 4 and 5 show the effect of changing rotation and macroturbulence on the diagnostic line O III λ5592, respectively.Most objects are included in Simón-Díaz & Herrero (2014), who provided an in-depth look at the contributions from macroturbulence and rotation in the overall line broadening. Generally we find that our confidence intervals include the values determined in Simón-Díaz & Herrero (2014), this is discussed on a case-by-case basis in Appendix A.

Fig. 4 Effect of changing rotation on the broadening of O III λ5592. GA best-fit model – red solid line, best fit with lower rotation – blue dashed line, best fit with higher rotation – green dotted line.

Fig. 5 Effect of changing macroturbulent velocity on the broadening of O III λ5592. GA best-fit model – red solid line, best fit with lower macroturbulence – blue dashed line, best fit with higher macroturbulence – green dotted line. Note that the effect is small in this case but increases for lower v sin i.

4.6 Radial velocities

We find radial velocities for these stars by fitting a 1D Gaussian profile to optical helium lines in the observed spectra and by measuring the offset between the mean of the Gaussian and the line centre wavelength used in FASTWIND. This offset correction is then applied to all FASTWIND models produced by the GA before the fitness of the model is evaluated.

5 Wind parameter diagnostics

When fitting with optically thin clumping we allow the mass-loss rate, β, and terminal wind speed to vary and use a fixed clumping factor corresponding to those found in BHL12. When running the GA with optically thick clumping the number of free parameters is increased. We include velocity filling factor, interclump density factor, and the clumping onset velocity alongside all other parameters. We also test varying a subset of parameters while fixing others to the values found in BHL12. Generally this leads to less accurate fitting, either due to systematic differences between CMFGEN (used by BHL12) and FASTWIND or due to optimisation on a case-by-case basis resulting from the difference in fitting techniques.

Fig. 6 Effect of changing mass-loss rate and clumping factor on the Hα recombination line and on the Si IV λλ1393–1402 resonance line, using the best-fit model including optically thick clumping. GA best-fit model – red solid line. Best fit with increased clumping and fixed mass-loss rate – blue dashed line. Best fit with reduced clumping and fixed mass-loss rate – yellow dotted line. Best fit with increased clumping and reduced mass-loss rate – magenta dashed line. Best fit with reduced clumping and increased mass-loss rate – green dotted line.

5.1 Mass-loss rate, clumping, and terminal wind velocity

In the optical, the mass-loss rate and clumping factor are highly degenerate. Fortunately, they affect resonance lines and recombination lines differently. The recombination lines are particularly sensitive to the clumping factor as the resultant increase in recombination rates is proportional to density squared and thus the ionisation balance changes. Clumping also has an impact on the optical depth, affecting the radiative transfer. The resonance lines are less sensitive as this process is proportional to density. This means the mass-loss rate should be constrained more so by the absorption components of resonance lines while the clumping factor is constrained by the recombination lines.

To illustrate this, Fig. 6 shows FASTWIND models for the Hα line with various combinations of clumping factor and mass-loss rate, applied to our best-fit model with optically thick clumping. In general mass-loss rates and clumping factors are degenerate in recombination lines; such that, if the best-fit model has a clumping factor of 20 and a mass-loss rate log(Ṁ) = −5.8, and we were to reduce the clumping factor to 10 and increase the mass-loss rate accordingly², we would produce a FASTWIND model that is statistically indistinguishable from our best fit within a reasonable error margin. This degeneracy is broken by the simultaneous fitting of resonance lines, which are less sensitive to the clumping factor (Fig. 6 bottom panel).

The terminal wind speed is usually determined by the extent of the bluest edge of the absorption in the saturated P-Cygni profiles of the UV resonance lines, of which C IV λλ1548–1550 is often a particularly good example (Prinja et al. 1990). This value can be determined almost independently from any other if one observes a saturated P-Cygni line such that the absorption flux is essentially zero up until the blue edge of the absorption trough. One should consider the microturbulence when fitting the blue edge of the P-Cygni line as we include a depth-dependent microturbulence v_turb defined in Sect. 3; however, the effect here is marginal, as discussed in Sect. 3.3.

Fig. 7 Effect of changing β on a recombination diagnostic line Hα. GA best-fit model – red, best fit with decreased β – blue, best fit with increased β – green.

Fig. 8 Effect of changing fvel on P V λλ1118–1128 resonance lines. GA best-fit model – red solid line, best fit with 0.3 higher fvel – blue dashed line, best fit with 0.5 higher fvel – green dotted line.

5.2 Beta wind parameter and clumping onset

These wind parameters control the wind acceleration and the location where clumping starts in the wind. A low β means the wind accelerates quickly thus casting an array of effects on line profiles, depending on which wind region dominates the line formation. In lines that are formed close to the stellar surface, for example Hα, a low β causes a significant reduction in emission due to the reduction in wind density that results from the high velocity gradient close to the wind onset. The clumping onset velocity sets the point at which clumps develop in the wind: if the clumping onset velocity is high the clumps begin to develop further away from the star. Within our formalism the clumping factor then increases linearly to the maximum clumping factor. For this analysis, we fix the velocity at which the wind reaches maximum clumpiness to twice the onset velocity. Figure 7 shows the effect of varying β on Hα.

5.3 Velocity filling factor and interclump density

These factors act as light blockers, as they increase so does the probability of light interacting with material in velocity space and between clumps, respectively. These parameters are constrained using the P-Cygni troughs in UV resonance lines. As f_vel and f_ic are reduced the troughs become less saturated. This effect is shown in Fig. 8.

Fig. 9 χ2 distributionper parameter for the stellar and wind parameters and chemical surface abundances of HD 14947. The shaded region shows the confidence interval, marking out the upper and lower bounds surrounding the best fitting model for each parameter. The colour scheme evolves with the generations to give one a sense of the convergence of the algorithm, the darker blue points are the earliest generations and lighter blue points correspond to the latest models. The red points show the model grids run afterwards to constrain errors. This figure is continued on the next page.

6 Results

For each object, we present best-fit model optically thick clumping parameters and further photospheric and wind parameters as found by the GA in Tables 4 and 5, respectively. In Table 5, there are two rows per object; these correspond to the best fit found by BHL12, who use optically thin clumping, and by the GA, which allows for optically thick clumping. We also produce a number of different GA best fits using different clumping assumptions and fitting strategies. These are discussed in Appendix A and presented in Table B.1. The best-fit model presented here is the model with the minimum global , where is the squared difference between model and data for the individual lines. The GA returns approximately 60 000 models in order to explore the parameter space. An example of these model distributions is shown in Fig. 9. For an in-depth analysis of the fit quality per star see the Appendix A.

For half the sample, our fits that allow for optically thick clumping agree with BHL12s effective temperature within 0.5 kK. For HD 210839 we find a temperature 1 kK higher than BHL12. For the three remaining objects (HD 16691, HD 190429A and HD 14947) we find clearly different solutions, with 1.5 kK lower than BHL12, 3 kK higher and 3.5 kK higher. These stars are all O4.5 or earlier but these results are not pointing towards systematic differences in fitting, rather these are different fitting solutions than found in BHL12. These are discussed in the Appendix A. Our surface gravity values again present no systematic differences, although for the two stars of latest type we find significantly larger log g values and higher β. Also, our surface gravity for HD 16691 is high, but this also appears too high in the best fit, there are a number of difficulties in fitting this object discussed in Appendix A.

The mass-loss rates are remarkably consistent with BHL12, with the largest difference being 0.1 dex, again in the two latest objects. In all objects we agree with BHL12 within 0.1 dex in mass-loss rate, despite including the effects of optically thick clumping and often finding significantly different clumping factors. We note that BHL12 had to reduce the phosphorus abundance in their models to find these mass-loss rates while we do not. We obtain a systematically higher broadening than BHL12 when considering the sum of rotation and macroturbulent broadening applied to each object. This can be explained by the implementation of macroturbulent broadening. In BHL12 an isotropic prescription is used while we implement a radial-tangential broadening. Three of the objects are included in Simón-Díaz & Herrero (2014), the broadening values found by these authors using a Fourier method are consistent with those found from the GA, the largest discrepancies found in v sin i and v_mac are 5% and 20%, respectively.

In terms of abundances, we find similar trends as BHL12, confirming the evolutionary status of these stars. The only clear difference is the difficulty in determining an oxygen abundance for ζ Puppis due to low signal-to-noise ratio (S/N) in the optical oxygen diagnostic. Abundances measurements are less consistent in the fits withoptically thin clumping, due to a combination of mainly using diagnostic lines dominated by wind processes and an insufficient wind prescription in the fitting. This shows that if one attempts to determine surface abundances using wind line profiles the wind prescription has a significant effect.

We find higher clumping onset velocities than BHL12, which may result from the different clumping laws used. This is discussed in detail in Sect. 7.4. For velocity filling, we find an average value throughout the sample of 0.44± 0.33 and an average interclump density of 0.13 ± 0.08. These findings are discussed in detail in Sect. 7.2.

The confidence intervals on the best-fit parameters are calculated as in Tramper et al. (2011, 2014), Ramírez-Agudelo et al. (2017) and Abdul-Masih et al. (2019): all χ² values are normalised such that the best fitting model has and the error bounds include all models within the 95% confidence interval and are considered statistically equivalent. In this scheme, the exact confidence boundaries depend on the sampling around the peak. Another usual limitation of this approach is the assumption that the best-fit model is a good replication of the observed spectrum.

To overcome the possibility of underestimated errors we run an orthogonal grid of models around the best fit, varying each free parameter independently, to complete any confidence interval coverage the GA may have missed. This method works well for parameters that have a fairly consistent effect on the global fitness and correlate little with other fitting parameters; for example, an abundance variation will have the same effect on the fitness of any relevant metal lines and other lines will be fairly unaffected, therefore giving a smooth fitness distribution around the global best fit. On the other hand, effective temperature has various effects on fitness depending on the line, an increase in temperature could well increase the fit to the Balmer series but reduce fitness to low ionisation helium and metal lines, thus the distribution of fitness of models around the best fit are less smooth. The errors on temperature quoted here are then extremely sensitive to the quality of fit, the line list used, and the weight of each line. We note also that these errors are statistical formal errors from the GA; the actual errors may be larger due to a number of effects. These effects include: using fixed luminosities, which affects the radii and thus mass-loss rates, potential shortcomings in the spectra such as normalisation, which has an impact on log g, and the limited line list used in the GA. In the case that errors remain unrealistically small we adopt a minimum error region, these are identifiable in Tables 4 and 5 as any symmetric error values around the best-fit value.

7 Discussion

7.1 Mass-loss rates and clumping factors

It is now well established that neglecting clumping may lead to erroneous estimates of mass-loss rates. Up until now, mostly models including optically thin clumping have been used to constrain spectroscopic mass-loss rates. The effect of optically thick clumping on various line diagnostics has been studied but the resultant systematic effects on empirical mass-loss rates across a sample of early-type stars remains unclear. It is thought that the additional leakage of light through porous channels between clumps, allowed by optically thick clumping, may have a significant effect on the derived mass-loss rate.

Overall our mass-loss rates are very similar to BHL12. The largest difference of 0.08 dex is only found in two stars while the average difference between the studies is 0.03 dex. For the present sample it seems that the mass-loss rates found with optically thick clumping and optically thin clumping do not diverge significantly. However, this is in the case of using a reduced phosphorus abundance to fit the unsaturated PV lines, in the rates determined using optically thin clumping. If a solar phosphorus abundance was used to determine the mass-loss rates with optically thin clumping we would see a systematic difference compared to the rates found with optically thick clumping. Instead, discrepancies are due more to goodness-of-fit on a case by case basis.

We find, for example, for HD 14947 that Vink et al. (2000) over-predicts the mass-loss rate by a factor of 3.9 compared to our results, which account for optically thick clumping, while the Björklund et al. (2021) rates (hereafter the Leuven rates) predict values within 10% of what is found by the GA. Such a reduction has been observed in a number of UV spectroscopic studies (Crowther et al. 2002; Hillier et al. 2003; Bouret et al. 2005; Fullerton et al. 2006), and also in mass-loss rates determined from X-ray diagnostics (Cohen et al. 2014). Most of the mass-loss rates found with the GA sit slightly above the Leuven predictions, with the Leuven rates predicting a mass-loss rate 0.8 times of that observed, on average (Fig. 10). The Vink et al. (2000) rates predict a mass loss that is on average a factor of 3.6 above the observed values. A similar trend can be seen in the modified wind-momentum (Fig. 10). It is worth noting that clumping is not the source of discrepancy between these two theoretical mass loss predictions: the effect of clumping is prominent when comparing observational diagnostics to theoretical predictions, whereas the differences between these theoretical predictions come from the different techniques used to calculate the radiation force and solving the equation of motion (see discussions in Sundqvist et al. 2019; Björklund et al. 2021).

Throughout the sample we find clumping factors in the range of 20–32 for any wind prescription with an average of clumping factor of 25,while BHL12 finds values in the range 17–33 with an average clumping factor of 22. This is consistent with a radially average clumping factor found in the LDI simulations by Driessen et al. (2019) for an O-supergiant with ζ Puppis like input parameters. We note that these are maximum clumping factors. There is clear evidence that the clumping factor may decrease in the outer wind (Puls et al. 2006; Najarro et al. 2011); however, it is difficult to diagnose this with the spectral coverage used for this study. Therefore, we use a near-constant clumping law and leave a further investigation into the effect of clumping stratification on observational diagnostics to a future study with infrared spectral coverage.

Fig. 10 Mass-loss rates (upper panel) and modified wind-momentum (lower panel) found by the GA best fits, compared to the predictions made by Leuven (Björklund et al. 2021) and Vink (Vink et al. 2000) for each star in our sample. The GA mass-loss rates are those found from the best fit with optically thick clumping. Black cross in bottom right shows the minimum error on mass-loss rate and wind momentum ± 0.05, from this analysis, and error on luminosity ± 0.1 from BHL12.

7.2 Velocity filling factors

For five stars in the sample the GA best fits have velocity filling factors in the range 0.2–0.6. These are the first systematic empirical constraints derived for this parameter. The exceptions are HD 192639 with a lower velocity filling and HD 16691/HD 163758, which show uncharacteristically high velocity filling factors, bordering on recovering a velocity profile reminiscent of a smooth wind outflow. On average, we find a velocity filling factor of 0.44 for the full sample. If the three stars that show uncharacteristically high or low velocity filling factors are excluded the average is reduced slightly to 0.32. Theoretical constraints for the velocity filling factor can be made using LDI simulations from Driessen et al. (2019). We calculate a clump velocity span and clump separation from a prototypical O-supergiant model using region cutoffs at v∕v_∞ = 0.05 or 1.1 R_* and consider a clump to be any region with a density higher than 2⟨ ρ⟩. These limits predict a velocity filling factor f_vel = 0.56, which is overall consistent with the empirical constraints. Altering these limits will influence the velocity filling estimate. Therefore, within reasonable density limits 2 − 5⟨ρ⟩ we find an estimate of 0.57 > f_vel > 0.44, as computed from the velocity span of the clumps. For the stars that do not agree with theoretical estimates we additionally attempt to fit the velocity filling parameters with fixed photospheric parameters, varying only the optically thick wind parameters. In these tests, it is possible to find fits of similar quality to the UV resonance lines with lower velocity filling but, to compensate, the interclump density converges to the upper limit of our parameter space. We are certainly limited by data quality, specifically S/N in this sample and so it is difficult to distinguish whether these stars do show uncharacteristically high velocity filling (for which the physical motivation is unclear) or if a higher S/N would allow fitting with lower velocity filling.

Similar parameters have been explored in works including Oskinova et al. (2007) and Šurlan et al. (2013) although the implementation of such a parameter is quite different; these authors use a clump separation parameter in physical space, with discrete clumps each providing independent contributions to line driving and absorption. Implementation in FASTWIND is a wind averaged velocity filling, so a direct comparison between parameters is difficult. In general our results are qualitatively consistent with these studies. We find that allowing optically thick clumping in the stellar wind has a significant impact on line profiles, and the resulting enhanced porosity reduces saturation in P-Cygni profiles.

7.3 Interclump density

This parameter is hard to constrain empirically, as can be seen throughout the GA fits for this sample. We generally find flat distributions around a best χ² with confidence intervals encompassing large portions of the parameter space. This parameter has a similar, although weaker, influence on the model as the velocity filling factor. Therefore, it may be the case that we are unable to significantly distinguish between values of interclump density in a global χ² metric, mostly due to relatively low S/N in our UV observations and the fact that only a few lines are impacted, decreasing the impact of f_ic on the global χ².

With this in mind, a loose constraint on the interclump density can be estimated by looking specifically at the fitness distribution for saturated P-Cygni line profiles. It is not possible for a line to saturate if the interclump medium is completely void and f_vel is significantly lower than unity, so these lines could act as a lower limit diagnostic (Sundqvist et al. 2011). The C IV λλ1548–1550 line is mostly saturated in our sample, and this line provides confidence intervals around 0.15–0.3 on the f_ic value, indicating a fairly high density in the interclump medium. Such a value is then higher than typically seen in 1D LDI models of clumping, which indicates that the lateral ‘filling in’ of radially compressed gas, observed in 2D LDI simulations (Sundqvist et al. 2018), likely plays an important role in setting the clump to interclump density contrast. While a global fit might point towards a more tenuous interclump medium, due to the degeneracy with the velocity filling factor, it is possible that actually a higher interclump density is preferred. Perhaps it will be possible to constrain this in other P-Cygni profiles with a higher S/N in the UV but for now saturated lines point towards a relatively high interclump density. Zsargo et al. (2008) also showed that a tenuous interclump medium is required in hot stars,in the context of reproducing the features of ‘super-ionised’ species such as OVI in UV and X-ray observations.

7.4 Beta and clumping onset

We find high velocities for clumping onset relative to previous spectroscopic studies. BHL12 has an average value of 60 km s⁻¹ for this sample, and our optically thick models give an average almost five times higher at 290 km s⁻¹. However, the clumping profile is different between our study and BHL12. In BHL12, an exponential increase in clumping factor with velocity is used while we implement a linear increase. In practice, the profile is similar between both studies but the nature of the clumping profile in BHL12 means that the onset velocity is not a strict onset, but the velocity at which the clumping factor reaches a threshold value. Also, the clumping development is more gradual for higher onset velocities. An example of the clumping profiles used in CMFGEN and FASTWIND are shown in Fig. 11. For most stars, BHL12 uses a clumping onset velocity of 30 km s⁻¹ and we find higher onset velocities, meaning we still have quite low clumping factors in regions where BHL12 is approaching maximum clumping. For the fast rotators where BHL12 finds higher onset velocities our clumping profiles are very similar.

For most lines, the clumping onset is unconstrained but it is clear that high onset velocities provide higher fitness in recombination lines throughout the sample. We note, however, that BHL12 had difficulty fitting the double peaked profiles that emerge in the emission components of the fast rotators in the sample. They were able to improve their fits by increasing the clumping onset velocity by around 100 km s⁻¹. In this study we find that the fits can indeed be further improved if the onset velocity is increased by an even greater extent, around 200 km s⁻¹. On the other hand, we point out that in this work we have convolved all lines with a rotational broadening profile assuming a fixed photospheric value for the projected rotation velocity v sin i. But due to angular momentum conservation in the wind, this value should really decrease as one moves away from the stellar surface. Bouret et al. (2008) shows that if a depth-dependant rotation is applied to the models the observational profiles can be better reproduced. We have not investigated this, but our findings are inline with the predictions made by these studies and so we can suggest thatdepth-dependent effects on the rotational broadening profiles of these objects are likely to have a significant effect on some of the emergent wind profiles.

In some stars, we observe very high clumping onsets and high β values. In these cases the GA tries to fit the core and blue edge of emission profiles, to varying degrees of success, generally not very well. The high clumping onsets might not have a huge impact on photospheric diagnostics but the resulting high β values may influence photospheric parameters.

Other studies have also noted that a high β mimics the effect of clumping as a higher β leads to a denser wind near the wind onset region, having a similar effect as an increased clumping factor or mass-loss rate (Petrov et al. 2014). This certainly appears to hold true, at least if one were to consider recombination lines as shown for Hα in Fig. 7. It is clear from this study that an increased clumping onset velocity can act to compensate for this increased emission, leading to no clear systematic impact on clumping or mass loss. However, it is possible that such an offset could occur if the clumping onset was fixed to lower velocities. We note that in our fits with optically thin clumping the clumping onset velocity is fixed to 0.05v_∞, which leads to onset speeds around 100 km s⁻¹. In these fitswith clumping onset fixed to a lower velocity we do not find the same high β values as derived with optically thick clumping, but values compatible with those from BHL12 in their analysis. However, in these stars we also find somewhat different mass-loss rates than obtained when fitting with optically thick clumping. An early onset of clumping has been found both from previous diagnostic studies (Puls et al. 2006; Cohen et al. 2011) and from theoretical LDI simulations (Sundqvist & Owocki 2013). Moreover, since for Galactic O supergiants radiative acceleration might exceed gravity even in deep sub-surface layers (at the so-called ‘iron-opacity bump’ at T ≈ 200 kK), it is possible that this might trigger a very turbulent atmosphere (Jiang et al. 2015); resulting, not only in the observed large values of ‘macroturbulence’, but also affecting clumping properties near the photosphere. Future multi-dimensional LDI simulations might examine such a potential connection by proper downward extension to those deep sub-surface layers.

It is still unclear whether the high onsets seen in this sample are truly high clumping onset velocities or an observational proxy of another physical process. Regardless, it is clear from our current dataset that it is difficult to distinguish between slow wind acceleration (high β) and an earlier onset of clumping. As such, the overall result here of generally quite slow acceleration in near photospheric layers, accompanied by a quite high velocity-onset of clumping, should be re-investigated in future work.

Fig. 11 Clumping profiles for CMFGEN (BHL12) and GA fits, showing clumping profiles with the same maximum clumping factor but varying the onset velocity. The velocities shown correspond to best-fit onset velocities found for HD 14947 (slow rotator – BHL12:vcl = 30, GA:vcl = 150) and HD 210839 (fast rotator – BHL12:vcl = 120, GA:vcl = 250).

7.5 Phosphorus

A common problem in spectral analyses of stars in this regime is the phosphorus abundance. Using models with optically thin clumping, up until this point, it has only been possible to reproduce the P-Cygni absorption by reducing the phosphorus abundance to sub-solar levels. This was first noticed by Pauldrach et al. (1994) and has been reobserved in a number of subsequent studies (see e.g. Pauldrach et al. 2001; Crowther et al. 2002; Hillier et al. 2003; Bouret et al. 2005). The phosphorus abundance reduction is generally thought to be unphysical as phosphorus should not be processed to this extent throughout the stellar lifetime, this is discussed in detail in Fullerton et al. (2006). To address this we exclude the phosphorus abundance as a parameter in our fitting, opting instead to keep phosphorus fixed to solar abundance. In a test run, we tried to fit the phosphorus abundance in the optically thin clumping case and indeed reproduced the problems discussed above. Therefore, we do not try to fit this line with optically thin clumping, we only include the phosphorus lines when modelling optically thick clumping to see if the added wind porosity in velocity space helps to reproduce the lines. Other studies (e.g. Oskinova et al. 2007; Sundqvist et al. 2010, 2011; Šurlan et al. 2013) have been able to reproduce these lines with solar phosphorus abundance by including the effects of wind velocity porosity, and we obtain the same result here, albeit using a quantitative fitting approach. Throughout our sample we are able to consistently fit the UV phosphorus resonance doublet as a result of our inclusion of wind velocity porosity without lowering the phosphorus abundance. An example of this effect on the emergent line profile is shown in Fig. 8.

8 Conclusions

We performed automated spectroscopic fitting of stellar atmosphere models using a genetic algorithm for a sample of O supergiants, simultaneously across optical and UV wavelengths, to investigate the effect of optically thick clumping and wind porosity on their wind properties. Our results are:

• We are able to resolve the PV problem by including optically thick clumping. By including the effects of wind porosity in velocity space it is possible to reproduce the phosphorus line profiles with solar phosphorus abundance;

• If the PV lines are excluded from the analysis, optically thick clumping does not significantly affect the mass loss measurements relative to those measured with optically thin clumping. Nonetheless, we stress the general importance of a multi-diagnostic approach to constraint this mass loss, where clumping must be considered in order to obtain consistent line-profile fits across the optical and UV ranges;

• The resulting empirical mass-loss rates are well constrained by our GA fitting. Overall, they agree well with the new theoretical predictions by Björklund et al. (2021) and are on average a factor of 3.6 lower than the models by Vink et al. (2000). Since the latter are often used as standard in evolutionary calculations, this strengthens earlier claims (e.g. Sundqvist et al. 2011, 2019; Najarro et al. 2011; Šurlan et al. 2013; Cohen et al. 2014; Krtička & Kubát 2017; Keszthelyi et al. 2017; Björklund et al. 2021) that line-driven mass-loss rates included in current models of massive-starevolution should be reduced. Although the empirical study here needs to be extended to a larger sample of stars (and the theoretical models by Björklund to a larger grid), our combined results suggest that simply scaling down the earlier rates by a factor of 3 is a quite reasonable first approximation for the considered spectral types;

• From our systematic study, we are able to derive empirical constraints on wind structure parameters associated with the effects of optically thick clumping. On average, we find that clumps cover roughly half of the wind velocity field and that inter-clump densities are around 10–30% of the mean wind density. While such clump velocity filling factors agree well with 1D LDI simulations of clumping, the rather high interclump densities suggest that the lateral filling-in of radially compressed material seen in corresponding 2D models might be critical for setting the density scales of the rarefied interclump medium. Regarding the clumping factor, we derive mean maximum values on order 20, which again agrees well with the results of current 1D LDI simulations for Galactic O supergiants;

• We notice our best fit models often have significantly higher clumping-onset speeds than what was found in previous studies; however, we stress both the importance of future studies accounting for a depth-dependent rotational broadening profile and current degeneracy-issues between the onset of clumping and the wind acceleration parameter β in near-photospheric layers;

We intend to extend this type of study to a wider range of spectral types across the upper Hertzsprung-Russell diagram and to lower metallicity environments. This method can also be applied within the framework of the ULLYSES and XShootU programmes which will provide large datasets of massive stars in the Magellanic Clouds.

Acknowledgements

This project has received funding from the KU Leuven Research Council (grant C16/17/007: MAESTRO), the FWO through a FWO junior postdoctoral fellowship (No. 12ZY520N) as well as the European Space Agency (ESA) through the Belgian Federal Science Policy Office (BELSPO). The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation – Flanders (FWO) and the Flemish Government. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

Appendix A Star-by-star and optically thin versus thick clumping

In this appendix we discuss the quality of the best fits found by the GA for each object in detail, highlighting specific shortcomings and discussing potential causes. We also introduce the initial round of fitting, for which we run the GA including the assumption that all clumps remain optically thin.

In the fits with optically thin clumping we use the same clumping profile as we use for optically thick clumping, diverging from BHL12 in the clumping law as discussed in Sect. 7.4. We also opt to fix the clumping consistently throughout the optically thin clumping fits, highlighting issues that may arise if a general assumption for the clumping factor is used and this parameter is not optimised. In BHL12 a clumping factor of 20 is used for 5 objects and so we use this value. Again, this is not optimised to test the effect and in 3 objects (HD14947, HD190429Aand HD16691) BHL12 adapts the clumping factor. In BHL12 the clumping onset velocity is fixed to 30 kms⁻¹ (and increased for 3 objects: HD16691, HD66811 and HD210839), whereas we opt to fix our clumping onset speed to 0.05v_∞, which results in an average clumping onset speed of 110 kms⁻¹. HD210839 and HD163758 display some effects of a clumping structure that is not optimised. Here, the surface gravities are significantly lower than what would be determined using Hγ more or less alone; the low surface gravity comes as the GA attempts to fit Hα and settles on a lower value as it cannot shape the profile with clumping. If one were to vary the clumping factor a more realistic surface gravity can be found. Another way to get around this is to increase the weight of the Balmer lines Hγ and Hδ to constrain surface gravity. In a sense, these optically thin clumping ‘best fits’ are not true best fits, but the best that can be achieved with the input assumptions and caveats. However, the optically thick clumping best fits allow more of these parameters to be free and result in more optimised fits.

We have found temperature discrepancies larger than 1kK (which is more than twice our minimum error) compared to BHL12 for three stars in the sample with similar helium and metal line strengths. This may indicate a difference between CMFGEN and FASTWIND as seen in Massey et al. (2013). However, it could also be that the GA converges on an alternate solution and other parameters are adjusted accordingly to compensate. It is difficult to make generalised statements across the sample when comparing to BHL12, a more complete picture is gained by looking at the fits to each object individually. We also notice that when fitting with optically thin clumping, we sometimes find significantly different surface abundances than found with optically thick clumping. This shows that the clumping assumption used can heavily influence abundance determination when fitting line profiles shaped by the winds.

A.1 HD16691

This is a difficult star to do a fit comparison to BHL12 as the GA best fit with optically thin clumping is particularly bad, largely due to difficulties reproducing Hβ and Hα. The inclusion of optically thick clumping was also unable to reproduce Hβ so we removed it from the global fit with optically thick clumping presented here. In this case, these lines drive the optically thin clumping fit to a vsini = 250 kms⁻¹ and vsini = 85 kms⁻¹. Such a high value is not supported by the metal absorption lines, which are all clearly overbroadened. The broadening parameters in the optically thin clumping fit are not consistent between UV and optical, it would be helpful to carry out a more extensive analysis of this star in the future.

Despite the difficulties in fitting, when comparing optically thin wind clumping fits between the two studies we find consistent values for T_eff, log g, β, v_∞, and surface abundances. We find a significantly lower mass-loss rate than BHL12, by a margin of 0.2 dex (four times the minimum error). A difference was indeed expected as we use a higher clumping factor than BHL12. The cause of the difference could lie in the clumping prescription. Given similar onset values are found in both fits, the clumping profile used in BHL12 has a much more gradual development. With a more rapid onset of clumping, here we require a lower mass-loss rate to reproduce the spectrum. We obtain noticeably poorer fits to He I λ4471 and He II λ4541 (our main temperature diagnostics), but this is due to the aforementioned overbroadening and helps to explain the discrepancy in helium abundance. We find significantly different best-fit parameters when allowing optically thick clumping. This is particularly notable as the maximum fitness found is significantly higher than that found with optically thin clumping, despite a peculiar velocity filling factor that borders on replicating a smooth wind. Our T_eff is 1.5kK lower than BHL12 (three times larger than minimum error) and log g 0.2 higher (four times minimum error). Comparing the relevant diagnostics, both studies are able to reproduce optical HeI lines while BHL12 has a more accurate optical HeII strength. The optically thick clumping finds a particularly good fit to the strengths of Hδ and Hγ compared to BHL12; however, the wing fits may be slightly better in BHL12, excluding the nitrogen line blend in Hδ, which BHL12 often has difficulty fitting. We note also that there is a sharp drop in fitness for both Hδ and Hγ at β lower than 1.6, while a β consistent with BHL12 is preferred by wind emission lines Hα and He II λ4686. The emission lines also drive a high clumping onset velocity, three times that found by BHL12, showcasing the issues discussed in Sect. 7.4.

We obtain a worryingly high velocity filling factor with a reasonable interclump density in our best fit. These parameters show no clear fitness distribution peaks in most lines except for N IV λ1718, which clearly prefers a high velocity filling. In another GA run we fixed the photosphere properties to those found in the optically thin clumping run and vary only the wind parameters. In this run we see a more reasonable velocity filling and an overestimated interclump density. In this run the velocity filling converges on a more physically motivated value of 0.6; however, the fitness distributions per line do not change considerably, nor does the best-fit spectrum. It is possible that a lower velocity filling could be compensating for the fixed higher temperature and lower log g. As a result, the interclump density is increased to fit the UV P-Cygni troughs, highlighting the degeneracy between these wind parameters.

In a further test run, we fix all photospheric parameters to those found by the GA optically thin fit and vary only wind parameters. In this test we allow for optically thick clumping and fit an even smaller subset of lines, only those heavily affected by the stellar wind. In this case our mass loss, β, and clumping onset velocity match BHL12, but f_cl = 10. We also find a reasonable f_vel = 0.5 but with an extremely high f_ic = 0.3, which is only limited by our upper parameter boundary. This test has a very similar quality of fit as found by the GA with optically thick clumping. With all of this is mind, and the peculiar Hβ profile, it is clear that this star has a unique wind profile worth further investigation.

A.2 HD66811 - ζ Puppis

HD66811 appears to have two solutions for T_eff and log g: after finding 41kK and 3.5 with optically thin clumping and 40kK and 3.5 with optically thick clumping, another GArun, which allows for optically thick clumping, was submitted that returned 41kK and 3.6. Comparing these to the 40kK and 3.6 found by BHL12, a safe conclusion is to say the solution lies within T_eff =40-41kK and logg =3.5-3.6, this fits comfortably within the errors estimated by BHL12 but is beyond the margin provided by our perhaps too optimistic minimum errors.

Most parameters agree between BHL12 and the GA with optically thick clumping, apart from the clumping onset velocity and some of the abundances. The onset velocity is again heavily restricted to higher values by Hα and He II λ4686, this is discussed in Sect. 7.4. When looking at all the GA runs, each abundance found by BHL12 is found at least once. The carbon abundance is matched by the fit with optically thin clumping (run 1), the nitrogen matched by the run with optically thick clumping, which finds a low T_eff and log g (run 3), and the oxygen abundance matches all but run 3. This is showing that there is an issue when not including enough lines in the GA to consistently constrain the abundances. For all GA runs we find a significantly higher helium abundance than BHL12; however, the run with optically thin clumping is close enough to cover the discrepancy within the error margin and this run has the best fit to He I λ4471 so indeed an abundance close to BHL12 is the more realistic and likely solution.

A.3 HD190429A

We find a significantly higher T_eff than BHL12 when fitting with optically thick clumping, with a consistent log g. Even with a large lower boundary margin on our best fit the temperature is still 3kK higher. The high temperature is being constrained by the optical hydrogen lines. Generally, the GA is able to reproduce these lines more accurately than BHL12; however, the He I λ4471 profile is too weak in the GA best fit and it is clear that if we were to focus on this line the temperature would be lowered. In our fit with optically thin clumping we see a slightly better fit to this profile with T_eff only 1kK higher than BHL12.

When fitting with optically thick clumping we match BHL12 for all other parameters except clumping factor and onset velocity. We find a difference of a factor of 1.2 in clumping factor and an onset velocity 50 kms⁻¹ higher than BHL12. Although, the differences in these parameters between studies are not so large considering the difference in clumping laws used.

The GA fit with optically thin clumping is clearly a poorer fit than that with optically thick clumping for all lines except for HeI lines. However, with the lower T_eff found with optically thin clumping also comes a low surface gravity, 0.2 dex lower than the other best fits. Another issue in the optically thin clumping fit for this star is the abundance measurements. As mentioned previously the subset of lines used here result in incorrect abundances when a limited prescription for the wind physics is used as, in the UV, the strength of the lines are generally dominated by the stellar wind. Therefore, when attempting to fit without a good enough wind implementation the abundances vary unrealistically in order to reproduce these profiles.

A.4 HD15570

In the best fit with optically thick clumping we match BHL12 for nearly all parameters within our error margins. There are also no clear differences in the best-fit spectra for the line list we use between the two studies, except the strength of N IV λ4058 appears to be overestimated in BHL12. As for parameter discrepancies, our β is larger by 0.1 but our clumping onset velocity is 200 kms⁻¹ compared to 30 kms⁻¹ in BHL12, again highlighting the issue discussed in Sect 7.4. Our macroturbulence is higher than BHL12 but this is to be expected given the difference between the isotropic and radial-tangential prescriptions. The discrepancy of 400 kms⁻¹ in terminal wind speed is interesting as both the GA and BHL12 achieve roughly the same quality of fit to the bluest edge of the P-Cygni trough in C IV λλ1548-1550.

Most of the best-fit parameters returned by the GA with optically thin clumping are not significantly different than with optically thick clumping, with the exception of the surface gravity and abundances. These issues have also been discussed in other objects. A low surface gravity can be limited by an attempt to fit emission features while using a fixed clumping prescription. Also, abundances can vary in an attempt to fit UV metal lines shaped by the winds when an insufficient wind prescription is used in the models. A way out of this dilemma would be to allow clumping to vary, and either allow optically thick clumping or add more photospheric metal lines into the fit.

A.5 HD14947

We find one of the largest discrepancies in temperature between the GA and BHL12 in this star, our best fit with optically thick clumping settling at 40.5kK. This is 3kK higher than BHL12, which is especially odd as the quality of fits of the main T_eff diagnostics between the two models are remarkably similar. This discrepancy is certainly aided by the offset in surface gravity; we find a log g almost 0.2 dex higher than BHL12 and again a by-eye comparison of the models in each study reveals little as both best fits replicate the wings of Balmer profiles very well. Abundances, broadening, mass loss, and β are all in good agreement. We find a significantly lower clumping factor, a difference of a factor of 1.5, but again the best-fit models do not appear to be significantly different. BHL12 finds a terminal wind speed 300 kms⁻¹ higher than the GA but the GA best fit does appear to fall short of the P-Cygni edge so the higher velocity is the more realistic result.

The best fit with optically thin clumping agrees well with the other best fits. The optically thin clumping best fit finds some middle ground in T_eff and log g with a surface gravity to match BHL12 and a temperature 1.5kK higher than BHL12. We notice again an issue in surface abundances with our optically thin clumping method, but in this case the differences are not too concerning.

A.6 HD210839 - λ Cep

Comparing T_eff and log g between BHL12 and the GA with optically thick clumping, we notice discrepancies slightly larger than can be accounted for by our errors; however, the two solutions are not too far apart, the T_eff certainly lies between 36-37kK and log g appears to be slightly weighted to lower values by optical recombination emission lines, so a higher value around 3.55 found by BHL12 is more likely the better solution. All other parameters agree quite well, our helium abundance is slightly higher and strongly weighted by the Balmer lines, and the GA does replicate the strengths of the Balmer series quite well here. A noteworthy feature of this spectrum is the prominent double peak shown in the Hα and He II λ4686 lines, this leads the GA to settle on a high clumping onset velocity as discussed in Sect. 7.4.

Our optically thin fit once again suffers some discrepancies due to our clumping prescription, a log g lower by 0.2 dex is found as well as a β value close to 0.6 with these lower values being preferred by the double peak optical recombination profiles, especially He II λ4686. The abundances found in our optically thin run are significantly lower than those found in BHL12, a common problem in this study with limited line list but again the problem is alleviated when including optically thick clumping.

The interclump density for this object is an intermediate value of 0.16, which is fairly common in this sample. The velocity filling is quite low but reasonable, driven by the PV lines, which have a weaker absorption depth than some of the other objects in this sample. The physical cause of the reduction in PV absorption depth is unclear.

A.7 HD163758

Bestfit temperatures agree comfortably within error regions for all fits, comparing the fit with optically thick clumping to BHL12 we have a slightly lower T_eff. Our fit with optically thick clumping finds a log g consistent with BHL12, the fit quality is very good in both studies, both have some trouble fitting Hδ due to a strong blend. We inspect the fitness for just the Hγ line, and it is clear that a slightly lower surface gravity is allowed, so the true value is somewhere between the bounds provided by both fits, therefore log g ≈ 3.4-3.55. Wefind a mass-loss rate consistent with BHL12 with a clumping factor 1.5 times higher, although a difference of 1.5 times in clumping factor only accounts for a 0.05 dex shift in mass-loss rate. From the models one can see that, in the GA, Hα is possibly underestimated in an attempt to fit the double-peak profile while He II λ4686 is replicated well. However, in BHL12 the strength of the synthetic Hα is closer to the observed spectrum but He II λ4686 is overestimated. We also find one of the highest clumping onset velocities in this star of 320 kms⁻¹, while BHL12 used 30 kms⁻¹ in their fitting. Our high onset causes β to converge to a high value of 2.3. This problem is discussed in Sect 7.4 and certainly can influence the mass loss.

Rotation rates and macroturbulence are fairly consistent in both fits. We find nitrogen and oxygen abundances that match BHL12 but our carbon abundance is lower. BHL12 also takes care to fit optical carbon lines so the carbon abundance of 8.6 estimated by BHL12 is likely more accurate. We also find a helium abundance significantly higher than BHL12 of 0.28, this is mainly driven by He II λ4686. Difficulty fitting this line likely has a negative impact on our helium abundance. Evidence for this comes from the fact that this line drives a high β value in our fit with optically thick clumping, while in our fit with optically thin clumping fit β remains low. In our fit with optically thin clumping the helium abundance matches BHL12 and He II λ4686 is poorly reproduced.

Our fit with optically thin clumping is very similar to that found in BHL12, with most best-fit parameters in agreement.

A.8 HD192639

For this object, we find a similar T_eff, 0.3 dex higher log g, 0.07 dex difference in mass-loss rate and a very high β and clumping onset velocity. Rotation rates are similar while we find a significantly higher macroturbulence. Carbon, nitrogen and oxygen abundances are in agreement and we find a slightly higher helium abundance. An additional problem present in our fitting of this object is that we find a much too high terminal wind speed, 800 kms⁻¹ higher than BHL12. The blue edge of the P-Cygni trough is not so well defined in the spectrum due to the iron forest. If we were to fit the spectrum by eye, we would aim to fit the point of bluest complete absorption as opposed to the point at which the edge meets the continuum, while the GA appears to fulfil the latter criteria. It is clear from our fit to Si IV λλ1393 − 1402 that our terminal wind speed is too high, this line reaches a maximum fitness at v_∞ ≃ 2100kms⁻¹, which is a much more consistent value for the full UV spectrum.

Appendix B Best fits

In this appendix we present models from our fits with optically thick clumping, showing only the line profiles used in our GA analysis. The observed spectrum is shown by the black points, the solid red line is our best fitting model, and green lines represent any models generated during the GA iterations that lie within the error regions.

Fig. B.1 Fitness distribution for HD14947 fit with optically thick clumping. The colour scheme of the points goes fromblue to red, corresponding to the earliest and latest models in the exploration, respectively. The region highlighted in red shows the confidence interval. This figure is included to demonstrate the good agreement between an assessment of fit quality by the fitness metric F as defined in Eq. 4 and χ2 as shown in Fig. 9.

Fig. B.2 Best fit for HD16691 from GA with optically thick clumping.

Fig. B.3 Best fit for HD66811 from GA with optically thick clumping.

Fig. B.4 Best fit for HD190429A from GA with optically thick clumping.

Fig. B.5 Best fit for HD15570 from GA with optically thick clumping.

Fig. B.6 Best fit for HD14947 from GA with optically thick clumping.

Fig. B.7 Best fit for HD210839 from GA with optically thick clumping.

Fig. B.8 Best fit for HD163758 from GA with optically thick clumping.

Fig. B.9 Best fit for HD192639 from GA with optically thick clumping.

References

All Tables

All Figures

	Fig. 1 Effect of the inclusion of X-rays on the resonance line N V λ1239–1243. GA best fit model without X-rays – red solid line, best fit with X-rays – blue dashed line. We note that N V λ1239–1243 is contaminated with Lyman-α in the blue part of the spectral line, which causes the discrepancy in width of the absorption trough between the spectrum and model.
In the text

	Fig. 2 Effect of changing effective temperature on diagnostic lines He I λ4471 and He II λ4541. GA best-fit model – red solid line, best fit with 1kK lower Teff – blue dashed line, best fit with 1 kK higher Teff – green dotted line.
In the text

	Fig. 3 Effect of changing surface gravity on Hγ. GA best-fit model – red solid line, best fit with 0.2 lower log g – blue dashed line, best fit with 0.2 higher log g – green dotted line.
In the text

	Fig. 4 Effect of changing rotation on the broadening of O III λ5592. GA best-fit model – red solid line, best fit with lower rotation – blue dashed line, best fit with higher rotation – green dotted line.
In the text

	Fig. 5 Effect of changing macroturbulent velocity on the broadening of O III λ5592. GA best-fit model – red solid line, best fit with lower macroturbulence – blue dashed line, best fit with higher macroturbulence – green dotted line. Note that the effect is small in this case but increases for lower v sin i.
In the text

Fig. 6 Effect of changing mass-loss rate and clumping factor on the Hα recombination line and on the Si IV λλ1393–1402 resonance line, using the best-fit model including optically thick clumping. GA best-fit model – red solid line. Best fit with increased clumping and fixed mass-loss rate – blue dashed line. Best fit with reduced clumping and fixed mass-loss rate – yellow dotted line. Best fit with increased clumping and reduced mass-loss rate – magenta dashed line. Best fit with reduced clumping and increased mass-loss rate – green dotted line.

In the text

	Fig. 7 Effect of changing β on a recombination diagnostic line Hα. GA best-fit model – red, best fit with decreased β – blue, best fit with increased β – green.
In the text

	Fig. 8 Effect of changing fvel on P V λλ1118–1128 resonance lines. GA best-fit model – red solid line, best fit with 0.3 higher fvel – blue dashed line, best fit with 0.5 higher fvel – green dotted line.
In the text

Fig. 9 χ2 distributionper parameter for the stellar and wind parameters and chemical surface abundances of HD 14947. The shaded region shows the confidence interval, marking out the upper and lower bounds surrounding the best fitting model for each parameter. The colour scheme evolves with the generations to give one a sense of the convergence of the algorithm, the darker blue points are the earliest generations and lighter blue points correspond to the latest models. The red points show the model grids run afterwards to constrain errors. This figure is continued on the next page.

In the text

Fig. 10 Mass-loss rates (upper panel) and modified wind-momentum (lower panel) found by the GA best fits, compared to the predictions made by Leuven (Björklund et al. 2021) and Vink (Vink et al. 2000) for each star in our sample. The GA mass-loss rates are those found from the best fit with optically thick clumping. Black cross in bottom right shows the minimum error on mass-loss rate and wind momentum ± 0.05, from this analysis, and error on luminosity ± 0.1 from BHL12.

In the text

	Fig. 11 Clumping profiles for CMFGEN (BHL12) and GA fits, showing clumping profiles with the same maximum clumping factor but varying the onset velocity. The velocities shown correspond to best-fit onset velocities found for HD 14947 (slow rotator – BHL12:vcl = 30, GA:vcl = 150) and HD 210839 (fast rotator – BHL12:vcl = 120, GA:vcl = 250).
In the text

Fig. B.1 Fitness distribution for HD14947 fit with optically thick clumping. The colour scheme of the points goes fromblue to red, corresponding to the earliest and latest models in the exploration, respectively. The region highlighted in red shows the confidence interval. This figure is included to demonstrate the good agreement between an assessment of fit quality by the fitness metric F as defined in Eq. 4 and χ2 as shown in Fig. 9.

In the text

	Fig. B.2 Best fit for HD16691 from GA with optically thick clumping.
In the text

	Fig. B.3 Best fit for HD66811 from GA with optically thick clumping.
In the text

	Fig. B.4 Best fit for HD190429A from GA with optically thick clumping.
In the text

	Fig. B.5 Best fit for HD15570 from GA with optically thick clumping.
In the text

	Fig. B.6 Best fit for HD14947 from GA with optically thick clumping.
In the text

	Fig. B.7 Best fit for HD210839 from GA with optically thick clumping.
In the text

	Fig. B.8 Best fit for HD163758 from GA with optically thick clumping.
In the text

	Fig. B.9 Best fit for HD192639 from GA with optically thick clumping.
In the text