© The Authors 2026

1. Introduction

O-type stars are key to understanding the evolution and fate of massive stars. These stellar objects, defined by the presence of He I and He II lines in their spectra (Morgan et al. 1943 and Sect. 4 in Maíz Apellániz et al. 2026 for a recent review) are located in the upper main-sequence region of the Hertzsprung-Russell diagram, corresponding to initial masses above ∼15 M_⊙ (Humphreys 1978; Martins et al. 2005; Holgado et al. 2025). As they represent the earliest and longest-lived phase of massive stars, it is essential to understand their evolution during this phase. This information is crucial to the study of a range of astrophysical phenomena in which they play a key role (e.g., core-collapse supernovae, gravitational waves, and galactic and extragalactic evolution; Langer 2012; Mandel & Farmer 2022; Massey 2013; Eldridge & Stanway 2022).

Nucleosynthesis is the main driver of stellar evolution. In massive main-sequence stars, the core nuclear processes are thought to be well understood, with hydrogen burning through the CNO cycle as the main energy source. As a result, the He and N abundances increase in the stellar core during this phase, while carbon and – at later stages – oxygen are depleted. However, how these internal abundance changes are reflected at the stellar surface remains an open question even for effectively single stars, as the nature and efficiency of mixing mechanisms are still not fully understood (Martins & Palacios 2013; Johnston 2021; Keszthelyi et al. 2022). In the last decades of the last century, some observational studies began to report surface overabundances of helium and nitrogen in OB stars, challenging the available schema of stellar evolution (e.g., Herrero et al. 1992; Gies & Lambert 1992; Schonberner et al. 1988).

One key process proposed to explain He and N enrichment of the stellar surfaces in O-type stars during the main sequence was rotationally induced mixing (e.g., Maeder & Meynet 1987; Meynet & Maeder 2000). In this framework, stellar rotation can give rise to instabilities that exchange material between the core and the surface, modifying the core size and composition while contaminating the stellar atmosphere with nucleosynthesis material. The adoption of physical prescriptions such as the advecto-diffusive transport (e.g., Maeder & Zahn 1998), magnetic dynamo (e.g., Spruit 2002), or internal gravity waves (e.g., Talon & Charbonnel 2003) in the different models and different mixing efficiencies coming from different calibrations lead to different predictions. In this regard, Martins & Palacios (2013) presented a systematic comparison of the outcomes from different state-of-the-art evolutionary models, including Brott et al. (2011) and Ekström et al. (2012).

In recent years, several studies have reported a high incidence of multiplicity among Galactic massive stars (e.g., Chini et al. 2012; Sana et al. 2012; Sota et al. 2014; Barbá et al. 2017). Others have claimed that binary interactions have a major role in their evolution, revealing a more complex scenario. In addition to the single-star physics, several processes involving the interaction of various components in a multiple system can take place. For example, tidal forces can regulate the spin rate and affect rotational mixing (Petrovic et al. 2005; de Mink et al. 2009, 2013; Song et al. 2013; Sciarini et al. 2026). Moreover, mass transfer can affect the surface abundances of both the donor, which may form a stripped star (Paczyński 1967; Kippenhahn 1969; Vanbeveren et al. 1998), and the gainer, which may become enriched in helium and nitrogen (Langer et al. 2008, 2020; Martínez-Sebastián et al. 2025; Jin et al. 2026). In fact, Bolton & Rogers (1978) first proposed mass transfer as the physical origin of ON stars (O stars with enhanced N absorption lines and deficient carbon and oxygen features; Walborn 1971; Sota et al. 2011). Despite subsequent observational studies having found evidence supporting binarity playing a major role in this spectral type (Boyajian et al. 2005; Li & Howarth 2020), their origin remains unclear (Martins et al. 2015b).

In this increasingly complex scenario, observations are essential to constraining stellar evolutionary models and theories. For this reason, the massive-star community has undertaken large spectroscopic surveys covering different regions and metallicities, such as FSMS (Evans et al. 2005), GOSSS (Maíz Apellániz et al. 2011), VFTS (Evans et al. 2011), MiMeS (Wade et al. 2016), IACOB (Simón-Díaz et al. 2011, 2020), OWN (Barbá et al. 2010, 2017), XShootU (Vink et al. 2023), BLOeM (Shenar et al. 2024), and upcoming projects as WEAVE-SCIP (Jin et al. 2024) and 4MIDABLE-LR (Chiappini et al. 2019). Although large samples are required to trace the evolutionary pathways of massive stars, such surveys inevitably comprise a mixture of objects, namely, some that effectively evolved in isolation and others that were shaped by binary interactions. It is therefore crucial to incorporate additional observational diagnostics that can help distinguish stars evolving effectively as single objects from binary-interaction products; otherwise, the inferred constraints for single-star evolutionary models may be contaminated by binary products mimicking some of the observed properties in isolated stars.

Jin et al. (2026) have shown that surface abundances – combined with other stellar parameters – constitute a powerful diagnostic to discriminate between single and binary evolution products. Certain relationships follow strict predictions from nuclear physics that must be satisfied regardless of other stellar processes. An example is the relationship between the surface enrichment of nitrogen and helium, which Martínez-Sebastián et al. (2025) used in order to argue that a group of O-type He-rich yet comparatively low N stars are likely binary products. Moreover, because chemical enrichment depends on the efficiency of internal mixing, studying the surface abundances of stars that evolved in isolation is essential to providing observational constraints on evolutionary models.

The determination of reliable abundances relies on quantitative spectroscopy in which observed spectra are compared with state-of-the-art atmosphere models. However, different approaches are subject to different limitations depending on the adopted atmospheric code and, more importantly, on the specific comparison techniques employed. To better understand the N abundance in O-type stars, Rivero González et al. (2011, 2012a,b) developed a detailed N atomic model, implemented it in FASTWIND (Santolaya-Rey et al. 1997; Puls et al. 2005), and investigated the formation of N lines in stellar atmospheres.

Several studies analyzing Galactic O-type stars (e.g., Martins et al. 2015a,b, 2017; Mahy et al. 2022) have used the stellar atmopshere code CMFGEN (Hillier & Miller 1998). Though the code provides a physically consistent treatment of the stellar atmosphere, the computation of line blanketing by solving the radiative transfer equation of all lines in the comoving frame makes CMFGEN models computationally expensive. As a consequence, the sample sizes analyzed in these studies have been effectively limited. Furthermore, the spectral synthesis methodology adopted in these works is sensitive to extrinsic broadening and is also slower than the curve-of-growth method.

Other studies have used FASTWIND. This code uses an approximation to treat line blanketing that enables faster model computation, allowing for the construction of larger and denser grids. Based on this code, Markova et al. (2018) used a visual fitting approach to determine abundances in a sample of Galactic O-type stars. This strategy implies some user-subjectivity and is also relatively time-consuming. Therefore, it is only suitable for relatively small samples (30 stars were used in their study). Carneiro et al. (2019) used a more objective semiautomatic approach based on the curve-of-growth method – comparing observed and synthetic equivalent widths (EWs) for multiple lines using an optimized grid of models. This approach allows contributions from lines not explicitly included in the synthetic spectra to be separated from the lines of interest, making it an ideal technique to use with FASTWIND. Although this approach enables the analysis of larger samples, it is limited to spectra in which reliable EWs can be measured.

All of the aforementioned studies have expanded the available information of N surface abundances in Galactic O-type stars, and some of them also include C and O abundances in large samples as well as all luminosity classes and rotational velocities in O-type stars (e.g., Martins et al. 2015a). Although their different approaches have been shown to be coherent (e.g., Markova et al. 2018), the diversity in methodologies, atmosphere codes, sample selections, and spectra makes it challenging to construct a uniform sample when combining them. Consequently, we still need to increase the number of systematically derived abundances from large, bias-controlled samples. In this context, we aimed to contribute to this effort by studying N abundances for a set of 120 Galactic O-type stars observed in the framework of IACOB project (Sect. 2). To achieve this, we applied the curve-of-growth method together with an optimized grid of FASTWIND models (Sect. 3). Our estimated N abundances are consistent with previous studies (Sect. 4), and by combining them with previously determined spectroscopic parameters and surface He abundances (Holgado et al. 2025; Simón-Díaz et al. 2026), they allowed us to assess the possible evolutionary origins of the sample (Sect. 5). We found evidence in line with a binary origin of He-rich O-type stars and a bimodal N abundance distribution within these objects. This highlights the interest of future work including fast rotators and supergiants to complement our understanding on massive-star evolution in their early stages (Sect. 6).

2. Sample

We have selected an optically curated sample in which N abundances can be reliably determined from the total 358 Galactic O-type stars for which the IACOB project has gathered high-quality spectra (last described in Holgado et al. 2025). Among them, we retained the 225 stars with a very good agreement between the observed spectra and the FASTWIND models for H and He lines. We excluded stars displaying P-Cygni profiles, double-peaked lines, or discrepancies in the fitting between Hα and He IIλ4686, to avoid possible effects of wind inhomogeneities. As modeling atmospheres with strong winds remains challenging, we restricted our sample to 210 stars with log Q ≤ −12.7 (wind-strength parameter, Q = Ṁ/(M_⊙ yr⁻¹)/(v_∞/(km s⁻¹)R/(R_⊙))^3/2; first described in Puls et al. 1996). Given the extremely low statistics resulting from the application of this criteria for stars with luminosity classes (LC) I and II (i.e., supergiants and bright giants), we decided to limit ourselves to LC III, IV, and V. We imposed two additional criteria on the signal-to-noise ratio (S/N > 50) and the apparent visual magnitude of the stars (m_V ≤ 10) to ensure the quality of the spectra. This results in a sample of 145 stars. Our sample contains SB1 systems while is curated to exclude detected SB2 binaries. Nonetheless, some stars may still harbor undetected companions¹ (Simón-Díaz et al. 2026).

After this initial selection, we conducted a visual pre-analysis of the resulting subset to ensure the availability of diagnostic lines for the determination of N abundances. We inspected two strong lines per ion (N II 3994.99 and 6610.56 Å, N III 4379.11 and 4514.86 Å, N IV 4057.76 and 6380.77 Ų, and N V 4603.73 and 4619.98 Å) to assess whether they are reliably detected in the available spectra (further details in Appendix B). Because N III lines are the most prevalent across the sample and are modeled with high reliability, we focused our abundance analysis on stars showing strong N III features (effective temperature T_eff ≤ 40 kK). Finally, to ensure an accurate EW measurement, we limited our sample to stars with low projected rotational velocities (v sin i ≤ 150 km s⁻¹).

Figure 1 depicts the distribution of the final sample of 117 stars in the spectroscopic Hertzsprung–Russell diagram (spectroscopic HR diagram or sHRD; Langer & Kudritzki 2014) with large filled symbols. The shape and color of each marker indicates the luminosity class of each star. For reference, we include the full set of 358 O-type stars analyzed in Holgado et al. (2025) and the non-rotating evolutionary tracks computed by Ekström et al. (2012) using the code GENEC³.

Fig. 1. Distribution in a spectroscopic HR diagram of the sample of Galactic O-type stars considered in this work (large filled symbols). For reference purposes, the full analyzed IACOB sample of O-type stars is included as small size circles, while red and yellow symbols are for LC I and II. The GENEC evolutionary tracks during the main sequence are plotted for vini/vcrit = 0 with solid lines. Dotted lines are used for gravities decreasing from log g = 4.2 dex (left) in steps of Δ log g = 0.2 dex.

The distribution of the analyzed stars is approximately uniform for log(ℒ/ℒ_⊙) < 4.2 and 30 ≤ T_eff ≤ 40 kK, roughly corresponding to luminosity classes III, IV, and V. Any discussion beyond this regime remains speculative and would require additional work.

3. Data and methods

The analyzed spectra have been gathered in the framework of the IACOB project, with resolving power in the range between 25 000 and 85 000 (Simón-Díaz et al. 2011, 2020). In each case, we use the spectra from the IACOB database with the best signal-to-noise ratio, with ∼80% of them having a S/N ≳ 100.

For each star, we benefited from information about its effective temperature (T_eff), surface gravity (log g), line-broadening parameters (v sin i and macroturbulent velocity, v_mac), and wind-strength parameter (log Q), and surface He abundance (Y_He = N(He)/N(H)) as provided in Holgado et al. (2025) and Simón-Díaz et al. (2026)⁴. We consider these parameters well determined and keep them fixed during the N abundance analysis. Nonetheless, we finetune some values when necessary (e.g., to correct inconsistencies in the nitrogen ionization balance).

For the estimation of the N abundance (ϵ_N = log(N/H) + 12), we follow an analysis strategy similar to Carneiro et al. (2019), based on the curve-of-growth method. We use the EWs of a suitable set of diagnostic lines (Table 1) to determine N abundances by comparing them with those predicted by the closest reference model from our FASTWIND grid (Sect. 3.2). This approach is largely insensitive to extrinsic broadening mechanisms (i.e., v sin i and v_mac) and to radial-velocity corrections. Consequently, it provides a more reliable determination of the microturbulent velocity (ξ). In addition, this method allows for a faster reanalysis once the EWs are measured.

3.1. Measurement of equivalent widths

We measured the EWs by fitting the observed lines with a modeled profile and integrating the resulting fits. We prefer this approach over direct integration because it allows us to treat blended lines separately through multiple-line fitting. Each profile is based on a Gaussian with an EW corresponding to the atomic profile and a full-width-half-maximum corresponding to the instrumental profile and a free component accounting for the intrinsic broadening (i.e., microturbulence). This is convolved with two broadening functions representing rotation and macroturbulence effects. To achieve the best possible fit, we let v sin i, v_mac, and the line center λ₀ vary as semi-free parameters. We treat the line depth and intrinsic broadening as free parameters, and fix the continuum level at unity.

For each star, we tried to measure a total of 24 N II, 12 N III, and 3 N IV lines⁵ (Table 1). We note that we never used all the listed lines, as their detection depends of the specific parameters of the stars, but we use as many as possible in each case. To account for the main source of uncertainty in the observed EW (i.e., the continuum normalization), we decided to fit each line using both the original and a locally renormalized continuum. This allowed us to estimate the EW uncertainty (ΔEW) and to select the optimal normalization as the one that minimizes the mean residual between the observed and fitted profiles.

3.2. Grid of FASTWIND models

We constructed a grid of synthetic spectra with FASTWIND v10.6.5 (Santolaya-Rey et al. 1997; Puls et al. 2005; Rivero González et al. 2012a), optimized to analyze N abundances in O-type stars (previously used in Martínez-Sebastián et al. 2025). We adopt parameter ranges similar to those of the grid described in Simón-Díaz et al. (2026), but including C, N, and O as explicit elements and only three values of Y_He⁶. We fix the exponent of the wind velocity law to β = 1 and neglect clumping – consistent to first order with our sample selection criteria. We run the models assuming solar metallicity (Asplund et al. 2009) and adopt the same atomic models as Carneiro et al. (2019), including the nitrogen model atom developed in Rivero González et al. (2011, 2012a,b). Table 2 summarizes the parameter ranges. As a sanity check, we confirm that all H I, He I and He II line profiles are consistent between our grid and that of Simón-Díaz et al. (2026).

3.3. Abundance and microturbulence determination

We determined, consistently and simultaneously, ξ and ϵ_N. We used the measured EW_l, obs to map the parameter space of ϵ_N and ξ for each line using a χ²:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{l}}^2={\left(\frac{E{W}_{\mathrm{l},\mathrm{obs}}-E{W}_{\mathrm{l},\mathrm{mod}}}{\mathrm{\Delta }E{W}_{\mathrm{l},\mathrm{obs}}}\right)}^2.\end{array}$$

To refine our mapping, we interpolated the model EWs (EW_mod) in steps of Δξ = 1 km s⁻¹ and Δϵ_N = 0.01 dex.

For each line, we selected the abundance that minimizes χ_l² for a given ξ. In Fig. 2 we show an example of the diagnostic diagram we used with all available N III lines for the O9.5 IV–V star HD 206 183 (the same star is used in subsequent figures within this section). For each ion, we then computed a weighted average minimum from the different lines with a factor w_l which average the different lines from the same multiplet to avoid giving extra weight to transitions between the same levels (∑_{l in mult}w_l = 1; dotted line in Fig. 2).

Fig. 2. Example of the mapping of N abundance as a function of microturbulence. Colored lines show the N abundance that minimizes Ml as a function of microturbulence for each of the analyzed lines in HD 206 183. The horizontal dashed line marks the ξ value determined with the N III multiplet. The colored X marker corresponds to the abundance of each line at this microturbulence. The dotted line indicates the average abundance, with the gray area representing the 1σ dispersion around this value.

We determined ξ using lines from the same multiplet whenever possible. We identified the ξ value that minimizes the abundance scatter among the lines of the same multiplet. We prefer this approach over the use of all the lines to minimize spurious results due to the scatter of different lines caused by the normalization of different regions. With this approach, we could ensure a consistent normalization over the used lines. We assumed a single ξ for all nitrogen ions, so this determination is valid for the element as a whole. In the example of Fig. 2, the resulting value is shown with the horizontal dashed line at ξ = 7 km s⁻¹. When this strategy was unfeasible, we estimated ξ through comparison with stars of similar LC and log g for which the microturbulence is well constrained.

In a subsequent step, we selected the abundance for each line corresponding to the adopted ξ (marked with X in Fig. 2), and we determined the surface N abundance as the total average of all lines. We estimated the uncertainty as the standard deviation of all studied lines, adopting a minimum value of Δϵ_N = 0.15 dex (justified in Appendix C). We then plotted the abundance of each line as a function of EW_obs to visually identify possible inconsistencies (Fig. 3). If we detected systematic discrepancies in the resulting line-by-line abundances that can be linked to a fundamental parameter (e.g., an incorrect ionization balance), we repeated the analysis after adjusting the reference model. In a final step, we assessed the quality of the analysis in a visual review by overplotting the closest synthetic spectra from the grid onto the observations for key diagnostic lines (Fig. 4).

Fig. 3. Example of the diagnostic plot used to correct possible systematic effects in the abundance determination. It shows abundances of the different analyzed lines in HD 206 183 for a fixed microturbulence velocity. Green and red symbols correspond to N II and N III lines, respectively. Solid lines in the same colors indicate the average abundance for each ion, with the shaded areas being their 1σ confidence intervals.

Fig. 4. Example of the evaluation of the N abundance fit in selected spectral regions. Orange vertical lines denote the nitrogen transitions employed in the abundance determination; dotted gray vertical lines indicate other N lines present in the model but excluded from the analysis. We note that at the temperature of this object, N III 4642 is dominated by an overlapping O II line (Rivero González et al. 2011).

We note that the parallel determination of microturbulence and N abundance in this work may lead to ξ values that differ from those derived from He in Simón-Díaz et al. (2026). To our knowledge, there is no evidence that the microturbulence associated with nitrogen and helium in O-type stars must be the same. Nevertheless, we consider the impact that ξ may have on the He abundance determination (e.g., McErlean 1998; Villamariz & Herrero 2000; Howarth & Smith 2001; Massey et al. 2013; Markova et al. 2020). Given the relevance of this effect for the following discussion, we repeated the He abundance analysis using the microturbulence derived in this work. Comparing both determinations, the results remain consistent within the uncertainties and do not affect the main conclusions of this study (Appendix C.2). Without any further clear physical argument to favor one set of microturbulence values over the other, we adopt our determinations, as they ensure internal consistency with the N abundance analysis.

4. Results

Table A.1 presents the N abundance estimates for the 117 analyzed stars, together with additional parameters relevant to our study. We also include the most up-to-date spectral classifications from GOSSS (Sota et al. 2011, 2014; Maíz Apellániz et al. 2016), as well as the runaway status from Carretero-Castrillo et al. (2023, 2026) and private communication. We organize the table into two groups according to luminosity class (IV + V and III), each sorted by spectral type. With respect to the quality of the results, we have 85 stars with a very good fitting of the diagnostic lines and 34 with a compromise solution, in which some lines would require different abundances. Nonetheless, all of them have good fittings considering uncertainties.

4.1. Comparison with the literature

This study includes 28 stars previously analyzed in the literature using different stellar atmosphere codes and methodologies. Martins et al. 2015a,b, 2017) and Mahy et al. (2022) used CMFGEN models and a spectral synthesis methodology. They assumed a microturbulence varying from 10 km s⁻¹ at the photosphere to 10% of the terminal velocity at the outer boundary of the atmosphere. Aschenbrenner et al. (2023) applied a hybrid non-local thermodynamic equilibrium (non-LTE) method, performing non-LTE spectral synthesis with DETAIL and SURFACE on LTE model atmospheres computed with ATLAS12, and derived abundances by fitting individual element lines. Markova et al. (2018) and Carneiro et al. (2019) used FASTWIND models but adopted different analysis strategies. Namely, Markova et al. (2018) derived N abundances by fitting a set of nine lines with fixed microturbulent velocities of ξ = 10 km s⁻¹ for cooler and ξ = 15 km s⁻¹ for hotter stars, whereas Carneiro et al. (2019) used the curve-of-growth. We use the results of N surface abundances of these studies to assess the consistency of abundance determinations obtained by different authors, with different model atmosphere codes, and analysis techniques. As a sanity check, we compared the observed spectra with the closest models from our grid using both our results and those from previous publications for all common stars. In general, for our grid, our results provide a better match to the diagnostic lines analyzed in this work. This confirms that our tools and methodology introduce no significant systematics.

In Fig. 5, we compare our abundances with previous determinations. Different colors and symbols identify the reference studies. We group the various atmosphere codes and abundance-determination approaches – summarized in the previous paragraph – into separate panels. We use solid symbols for stars whose ϵ_N values agree within the uncertainties and transparent symbols for inconsistent determinations. For reference, we include the cosmic abundance standard (CAS) for nitrogen from Nieva & Przybilla (2012), hereof considered the reference baseline abundance for our sample.

Fig. 5. Comparison of N abundances derived in this work with values from previous studies. Different panels present different abundance analysis methods and atmosphere codes (see text). Transparent symbols correspond to stars with discrepant comparisons (beyond uncertainties). Vertical and horizontal dotted lines indicate the CAS values at the solar position (Nieva & Przybilla 2012). The solid and dashed diagonal lines respectively represent the one-to-one relation and the mean relation with respect to the literature (obtained by fitting the consistent stars considering a fixed slope of unity), with the offset between both relations indicated at the top of the figure.

The comparison sample spans in ϵ_N values between 7.6 and 8.5 dex, which broadly overlaps with our covered range. The distribution, however, is not uniform: 21 stars cluster between 7.7 < ϵ_N < 8.2 dex, while only seven⁷ have ϵ_N ≥ 8.2 dex. Despite this uneven sampling, the number of comparable stars is sufficient across the full range to draw meaningful conclusions. We find overall good agreement with previous studies, with a small systematic shift toward lower abundances in our results. The offset amounts to 0.06 dex, below the typical abundance uncertainty (0.15 dex) and the intrinsic scatter in the comparison (0.11 dex).

In the first panel (literature abundances derived with spectral synthesis using CMFGEN), we find good agreement for stars with ϵ_N ≲ 8.1 dex. At higher abundances (ϵ_N ≳ 8.1 dex), our results show a slight systematic offset toward lower values (Δϵ_N ∼ 0.1 dex). However, the relatively large uncertainties in this regime make most of the common stars consistent within errors across the different studies.

In the second panel, we include the two studies that use the same atmospheric code as in this work. In the case of Carneiro et al. (2019) – whose methodology most closely resembles ours – we find very good agreement for stars with ϵ_N < 8.15 dex. However, two objects (HD 12 993 and HD 303 311, shown as semi-transparent gray star markers in the middle-right and lower-left regions of Fig. 5, respectively) have significantly higher abundances in our analysis. For the O6.5V((f)) Nstr star HD 12 993, our fit reproduces the observed spectrum more accurately, while for the O6V((f))z star HD 303 311 we obtain a compromise solution, with some lines fitted better and others worse. This highlights that the line selection for abundance determinations directly affects the results. We share only three stars in common with Markova et al. (2018). Two display good agreement between both studies, while for the third one they obtain a higher abundance. Given the limited statistics and the different methodologies – ours being more systematic – we do not consider this significant. Finally, in the third panel, we find very good agreement with Aschenbrenner et al. (2023), with no significant systematic differences across the full range of abundances covered.

4.2. Overall distribution

In Fig. 6, we present the relation between N and He abundances (the latter coming from Simón-Díaz et al. 2026) for the analyzed stars, using the same color and symbol coding as in Fig. 1. As a reference, we include predictions from single-star evolutionary models computed with GENEC for initial masses between 20 and 60 M_⊙ and rotational velocities v_ini/v_crit = 0.2–0.4 (gray shaded region, Ekström et al. 2012, and private communication). We also mark the CAS values for the solar neighborhood from Nieva & Przybilla (2012) in light brown. The cross in the bottom right corner indicates the standard uncertainties for both quantities. Following Simón-Díaz et al. (2026), we divide the sample into three groups based on their He abundances divided by the vertical dotted lines. These groups – defined as He-low, He-normal, and He-rich by Simón-Díaz et al. (2026) – correspond to Y_He ≤ 0.08, 0.08 < Y_He ≤ 0.12, and Y_He > 0.12 by number (i.e., Y ≤ 0.24, 0.24 < Y ≤ 0.32, and Y > 0.32 by mass), respectively. The horizontal dotted line identifies the threshold in N abundance below which the estimated abundance is clearly incompatible with the CAS baseline, even when we take into account the typical uncertainties of the determinations (Δϵ_N ∼ 0.15 dex; Sect. 4.3).

Fig. 6. Distribution of the sample in the N against He abundance plane (cf. Fig. 1 in Martínez-Sebastián et al. 2025). The gray area marks the region covered by Ekström et al. (2012) models for 20–60 M⊙ and vini/vcrit = 0.2–0.4 during the main sequence; the light brown area indicates CAS values from Nieva & Przybilla (2012). Symbols are the same as in Fig. 1; red crosses in the background mark SB1 systems. Vertical dotted lines mark the He low and rich limits (YHe ≤ 0.08 and YHe > 0.12); the horizontal dotted line separates N-low stars. In the upper X-axis, He is shown by mass; in the lower X-axis, He is shown by number.

From a theoretical standpoint, we expect a strong correlation between Y_He and ϵ_N for single stars; significant surface helium enrichment (i.e., Y_He ≥ 0.12) should occur only after it has been clearly enriched in nitrogen (i.e., ϵ_N ≳ 8.5 dex for CAS initial abundances). This strong correlation was observed in LMC stars by Rivero González et al. (2012b) and Grin et al. (2017), and also including the SMC in Martins et al. (2024). In our sample, most stars (∼65%) agree with the predictions of single-star evolution models within their uncertainties in this projection. However, several outliers remain. Namely, stars with (i) unexpectedly low He abundances (left region; discussed in Simón-Díaz et al. 2026, and later in this work); (ii) N abundances below the CAS reference (under the horizontal dotted line); and (iii) relatively low N enrichment (ϵ_N ≲ 8.4 dex) compared to their clear He enhancement (Y_He ≳ 0.12).

Thanks to the larger size of the analyzed sample compared to previous works, Martínez-Sebastián et al. (2025) were able to identify, for the first time, a non-negligible group of stars with high He but mild N enrichment (corresponding to the third group of outliers previously indicated). Notably, the sample considered in the aforementioned work differed from the one in this study in two key aspects. First, the selection criteria for targets to be considered in our previous study were less restrictive. In particular, we included stars classified as LC I and II, and some other cases which has been discarded for the present work (Sect. 2). Second, the group identified as He-low – which comprises ∼28% of the current sample – was excluded in our previous work, since their low He abundances lacked a clear physical explanation and their omission did not impact the main conclusions of that study. In the present work, we identify 11 out of the 19 stars with Y_He ≥ 0.12 with N abundances under the expected values according to Ekström et al. (2012) single-star evolutionary models even considering their uncertainties (i.e., about 60% of the He-rich subsample). Notably, this fraction is somewhat lower than the ∼80% reported in Martínez-Sebastián et al. (2025). The difference arises primarily from the reanalysis we have performed, which incorporates additional diagnostic lines and yields slightly revised abundances. While these variations remain within the expected uncertainty range, they reconcile part of the He-rich population with predictions from the considered models in this diagram. Nonetheless, the presence and frequency of outliers among the He-rich stars remain statistically significant and in line with the interpretation proposed in Martínez-Sebastián et al. (2025) that binary interaction is a likely origin for these objects.

Complementary to Figs. 6, 7 presents the observed N-abundance distributions for stars in the three regimes of He abundances (Sect. 5). From top to bottom, the panels correspond to the He-low, He-normal, and He-rich subsamples. Each panel shows the histogram of the sample, including the uncertainty of each bin, with the number of sources indicated on the left y-axis. The empirical probability density function – computed by summing the individual density functions of all measurements – is overplotted, with its scale shown on the right y-axis. The CAS value is indicated by the gray vertical band in all panels. Vertical dotted lines separate the N-low regime. Bins corresponding to ON stars are highlighted in black.

Fig. 7. Observational histograms (steps) and probability density estimates (solid lines) of ϵN for the sample, separated – from top to bottom – into He-low, He-normal, and He-rich (red, green, and blue, respectively). The CAS ϵN value from Nieva & Przybilla (2012) is indicated by the shaded vertical band. Vertical dotted lines mark the threshold separating N-low stars. The area with diagonal black lines highlights the ON stars.

The three distributions differ. An Anderson–Darling test (Anderson & Darling 1954; Scholz & Stephens 1987) between the He-low and He-normal groups, and between the He-normal and He-rich groups, yields p-values of 5 ⋅ 10⁻⁵ and 3 ⋅ 10⁻⁴, respectively. Therefore, we could reject the null hypothesis that the samples are drawn from the same parent distribution. We discuss each subsample separately in the following paragraphs.

Roughly ∼27% of the analyzed sample (32 stars) show Y_He ≤ 0.08, lower than the expected birth abundance even within uncertainties (left region in Fig. 6). Simón-Díaz et al. (2026) propose that the estimated abundances for this group of stars is likely a consequence of the impact of additional light from a nearby source contaminating the analyzed spectra (also Sect. 4.3 and Appendix C.4). This group has an SB1 incidence of ∼48% and spans ϵ_N values from ∼7.4 to ∼8.4 dex. Its N abundance distribution is asymmetric, with a global maximum at ϵ_N ∼ 7.6 dex – significantly below the CAS value for nitrogen – and a secondary peak around ϵ_N ∼ 7.95 dex. A total of 15 stars have ϵ_N below CAS, 11 of which remain inconsistent even within uncertainties (∼34% of the subset against ∼6% in the rest of the sample).

The majority of the analyzed sample (∼56%; 65 stars) show normal He abundances. This group has an SB1 incidence of ∼42% and spans ϵ_N values from ∼7.4 dex to ∼8.8 dex. The N abundance distribution is more symmetric than that of the He-low sample, with a maximum around ϵ_N ∼ 7.9 dex, coinciding with the secondary peak of the He-low group. The decline at ϵ_N values below CAS is steeper than at higher abundances, with only ∼24% with lower values. Moreover, only five stars show ϵ_N values inconsistent with CAS within uncertainties, all of them SB1. Overall, ϵ_N ≈ 8.45 dex emerges as an empirical upper limit for nitrogen enrichment without accompanying helium enrichment during the main-sequence phase for dwarfs and giants.

The remaining 19 stars (∼16% of the analyzed sample) are He-rich, with an overall SB1 incidence of ∼33%. The N abundance distribution of this latter group is bimodal, with a first peak around ϵ_N ∼ 8.1 dex and a second around ϵ_N ∼ 8.5 dex, separated by a gap at ϵ_N ∼ 8.4 dex. Among the tail of stars with lower N abundances (11 stars, none of them with an abundance below the CAS value), the SB1 fraction is ∼40%, while in the higher-ϵ_N group (eight stars) it decreases to ∼33%. Five out of the 9 stars comprising the high N abundance peak are classified as ON (a further discussion on the others is presented in Appendix D).

4.3. Low nitrogen sample

Within the current theoretical framework, we cannot explain N abundances under birth values, as nitrogen is expected to be enhanced during the main sequence. Therefore, we assess the reliability of the N abundance estimates for the 10 stars (75% SB1) with ϵ_N values below the CAS⁸. Simón-Díaz et al. (2026) found a similar result for stars with helium under CAS values (Group #2 in their work). Similar to them, we test two possible explanations: (i) an overestimated microturbulence, or (ii) contamination of the continuum by an external source.

We test the first hypothesis in Appendix C.3, where we show that this subsample cannot be explained by microturbulence effects. Alternatively, we examined the potential impact of the presence of light-pollution by a companion star on the continuum normalization of the analyzed spectrum in Appendix C.4. Our analysis shows that such dilution effect can cause spurious determinations of N abundances which might be lowered by ∼0.15 dex or even more in the case of undetected SB2 systems. This would be coherent with the high binarity incidence in this group (text and third column in Table C.3). Notably, the same effect similarly influences the derived Y_He abundances, consistent with the higher frequency of N-low stars among the He-low subsample – both effects likely share a common origin.

4.4. Distribution of stars in the Hunter diagram: General considerations

Although our study focuses primarily on the regime of slow rotators, we can still draw meaningful conclusions regarding the role of rotational mixing in the chemical enrichment of the stellar surface. To this end, we present in the upper row of Fig. 8 the so-called “Hunter diagram” (introduced by Hunter et al. 2008), displaying the relation between N abundance and v sin i for the sample. We separated the sample in the three groups according to their surface He abundance. This provides a complementary view to Figs. 6 and 7. From left to right, the panels correspond to the He-low, He-normal, and He-rich subsamples. The symbols indicate luminosity class (same code as in Fig. 1). We note that in this diagram we display the present-day v sin i, which may differ significantly from the rotational velocity at birth due not only to projection effects, but to the braking during the evolution (depending mainly on the mass loss, magnetic field and angular momentum transport). In fact, the initial rotational velocity should be higher than the present-day v_rot to enrich the surface with nitrogen in single-star evolution.

Fig. 8. Top: Surface N abundance versus projected rotational velocity. From left to right, He-low, normal, and rich sample. For reference purposes, we include evolutionary tracks for vini/vcrit = 0.2 and 0.4 from GENEC (Ekström et al. 2012, and priv. comm.; dotted and dashed lines; Mini = 20, 25, 32, and 40 M⊙) and MESA models, adopting the Mix1 prescription from Keszthelyi et al. (2022) (equivalent to Brott et al. 2011, solid and dashed dotted lines: Mini = 20, 25, 30, and 40 M⊙). Model segments with YHe ≥ 0.12 are highlighted in blue. The same markers as in Fig. 1 are used. Bottom: For the subsets in the upper panels, v sin i distributions separated in low, normal, and high nitrogen (red, green, and blue histograms, respectively).

For reference, we compare our observations with evolutionary tracks computed with the GENEC code (Ekström et al. 2012, and private communication) for main-sequence stars with initial masses between 20 and 40 M_⊙ and two initial rotational velocities (0.2 and 0.4 v_ini/v_crit; dotted and dashed lines, respectively). We also include MESA models using the Mix1 prescription from Keszthelyi et al. (2022) (equivalent to the models by Brott et al. 2011, hereafter referred as “Mix1 MESA models”) for the same initial masses, initial rotational velocities and initial chemical composition (solid and dash-dotted lines, respectively; described in Appendix A in Martínez-Sebastián et al. 2025). We select these models as they include different prescriptions, particularly of convection (Schwarzschild criterion with an overshooting edge vs. Ledoux criterion with semi-convective area) and rotational mixing (transport of angular momentum and nuclides due to meridional circulation and turbulent shear vs. diffusive formalism including magnetic fields); both works consider the mass loss prescriptions from Vink et al. (2000, 2001). Those represent two standard yet different assumptions in single-star evolutionary models, leading to a more efficient mixing and surface braking in Ekström et al. (2012). The models are color-coded by He abundance: gray for Y_He ≤ 0.12 and blue for Y_He > 0.12. Under each panel, we included the v sin i distribution separated into three N abundance subgroups (ϵ_N ≤ 7.7 dex, 7.7 < ϵ_N ≤ 8.45 dex, and ϵ_N > 8.45 dex; red, green, and blue histograms, respectively).

The two main N abundance groups identified in Fig. 7 among He-low stars (Y_He ≤ 0.08) are also evident in the left panel of Fig. 8. Approximately 30% of them correspond to low ϵ_N stars (discussed in Sect. 4.3) and have a fairly flat v sin i distribution. The remaining ∼70% of these stars show normal to moderately enriched nitrogen (ϵ_N < 8.4), and their v sin i distribution is concentrated in the range between ∼40 and 85 km s⁻¹. The distribution of stars from this latter group in the Hunter diagram is quite similar to the corresponding one in the He-normal group, but lacking stars in the low and high tails of the v sin i distribution.

The same two groups of N abundances are identified in stars with 0.08 < Y_He ≤ 0.12 (middle panels) but, as discussed in Sect. 4.2 (also Fig. 7), the percentage of stars with abnormally low N abundances is much smaller. Again, the group of stars with normal to moderately enriched N abundances (up to ϵ_N ∼ 8.5, group 2 in Markova et al. 2018) are mostly concentrated at v sin i ≲ 60 km s⁻¹ (∼80%); however, four of them are also found at larger v sin i values (≥100 km s⁻¹). Similarly to the He-low group, the v sin i distribution is fairly flat for the N-low stars (ϵ_N ≤ 7.7 dex; red histogram) in the He-normal sample.

In the He-rich regime (Y_He > 0.12; right panel), two clearly differentiated groups are observed, corresponding to the bimodal distribution peaks in Fig. 7. The very N-rich group (upper part) comprises eight stars with ϵ_N ∼ 8.6 dex, five of them classified as ON (all the stars with this qualifier in the sample). They show an approximately flat distribution in v sin i, ranging from ∼30 to ∼130 km s⁻¹. The second group consists of nine stars centered around ϵ_N ∼ 8.1 dex, but with a larger internal scatter (Δϵ_N ∼ 0.25 dex). These stars span the full range of projected rotational velocities covered in this study and show no evident correlation between ϵ_N and v sin i. Finally, two stars with very low projected rotational velocities exhibit CAS nitrogen abundances; however, their helium content (Y_He = 0.12) places them at the boundary of the He-rich regime, making them also consistent with the He-normal population.

5. Comparison with evolutionary models

5.1. Helium rich sample

Among the two sets of evolutionary models considered in Fig. 8, only GENEC models with v_ini/v_crit > 0 reproduce surface He abundances above Y_He ≥ 0.12 during the main sequence, and these cases are always accompanied by significant N enrichment (ϵ_N ≥ 8.45 dex; Fig. 6). Consequently, all stars in our He-rich subsample (right panel of Fig. 8) would be expected to display at least this level of nitrogen enrichment if they followed single-star evolutionary paths. However, several objects deviate from this prediction, showing enhanced helium but comparatively low N abundances. These inconsistent stars correspond to the lower-N peak in the He-rich distribution of Fig. 7 and to the lower region of the right-hand side of Fig. 6, which we have argued to be binary products (Martínez-Sebastián et al. 2025, and Sect. 4.2).

In the framework of single-star evolution, the eight He-rich stars that are also N-rich can only be reproduced by models assuming very high initial rotational velocities (v_ini/v_crit > 0.2). Five of these stars are classified as ON, while the others are justified to be so in Appendix D and from here on are referred as so. They rotate, on average, at v sin i ∼ 80 km s⁻¹, compared to v sin i ∼ 60 km s⁻¹ for the He-normal sample. Also, the associated v sin i distribution seems to be fairly flat. However, this can be a low number effect, highlighting again the need to increase the statistics for the He-rich group.

We highlight these eight targets with black-edged squares in Fig. 9, where the color (and marker size) represents the surface N abundance. As a reference, we include GENEC evolutionary tracks with v_ini/v_crit = 0, 0.2, and 0.4 (solid, dotted, and dashed lines, respectively), marking in blue the regions where ϵ_N ≥ 8.45. The stars under discussion appear more enriched than their neighbors, yet they do not occupy any particular region of the sHRD. Moreover, when compared with predictions from Ekström et al. (2012) models with high initial rotational velocities, five of them appear too close to the ZAMS to be He enriched through single-star evolution. This is consistent with previous findings (e.g., Martins et al. 2015b; Simón-Díaz et al. 2026), which showed that ON stars cannot be systematically reproduced by single-star evolutionary models. Rivero González et al. (2012b) found a similar N and He surface overenrichment in the bulk of O-type stars in the LMC, and interpreted it as an indication of efficient mixing during the very early phases of stellar evolution of these stars.

Fig. 9. Spectroscopic HR diagram of the analyzed sample with color and marker size scaled to N abundance. Solid black squares highlight stars with YHe > 0.12 and ϵN > 8.45 dex. As a reference, we show GENEC tracks with vini/vcrit = 0, 0.2, and 0.4 (solid, dotted, and dashed lines, respectively.) Blue segments correspond to phases with ϵN ≥ 8.45 dex.

Alternatively, Bolton & Rogers (1978) proposed a binary-interaction origin for ON stars, which subsequent observational studies have reinforced (e.g., Boyajian et al. 2005; Li & Howarth 2020). In this framework, mass transfer can reproduce the observed surface abundances (Jin et al. 2026) while changing the T_eff and log g of the star, mimicking the position of younger objects in the sHRD (McCrea 1964). This scenario also provides a straightforward explanation for the higher average rotational velocities observed among ON stars, as spin-up is an expected outcome of mass accretion (de Mink et al. 2013; Benvenuto et al. 2025). Moreover, the fact that they are far from critical rotation would point to case A mass transfer preventing the gainers spinning up to critical velocities (Langer 2012; Simón-Díaz et al. 2026).

The subsequent evolution of the binary system may result in a runaway star following the supernova explosion of the initially more massive component (Blaauw 1961; Leonard & Duncan 1988), potentially leading to a higher fraction of runaways among ON stars. Using complementary information from Carretero-Castrillo et al. (2023, 2026), and private communication, we find 50% of runaways among the stars which are simultaneously He- and N-rich, compared to only 15% in the He-normal sample. We note that this process can also account for the three stars that could, in principle, be reproduced by single-star evolution (located in the upper left corner of Fig. 9). Binary interaction therefore emerges also as the most plausible scenario for the physical origin of He-rich, N-rich stars; however, a single star origin cannot be rejected for some of them. This physical interpretation for the groups comprising all ON in our sample, but beyond the morphological classification, supports the discussion in Simón-Díaz et al. (2026).

To summarize, we have presented observational evidence challenging the single star origin of He-rich O-type stars under the current paradigm. However, these characteristics are naturally explained if they are binary interaction products. This agrees with findings in Simón-Díaz et al. (2026) regarding He abundances in the full IACOB sample of Galactic O-type stars (not limited in v sin i). Consequently, we argue that binary interaction is the dominant channel for surface helium enrichment in O-type stars. In line with this, Cazorla et al. (2017) found that approximately half of their sample of 40 fast-rotating Galactic massive stars could not be reproduced by single-star evolutionary models. In particular, they reported unexpectedly high He enrichment and the presence of ∼20% of stars with low nitrogen abundances. This behavior is consistent with the binary-interaction scenario proposed in this work. Nevertheless, this interpretation should be further tested through dedicated analyses of N abundances in larger samples of fast-rotating Galactic O-type stars, explicitly considering the potential role of binary interaction.

5.2. Normal helium sample

We argue that the He-low sample is dominated by spurious abundance determinations and hinders a meaningful comparison with single-star evolutionary models. In addition, the star comprising the He-rich sample are most plausibly binary products. Therefore, the He-normal stars (middle panel in Fig. 8) constitute a more optimal subsample for tracing single-star evolutionary trends. The removal of the two aforementioned groups from this sample produces a cleaner set of data to confront with single-star evolutionary models. Nonetheless, it can still contain binary interaction products lacking clear observational signatures (e.g., Jin et al. 2026). Moreover, the presence of binaries, even preinteraction, can play an unclear yet important role in the internal transport through tides (Sciarini et al. 2026).

Within the sample of He-normal stars, the objects with ϵ_N ≥ 7.9 dex remain unreproducible by models without efficient rotational mixing (i.e., Mix1 MESA models). Although these models can reach N abundances up to ϵ_N ∼ 8.25 dex for v_ini/v_crit = 0.4, the bulk of our observed stars exhibit substantially lower projected rotational velocities, which would require a strong surface braking mechanism than observed (Holgado et al. 2022).

Models including more efficient rotational mixing (Ekström et al. 2012) predict higher N abundances than observed (ϵ_N > 8.4 dex compared to the measured range 7.7 dex < ϵ_N ≤ 8.4 dex) for the group of stars⁹ with v sin i ≲ 50 km s⁻¹. Considering projection effects, these stars could be rotating up to ∼50% faster on average, which would bring them closer – still slightly slower – to the v_ini/v_crit = 0.2 evolutionary tracks. We expect that models computed with lower initial rotational velocities but adopting the same mixing prescriptions as Ekström et al. (2012) would better reproduce the bulk of our observations.

To further investigate surface enrichment throughout stellar evolution, we use log g as a proxy for evolutionary status in left and central panels of Fig. 10. We include evolutionary tracks with initial masses between 20 M_⊙ and 40 M_⊙ and v_ini/v_crit = 0.2 and 0.4. In the left pannel, models from GENEC (only extreme cases for v_ini/v_crit = 0.4 are included); in the middle panel; Mix1 MESA models. In both cases, we use the same coding as in Fig. 8). The right panel presents the corresponding ϵ_N histogram, separated by luminosity class using the same color scheme.

Fig. 10. Left: Nitrogen abundance as a function of surface gravity for the He-normal subsample. The dotted and dashed lines correspond to GENEC evolutionary tracks with vini/vcrit = 0.2 and 0.4, respectively (the latter only include tracks for Mini = 20 M⊙ and Mini = 40 M⊙). The blue segments indicate regions where models reach YHe ≥ 0.12. Middle: Same targets as in the left panel but using Mix1 MESA models with vini/vcrit = 0.2 and 0.4 (solid and dotted dashed lines, respectively) as reference. Right: Density histogram of N abundance grouped by luminosity class using the same color scheme as in the left panel.

Stars labeled as a LC V and IV occupy a similar gravity range (log g ∼ 4.15 − 3.5 dex) and exhibit comparable abundance distributions, both peaking around ϵ_N ∼ 7.9 dex, consistent with GENEC evolutionary tracks with v_ini/v_crit = 0.2. LC III stars occupy lower gravities (log g ∼ 3.6 − 3.3 dex, with one outlier at log g ∼ 3.8 dex) and are slightly more N enriched than LC IV and V, with their distribution peaking at ϵ_N ∼ 8.1 dex and reaching values up to ϵ_N ∼ 8.4 dex. However, their enrichment level is not as pronounced relative to dwarfs as expected from GENEC evolutionary models with the considered initial rotations. In line with the previous discussion, we expect that lower initial rotational velocities better reproduce the stars with modest enrichment and low v sin i. Nonetheless, given the initial rotation distribution reported by Holgado et al. (2022), we would expect our sample to be dominated by stars born with v_ini/v_crit ∼ 0.2. If this were the case, and the stars in our sample are the result of subsequent braking, their present-day surface abundances would show lower enrichment than predicted by these models.

In this projection, the data for LC III stars seem to be better reproduced by the Mix1 MESA models with 0.2 ≲ v_ini/v_crit ≲ 0.4 (models in the middle panel of Fig. 10). However, this is difficult to reconcile with the low projected rotational velocities observed in Fig. 8. Although we expected models similar to Ekström et al. (2012) with lower initial rotational velocities to better reproduce LC III stars in both projections, they would face difficulties reproducing at the same time LC V objects with the largest log g. Therefore, the difficulties faced by both of the analyzed set of models – two extreme cases of enrichment – to reproduce the observed distribution of ϵ_N, indicates that rotational mixing alone is not enough to account for N enrichment in this sample.

6. Conclusions

In this paper we have present the results from an N abundance analysis of a curated sample of 117 Galactic O-type stars with v sin i ≤ 150 km s⁻¹ from the IACOB project. The sample was selected to maximize the reliability of the N abundance analysis, ensuring homogeneity and quality results. This work provides a uniform and statistically significant dataset and thus allows for a robust discussion of N enrichment in slow rotating O-type stars.

From the total of 117 analyzed stars, 28 are in common with previous studies that also derived N abundances (Section 4.1), and we found an overall good agreement with them. Though our work has a small systematic offset toward lower abundances (ϵ_N = 0.06 ± 0.11 dex), it can still be considered consistent since it is within the typical uncertainty (Δϵ_N ∼ 0.15 dex). This comparison highlights the strong consistency of results when comparing large samples.

We found distinct N abundance distributions across the three helium regimes defined in Simón-Díaz et al. (2026) (Sect. 4.2). Stars with Y_He ≤ 0.08 exhibit a higher fraction of objects with abnormally low N abundances compared to the rest of the sample (Sect. 4.3). We argue that the origin of both the low N and He abundances is the same: spurious determinations caused by light pollution of the analyzed spectra by a close-by object.

On the opposite side, stars defined as He-rich (Y_He > 0.12) display systematically higher N abundances and projected rotational velocities than the rest of the sample (Section 5.1). Within this group, we identified two distinct subpopulations based on their N abundance, and neither is totally compatible with single-star evolution. The first consists of moderately N-rich stars whose properties are consistent with the outcomes of past mass-transfer events (Martínez-Sebastián et al. 2025). The second comprises a much more enriched group (ϵ_N ≥ 8.45 dex) whose locations in the sHRD are inconsistent with the high He abundance predicted by single-star evolution, even in the most mixing efficient models (i.e., Ekström et al. 2012). We interpret these stars as mass gainers mimicking younger objects after mass accretion.

Given that the two He-rich subgroups in our sample – highly or comparatively low N enriched – are better explained as binary products, we propose that binary interaction is the dominant mechanism driving surface helium enrichment in slow rotating O-type stars. The fact that the faster mean rotation in the He-rich sample is still far from critical points to a case A mass transfer that prevents critical rotation. In this framework, helium overabundance in main-sequence massive stars could serve as a relatively simple yet powerful indicator of past binary interaction. A complementary study in this line is presented in Simón-Díaz et al. (2026).

As a corollary of the previous discussions, the He-normal stars (0.08 < Y_He ≤ 0.12) emerge as the most suitable – yet contaminated – subsample for investigating single-star evolution (Section 5.2). Within this group, the N abundances of most LC IV and V are consistent with stars born as slow rotators (v_ini/v_crit ≲ 0.2) and evolving as single stars with efficient mixing. For giants, however, we did not observe the expected abundance pattern (i.e., systematically higher ϵ_N), although their distribution is shifted toward higher nitrogen values. In their case, models with less efficient mixing and higher initial rotation rates could better reproduce the observed surface abundances but would fail to match the low projected rotational velocities. Therefore, we claim that rotational mixing alone is not enough to explain the N enrichment pattern in this group.

In conclusion, this study (together with Simón-Díaz et al. 2026) provides observational evidence supporting the idea that the helium enrichment observed in main-sequence massive stars is better explained in the current paradigm as a consequence of past mass-transfer events. Consequently, we caution against the use of samples including He-rich O-type stars to calibrate single-star evolutionary models. In addition to He abundances, the analysis of N and, in the future, C and O abundances for large samples of O-type stars is essential to drawing further constraints in the internal mixing during the main sequence.

The selection criteria imposed on our sample in terms of effective temperatures and wind strengths helped minimize potential biases. However, they also make further work of including the so-far excluded stars necessary in order to investigate the late main-sequence phase of the most massive stellar objects. In the same vein, the methodology followed imposed a limit to only study slow rotators. A dedicated abundance study focused on fast rotators (which have been proposed to be binary products; de Mink et al. 2013; Britavskiy et al. 2023) could provide additional evidence to test our hypothesis.

Acknowledgments

We thank the anonymous referee for their careful review of the manuscript and the consequent comments and suggestions provided. C.M-S, G.H, and S.S-D acknowledges the support from the State Research Agency (AEI) Spanish Ministry of Science and Innovation (MICIN) and the European Regional Development Fund (FEDER) under grant PID2021-122397NB-C21 and grant PRODUCTOS DE LA INTERACCION DE ESTRELLAS MASIVAS REVELADOS POR GRANDES SONDEOS ESPECTROSCOPICOS, with references PID2024-159329NB-C21. This publication made use of the IAC HTCondor facility (http://research.cs.wisc.edu/htcondor/), partly financed by the Ministry of Economy and Competitiveness with FEDER funds, code IACA13-3E-2493. The project leading to this application has received funding from European Commission (EC) under Project OCEANS – Overcoming challenges in the evolution and nature of massive stars, HORIZON-MSCA-2023-SE-01, No G.A 101183150 Funded by the European Union. C.M.S received funding from the ERASMUS+ Program – student mobility for traineeship from the University of La Laguna (course 2024-2025). G.H. received the support from the “La Caixa” Foundation (ID 100010434) under the fellowship code LCF/BQ/PI23/11970035.

References

Appendix A: Table of resulting stellar parameters

Table A.1 summarize the parameters of the sample. The two last columns correspond to the results of this work. Spectral Classification from GOSSS, when possible.

Appendix B: Selection criteria

Figure B.1 summarizes the results of the visual preanalysis described in Sect. 2, shown in the T_eff – v sin i plane. As discussed there, we only include stars for which we have achieved a good fitting of the H and He I-II diagnostic lines when performing the IACOB-GBAT analysis (Simón-Díaz et al. 2026). Black filled symbols indicate stars in which we clearly detected at least one of the indicated diagnostic lines. Gray filled symbols mark stars with detections of at least one line of the ion, although none of them strong. Open symbols represent stars where none of the visualized lines of the ion were found (Sect. 2). Stars showing emission lines – excluded from the analysis – are marked in blue. For completeness, we also include stars with winds stronger than log Q =  − 12.7, indicated by star-shaped symbols and following the same color scheme. Only six such stars are present in the visualized sample of 145, proving that excluding them has a negligible effect on the overall sample statistics.

Fig. B.1. Summary of detection of key diagnostic lines from different N ions in the stars from the sample (Sect. 2). From left to right, and top to bottom: detections of key diagnostic lines from N II, N III, N IV, and N V, respectively. Filled black, gray, and empty gray symbols for stars displaying strong, weak and no lines, respectively. Circle and star symbols for stars with log Q ≤ 12.7 and log Q > 12.7, respectively. Dotted lines indicate the final selection criteria; area marked as "C" corresponds to the final analyzed stars.

Our analysis reveals that, with the exception of two objects, all stars in the initially considered sample displaying nitrogen lines include N III among them. Due to the very limited number of N IV and N V lines, we restricted to stars also exhibiting N III lines (T_eff ≲ 40 kK). For v sin i > 150 km s⁻¹, nitrogen lines become so broadened that they are only clearly detected in only about half of the stars, whereas for lower v sin i values the detection rate exceeds 85 %. To avoid larger biases toward higher abundances, we limited the sample to slow rotators (v sin i ≤ 150 km s⁻¹). We note that a similar argument applies to stars with v sin i > 100 km s⁻¹, as five targets with 100 < v sin i < 150 km s⁻¹ show no detectable lines. However, as they represents only ∼15% of this subsample, we decide to extend our analysis beyond v sin i = 100 km s⁻¹ in order to maximize our statistics while introducing only an acceptable bias. The final analyzed sample corresponds to region C in Fig. B.1, since regions A and B contain very few stars, and the analysis of region D would require a different approach (i.e., spectral synthesis instead of the curve-of-growth method).

Appendix C: Methodological limitations

In this appendix, we examine the limitations of our analysis and identify the dominant sources of uncertainty.

C.1. Abundance determination consistency

First, we performed a consistency test to estimate the intrinsic uncertainty of our methodology. We generated a synthetic spectrum from a FASTWIND model and convolved it with instrumental, rotational, and macroturbulence profiles. The original model parameters appear in the second column of Table C.1. We then added Gaussian noise centered at zero, with a standard deviation chosen to achieve S/N = 100.

We conduct a full blind analysis over this synthetic spectrum. To do so, we determine v sin i and v_mac with IACOB-BROAD (Simón-Díaz & Herrero 2014), the fundamental parameters with IACOB-GBAT (Simón-Díaz et al. 2011) and the abundance following the methods from Sect. 3. To estimate the effect the fundamental parameters determination itself, we repeat the abundance analysis with the original model parameters.

The results of both analyses, together with the original model parameters, are summarized in Table C.1. The second column lists the original parameters used to generate the synthetic spectrum. The third column presents the parameters recovered in the blind analysis. The fourth column shows the results obtained when we fixed the stellar parameters (except for microturbulence and N abundance) to their original model values.

We find good agreement between the two analyses and the model values. In both cases, the derived abundances are consistent with the nominal value within the associated uncertainties. The differences between the nominal abundance and the recovered ones remain well below the typical minimum uncertainty adopted in our study.

The scatter in the abundance when fixing the parameters to the nominal values arises from the uncertainty in the EW measurements. In contrast, in the blind analysis, we obtain a slightly higher abundance and a larger scatter, as expected from the propagation of uncertainties in the fundamental parameter determination. However, this effect on the final abundance remains smaller than the dispersion produced by the individual-line measurements.

We also investigated which is the dominant source of error in the abundance determination. To this end, we repeated the N abundance analysis on the same synthetic spectrum with different combinations of T_eff, log g, and ξ within typical uncertainties. In all runs, we set v sin i = 29 km s⁻¹, v_mac = 50 km s⁻¹, Y_He = 0.10, and log Q = −13.5. To evaluate the impact of EW uncertainties, we repeated each analysis using two sets of EW: (i) the measured values from the synthetic spectrum, and (ii) the EW from the original FASTWIND model (the corresponding results are shown in the fourth and fifth columns of Table C.2, respectively). The first three columns of Table C.2 list the parameters of the reference model used for the abundance determination. The uncertainty reported for each abundance corresponds to the scatter among the individual-line measurements.

Overall, we find a higher sensitivity to uncertainties in T_eff and ξ than in log g. However, the abundance variations introduced by these parameters remain smaller than the scatter produced by the individual-line measurements. This scatter therefore represents the dominant source of uncertainty in our analysis and covers those caused by the uncertainty in the fundamental parameters.

We note that the scatter increases significantly when the determined T_eff varies enough to affect the ionization balance. This situation is hardly rectifiable in observations where only one ionization stage is available, and it is mitigated when multiple stages are present, since we adjust T_eff to recover the correct ionization balance in nitrogen. In any case, this effect is expected to be smaller than the adopted minimum uncertainty of Δϵ_N = 0.15 dex. This value arises from typical uncertainties observed in stars for which several lines and ionization stages are available. Therefore, and as shown in this appendix, it supposes a conservative uncertainty.

C.2. Impact of microturbulence in the estimation of helium abundances

Despite its relevance in quantitative spectroscopy based on 1D model atmosphere calculations, the physical origin of microturbulence in massive stars remains poorly understood (Markova et al. 2025). As in previous works, we treat it as a free parameter and determine it simultaneously with the N abundance (ξ_N; e.g., Carneiro et al. 2019; Rivero González et al. 2012b). A similar approach was adopted by Simón-Díaz et al. (2026) in their analysis of helium, deriving ξ_He. However, as illustrated in Fig. C.1, the resulting ξ values can differ by up to ∼15 km s⁻¹. This difference might have a measurable impact in the estimated He abundances (e.g., McErlean 1998; Villamariz & Herrero 2000; Howarth & Smith 2001; Massey et al. 2013; Markova et al. 2020), potentially affecting the conclusions of this work. To assess this effect, we performed a consistency check by reanalyzing our sample with IACOB-GBAT (Simón-Díaz et al. 2011; Sabín-Sanjulián et al. 2014; Holgado et al. 2018), fixing the microturbulence to the values derived in this study with N lines. The spectroscopic parameters (i.e., effective temperature, surface gravity, and wind-strength parameter) remain consistent between both analysis – within the quoted uncertainties – for the entire sample.

Fig. C.1. Left: Difference in He abundance obtained using microturbulence derived from N and He abundance analysis, respectively, as a function of the difference between these microturbulence values. Right: Histogram of the He abundance differences derived with the two ξ determinations.

In Fig. C.2, we compare the results of the reanalysis (Y_He(ξ_N)) with those from the reference works (Y_He(ξ_He)). As a reference, we include the one-to-one relation (solid diagonal line) and the limits defining He-low and He-rich regimes (dotted lines). Open symbols denote stars for which both determinations agree within the uncertainties. Notably, only one object shows a significant discrepancy (black point in the lower-left corner). In both studies, however, this star is classified as He-low.

Fig. C.2. Comparison of the He abundances obtained using the microturbulence derived from N lines (this work) and from He analysis performed in (Simón-Díaz et al. 2026).

Complementarily, in Fig. C.1 we show the difference in He abundance obtained with both microturbulence values (ΔY_He  =  Y_He(ξ_N)  −  Y_He(ξ_He)) as a function of the difference in microturbulence (Δξ  =  ξ_N  −  ξ_He); same marker coding as in Fig. C.2). The right panel presents the distribution of ΔY_He for the full sample. As expected, we find a clear correlation between both quantities. However, even for |Δξ|∼ 15 km s⁻¹, the resulting differences in He abundance remain comparable to the typical uncertainty (ΔY_He ∼ 0.25). Notably, the only star showing inconsistency between both determinations corresponds to the largest difference in ξ.

As a general trend, we find that some stars located near the boundaries between helium groups may shift classification depending on the adopted ξ. However, no systematic trend emerges when using either set of microturbulence values. As a final check, we include Fig. C.3 (equivalent to Fig. 6, but using He abundances derived with ξ_N). The overall distribution remains consistent with that shown in Fig. 6. To quantify possible differences, we compare the number of stars in each helium group, as well as the fraction of SB1 systems, using both Y_He(ξ_He) and Y_He(ξ_N) (Table C.3). We find that the differences between groups remain within the expected uncertainties and do not affect the conclusions of this work.

Fig. C.3. Distribution of the sample in the N versus He abundance plane (equivalent to Fig. 6) for He abundances derived with ξN.

Given these consistency checks, and in the absence of a physical argument favoring one microturbulence determination over the other, we adopt the He abundances derived using the microturbulence values obtained from the parallel analyses (i.e., Simón-Díaz et al. 2026) due to the self consistency of the methodology.

C.3. Impact of microturbulence on the estimation of nitrogen abundances

We have shown that the impact of a microturbulence uncertainty of ∼3 km s⁻¹ is smaller than the typical abundance uncertainty (App. C.1). This effect may become relevant only when the microturbulence error is large or when it systematically biases the abundance. To test this possibilities, we repeated the analysis of the full sample by fixing ξ  =  10 km s⁻¹, a typical value for dwarfs. Figure C.4 compares the abundances derived with fixed (ξ = 10 km s⁻¹) and free microturbulence. Dashed horizontal and vertical lines indicate the CAS values (Nieva & Przybilla 2012), while the solid and dotted diagonal lines mark the one-to-one relation and the expected agreement within the typical uncertainty (Δϵ_N = 0.15 dex), respectively. Stars with Δξ  >  5 km s⁻¹ are highlighted with red edges.

Fig. C.4. Comparison between N abundances derived fixing microturbulence to ξ  =  10 km s−1 and those obtained following the procedure in Sect. 3. The solid line marks the 1-to-1 relation. The dotted lines indicate the typical uncertainty range. Stars with Δξ  >  5 km s−1 are highlighted with red edges.

We find good agreement between both sets of results, with abundance differences smaller than the typical uncertainty for most stars. The six objects showing larger differences remain consistent once we consider their individual abundance uncertainties rather than the general value.

C.4. Continuum dilution effect in He and N abundance determinations

In quantitative spectroscopy, additional light sources can significantly affect the analysis as they increase the continuum level, which results in the dilution of the spectral lines after normalization. This effect may impact part of our sample. In this work, we classify as SB1 those systems that show no evidence of an additional component in the visual inspection of He I 5875 or the O I 7772–7775 triplet in any available spectrum. This criterion is consistent with previous studies in the IACOB series (e.g., Holgado et al. 2018; de Burgos et al. 2023), but differs from the approach adopted in other recent works, which rely on spectral disentangling (e.g., Mayer et al. 2017; Mahy et al. 2022). This is because we noted two main limitations of the latter method for our purposes: (i) disentangling is not available for the full sample, and (ii) the results can depend on the initial assumptions, leading to different classifications across studies (e.g., HD 163 892, for which Barbá et al. 2026 find no companion, while Mahy et al. 2022 classify it as an SB2). Since no disentangling has been performed for this work, we adopt the classical criterion based on visual inspection. In this context, some stars are known to host companions that contaminate the continuum of the primary (e.g., HD 91 824 or the eclipsing binary HD 152 590), while others may remain undetected. Even in the absence of double-lined signatures, the companion may still contribute up to ∼10 % of the continuum – larger flux ratios would typically reveal a secondary component in the spectra. Additional contamination may also arise from unresolved nearby stars along the line of sight or from undetected SB2 systems, which could introduce even larger contributions to the continuum.

To analyze this effect, we followed a similar strategy to the one described in App. C.1. Namely, we used six FASTWIND models to generate synthetic spectra (circles in Fig. C.5), covering three regions of the sHRD (T_eff–log g combinations of 32 kK–3.5 dex, 35 kK–3.5 dex, and 40 kK–3.9 dex; red, green, and blue symbols, respectively) and two chemical compositions representative of ZAMS stars (Y_He = 0.10 and ϵ_N = 7.8 dex; filled symbols) and more evolved objects (Y_He = 0.15 and ϵ_N = 8.6 dex; open symbols).

Fig. C.5. Results of the analysis of synthetic spectra from FASTWIND models. Six parameter combinations: Teff–log g = 32 kK–3.5 dex, 35 kK–3.5 dex, and 40 kK–3.9 dex (red, green, and blue symbols, respectively), each with YHe–ϵN = 0.10–7.8 dex and 0.15–8.6 dex (filled and open symbols, respectively). For each model, we show the input values (circles) and the results obtained from the synthetic spectra with 0%, 5%, and 10% continuum dilution (squares, stars, and triangles, respectively). Top: Spectroscopic HR diagram. For reference, we include the non-rotating evolutionary tracks from Ekström et al. (2012). Dotted lines indicate constant gravities, starting at log g = 4.2 dex on the left and decreasing in steps of 0.2 dex. Bottom: Nitrogen abundance versus He abundance (as in Fig. 6). For reference, we include GENEC tracks for initial masses between 20 and 60 M⊙ and vini/vcrit  =  0.2 − 0.4.

For each case, we produced three synthetic spectra: one without dilution and two diluted to 95% and 90% of the original (square, star, and triangle symbols, respectively). We then analyzed these spectra —including the determination of the fundamental parameters— to quantify the impact of continuum dilution on the derived abundances.

The results are presented in Fig. C.5. The top panel shows the sHRD with the derived fundamental parameters. For reference, we include the non-rotating evolutionary tracks from Ekström et al. (2012) (solid gray lines) and lines of constant surface gravity (gray dotted lines), starting at log g = 4.2 dex on the left and decreasing toward the right in steps of 0.2 dex. The bottom panel displays N versus He abundances, together with GENEC tracks for initial masses between 20 and 60 M_⊙ and v_ini/v_crit  =  0.2 − 0.4.

We do not find any significant effect of continuum dilution on the determination of T_eff from the He I–He II ionization balance. In contrast, we detect an effect on log g, which we associate with distortions in the wings of the Balmer lines. These distortions can mimic higher surface gravities. This effect is not uniform across the sHRD, becoming more pronounced at higher gravities.

For the models with T_eff = 40 kK, log g = 3.9, Y_He = 0.10, and ϵ_N = 7.8 dex, we could not determine ϵ_N due to the weakness of the N lines, which prevents a reliable analysis. For the remaining models, we represented the uncertainty in ϵ_N as the scatter among the different analyzed lines, without imposing the minimum uncertainty of Δϵ_N = 0.15 dex —adopted elsewhere in this work. In these cases, the continuum dilution leads to lower derived ϵ_N and Y_He values, mimicking the effect of lower surface abundances.

As discussed above, this effect is expected to be more pronounced in undetected SB2 systems, where the contribution of the companion to the continuum is larger. To test this effect, we analyzed HD 93 028 and HD 101 413, two known SB2 systems composed of an O-type primary and an early B-type secondary, where the latter rotates significantly faster, making its detection more difficult (Putkuri et al. 2026). For these stars, we derive Y_He = 9.0 ± 3.4 and ϵ_N = 7.60 ± 0.15, and Y_He = 7.3 ± 1.5 and ϵ_N = 7.51 ± 0.15, respectively. While these results are consistent with the expected dilution effect, the larger continuum contribution from the rapidly rotating companion leads to lower derived He and N abundances than in the synthetic model tests, where the dilution was limited to 10 % (maximum).

It is important to note that none of the studied changes shift the derived abundances toward the region of He-rich and relatively low N. This rules out continuum dilution as a plausible explanation for the group of stars discussed in Martínez-Sebastián et al. (2025).

Appendix D: Revisiting the ON qualifier

A particularly relevant stellar group in the context of our study are the ON stars. This morphological classification correspond to O-type stars with enhanced N absorption, and C and O deficient compared to their spectral type and luminosity class (Walborn 1976, 2003). In particular, ON stars are delineated by the relative strength of N III 4634, N III 4640, and C III 4650 (Sota et al. 2011). Their enhanced nitrogen lines may indicate N over-enrichment at the surface (Martins et al. 2015b). Table D.1 summarizes the relevant information of the five objects of this type identified in our sample. Interestingly, all of them are He-enriched, consistent with the results of Martins et al. (2015b) and Simón-Díaz et al. (2026). They are also systematically N-enriched (ϵ_N ≥ 8.45 dex; most enriched bins in Fig. 7); however, not all N-enriched objects are classified as ON. Since this classification is purely morphological, CNO anomalies should be interpreted in the context of the corresponding spectral type.

As a visual approach to place N-rich stars in an evolutionary context, Fig. 9 shows the sample of stars in the sHRD, with color and symbol size indicating ϵ_N. As reference, we include GENEC evolutionary tracks for initial rotational velocities v_ini/v_crit = 0, 0.2, and 0.4 (solid, dotted, and dashed lines, respectively), distinguishing between models with ϵ_N < 8.45 dex (gray) and ϵ_N ≥ 8.45 dex (blue). Stars with Y_He ≥ 0.12 and ϵ_N ≥ 8.45 dex are highlighted with black squares.

All these stars display comparable abundances and are more enriched than other objects with similar spectral type, but not all of them are classified as ON. To investigate this, we revisited the morphological classification of those stars by directly comparing their spectra with those of reference ON stars of similar spectral type and luminosity class (Table D.2) and revisiting the relative strength of N III 4634, N III 4640 and C III 4650. Based on these comparisons, we propose stars in the firs column of Table D.2 to be reclassified as ON.

Our proposed additions are consistent with a correlation between line strengths and nitrogen enrichment, reinforcing the physical meaning of this morphological classification. Moreover, this exercise illustrates the diagnostic potential of quantitative analyses in spectral classification, as EW remain unaffected by extrinsic line broadening.

One particular case deserves additional discussion, as it is already classified as Nstr but we propose it should be reclassified as ON. For HD 12 993, we derive a N abundance comparable to that of HD 90 273, whereas a lower ϵ_N would be expected given its spectral type. Moreover, when compared to other O6.5 V stars in the sample, HD 12 993 shows a markedly higher nitrogen enrichment relative to its classification, further supporting its reclassification as ON. Consequently, CPD -58 2611 should also be included in the ON group. According to our revised classification, ON stars form a clearly distinct group within our sample, characterized by ϵ_N > 8.45 dex.

All Tables

All Figures

Fig. 1. Distribution in a spectroscopic HR diagram of the sample of Galactic O-type stars considered in this work (large filled symbols). For reference purposes, the full analyzed IACOB sample of O-type stars is included as small size circles, while red and yellow symbols are for LC I and II. The GENEC evolutionary tracks during the main sequence are plotted for vini/vcrit = 0 with solid lines. Dotted lines are used for gravities decreasing from log g = 4.2 dex (left) in steps of Δ log g = 0.2 dex.

In the text

Fig. 2. Example of the mapping of N abundance as a function of microturbulence. Colored lines show the N abundance that minimizes Ml as a function of microturbulence for each of the analyzed lines in HD 206 183. The horizontal dashed line marks the ξ value determined with the N III multiplet. The colored X marker corresponds to the abundance of each line at this microturbulence. The dotted line indicates the average abundance, with the gray area representing the 1σ dispersion around this value.

In the text

Fig. 3. Example of the diagnostic plot used to correct possible systematic effects in the abundance determination. It shows abundances of the different analyzed lines in HD 206 183 for a fixed microturbulence velocity. Green and red symbols correspond to N II and N III lines, respectively. Solid lines in the same colors indicate the average abundance for each ion, with the shaded areas being their 1σ confidence intervals.

In the text

Fig. 4. Example of the evaluation of the N abundance fit in selected spectral regions. Orange vertical lines denote the nitrogen transitions employed in the abundance determination; dotted gray vertical lines indicate other N lines present in the model but excluded from the analysis. We note that at the temperature of this object, N III 4642 is dominated by an overlapping O II line (Rivero González et al. 2011).

In the text

Fig. 5. Comparison of N abundances derived in this work with values from previous studies. Different panels present different abundance analysis methods and atmosphere codes (see text). Transparent symbols correspond to stars with discrepant comparisons (beyond uncertainties). Vertical and horizontal dotted lines indicate the CAS values at the solar position (Nieva & Przybilla 2012). The solid and dashed diagonal lines respectively represent the one-to-one relation and the mean relation with respect to the literature (obtained by fitting the consistent stars considering a fixed slope of unity), with the offset between both relations indicated at the top of the figure.

In the text

Fig. 6. Distribution of the sample in the N against He abundance plane (cf. Fig. 1 in Martínez-Sebastián et al. 2025). The gray area marks the region covered by Ekström et al. (2012) models for 20–60 M⊙ and vini/vcrit = 0.2–0.4 during the main sequence; the light brown area indicates CAS values from Nieva & Przybilla (2012). Symbols are the same as in Fig. 1; red crosses in the background mark SB1 systems. Vertical dotted lines mark the He low and rich limits (YHe ≤ 0.08 and YHe > 0.12); the horizontal dotted line separates N-low stars. In the upper X-axis, He is shown by mass; in the lower X-axis, He is shown by number.

In the text

Fig. 7. Observational histograms (steps) and probability density estimates (solid lines) of ϵN for the sample, separated – from top to bottom – into He-low, He-normal, and He-rich (red, green, and blue, respectively). The CAS ϵN value from Nieva & Przybilla (2012) is indicated by the shaded vertical band. Vertical dotted lines mark the threshold separating N-low stars. The area with diagonal black lines highlights the ON stars.

In the text

Fig. 8. Top: Surface N abundance versus projected rotational velocity. From left to right, He-low, normal, and rich sample. For reference purposes, we include evolutionary tracks for vini/vcrit = 0.2 and 0.4 from GENEC (Ekström et al. 2012, and priv. comm.; dotted and dashed lines; Mini = 20, 25, 32, and 40 M⊙) and MESA models, adopting the Mix1 prescription from Keszthelyi et al. (2022) (equivalent to Brott et al. 2011, solid and dashed dotted lines: Mini = 20, 25, 30, and 40 M⊙). Model segments with YHe ≥ 0.12 are highlighted in blue. The same markers as in Fig. 1 are used. Bottom: For the subsets in the upper panels, v sin i distributions separated in low, normal, and high nitrogen (red, green, and blue histograms, respectively).

In the text

	Fig. 9. Spectroscopic HR diagram of the analyzed sample with color and marker size scaled to N abundance. Solid black squares highlight stars with YHe > 0.12 and ϵN > 8.45 dex. As a reference, we show GENEC tracks with vini/vcrit = 0, 0.2, and 0.4 (solid, dotted, and dashed lines, respectively.) Blue segments correspond to phases with ϵN ≥ 8.45 dex.
In the text

Fig. 10. Left: Nitrogen abundance as a function of surface gravity for the He-normal subsample. The dotted and dashed lines correspond to GENEC evolutionary tracks with vini/vcrit = 0.2 and 0.4, respectively (the latter only include tracks for Mini = 20 M⊙ and Mini = 40 M⊙). The blue segments indicate regions where models reach YHe ≥ 0.12. Middle: Same targets as in the left panel but using Mix1 MESA models with vini/vcrit = 0.2 and 0.4 (solid and dotted dashed lines, respectively) as reference. Right: Density histogram of N abundance grouped by luminosity class using the same color scheme as in the left panel.

In the text

Fig. B.1. Summary of detection of key diagnostic lines from different N ions in the stars from the sample (Sect. 2). From left to right, and top to bottom: detections of key diagnostic lines from N II, N III, N IV, and N V, respectively. Filled black, gray, and empty gray symbols for stars displaying strong, weak and no lines, respectively. Circle and star symbols for stars with log Q ≤ 12.7 and log Q > 12.7, respectively. Dotted lines indicate the final selection criteria; area marked as "C" corresponds to the final analyzed stars.

In the text

	Fig. C.1. Left: Difference in He abundance obtained using microturbulence derived from N and He abundance analysis, respectively, as a function of the difference between these microturbulence values. Right: Histogram of the He abundance differences derived with the two ξ determinations.
In the text

	Fig. C.2. Comparison of the He abundances obtained using the microturbulence derived from N lines (this work) and from He analysis performed in (Simón-Díaz et al. 2026).
In the text

	Fig. C.3. Distribution of the sample in the N versus He abundance plane (equivalent to Fig. 6) for He abundances derived with ξN.
In the text

	Fig. C.4. Comparison between N abundances derived fixing microturbulence to ξ  =  10 km s−1 and those obtained following the procedure in Sect. 3. The solid line marks the 1-to-1 relation. The dotted lines indicate the typical uncertainty range. Stars with Δξ  >  5 km s−1 are highlighted with red edges.
In the text

Fig. C.5. Results of the analysis of synthetic spectra from FASTWIND models. Six parameter combinations: Teff–log g = 32 kK–3.5 dex, 35 kK–3.5 dex, and 40 kK–3.9 dex (red, green, and blue symbols, respectively), each with YHe–ϵN = 0.10–7.8 dex and 0.15–8.6 dex (filled and open symbols, respectively). For each model, we show the input values (circles) and the results obtained from the synthetic spectra with 0%, 5%, and 10% continuum dilution (squares, stars, and triangles, respectively). Top: Spectroscopic HR diagram. For reference, we include the non-rotating evolutionary tracks from Ekström et al. (2012). Dotted lines indicate constant gravities, starting at log g = 4.2 dex on the left and decreasing in steps of 0.2 dex. Bottom: Nitrogen abundance versus He abundance (as in Fig. 6). For reference, we include GENEC tracks for initial masses between 20 and 60 M⊙ and vini/vcrit  =  0.2 − 0.4.

In the text