© ESO, 2014

1. Introduction

The chemical evolution of galaxies is dominated by the evolution of early-type stars, since these objects are the progenitors of core-collapse supernovae, and therefore contribute to stellar nucleosynthesis in a pronounced way. In this context, important issues are the effects of rotation, especially for that of rotational mixing, on the evolution of massive stars (e.g., Heger & Langer 2000; Meynet & Maeder 2000), as well as spatial and temporal variations of the chemical composition within the Galactic disk (e.g., Fuhrmann 2008; Przybilla 2008). Quantitative spectroscopic analyses of B- and late O-type stars allow for atmospheric parameters and chemical elemental surface abundances to be inferred with high precision, which directly addresses both of the aforementioned topics (see Nieva & Przybilla 2012). Due to the high frequency of binary systems among early-type stars (see, e.g., Sana et al. 2012; Chini et al. 2012), analysis techniques, which are also able to deal with spectra of double-lined spectroscopic binary systems (SB2) are desirable.

Quantitative spectroscopy is based on the comparison of synthetic and observed spectra. Owing to the multi-dimensionality of the parameter space involved (which include the following for B- and late O-type stars: effective temperature T_eff, surface gravity log (g (cm s^-2)), microturbulence ξ, macroturbulence ζ, projected rotational velocity ν sin (i), radial velocity ν_rad, metallicity Z, and elemental abundances {n(x)}), investigations are time-consuming since an iterative approach is required (for details see Nieva & Przybilla 2010). Starting from initial estimates for the entire set of parameters, individual variables are refined by using spectral indicators that are sensitive to as few parameters as possible to reduce the complexity of the problem. In early-type stars of solar metallicity, for instance, Stark-broadened hydrogen and helium lines are primarily affected by changes in T_eff, log (g), and n(He) while they are comparatively insensitive to all others. Consequently, these features allow for the temperature, surface gravity, and helium abundance to be constrained. The use of multiple ionization equilibria, which requires that spectral lines of different ionization stages of the same element indicate equal abundances, yields further constraints on T_eff and log (g) but also on ξ and n(x). Matching the strength of spectral lines and their shape allows for ζ and ν sin (i) to be derived. Because of the highly non-trivial coupling of different parameters, the adjustment of individual variables involves re-evaluation of most of the available indicators, leading to an iterative procedure.

Although those iterative steps can be automated to speed up the investigation (see Lefever et al. 2010), the underlying strategy is still prone to miss the global best solution. On the one hand, not all parameters are varied at the same time but many of them separately so that correlations between them are neglected (see Mokiem et al. 2005). On the other hand, parameters are constrained from selected spectral indicators or windows instead of exploiting the information encoded in the entire spectrum. A global analysis method, which is a method simultaneously probing all parameters while considering the maximum useful spectral range, is therefore our goal.

Automated fitting techniques are suitable for this purpose. Moreover, automation is far less subjective, since the matching of theory and observation is based on a mathematical measure, such as a χ² criterion instead of visual inspection. This is particularly important when one wants to analyze larger samples in a homogeneous manner. The size of the corresponding multi-parameter space, however, requires calculations of numerous synthetic spectra, which is computationally expensive and, therefore, a major obstacle for automated fitting. To minimize the number of calculations involved, synthetic atmospheres may be computed on demand in the course of the fitting process, as realized by Mokiem et al. (2005). In this way, spectra are computed only if they are actually used. Nevertheless, even very efficient fitting algorithms can take from several dozens to hundreds of iterations to find the best solution, which implies a non-negligible run-time of the fitting process. This drawback can be overcome by making use of pre-calculated model grids in which interpolation between grid points can be used to evaluate the fitting function. Unfortunately, sufficient sampling of the whole multi-parameter space is typically not possible, given its large dimension. Consequently, grid-based fitting methods are usually restricted to small subspaces by either keeping some parameters fixed when computing tailored grids (see Castro et al. 2012) or limiting the allowed parameter range, thus reducing the advantages of global automated fitting.

However, due to some unique properties of spectra of early-type stars, such as the low density of spectral lines and the continuous opacity that is dominated by hydrogen and helium, many parameters – in particular, the elemental abundances of the trace elements – are independent of each other. Exploiting this fact, it is possible to probe the entire parameter space by computing only a tiny fraction of it. Based on this idea, we have developed a grid-based global fitting method that facilitates quick and precise determinations of the atmospheric parameters of B- and late O-type stars, which takes non-local thermodynamic equilibrium (non-LTE) effects into consideration (Sect. 2). Furthermore, it is shown that the accuracy of the analysis is generally not limited by statistics, such as the signal-to-noise ratio (S/N) of the observed spectrum, but rather by systematics, such as the uncertainties in atomic data (Sect. 3). For demonstration purposes, the method is then applied to three well-studied early-type stars in Orion and to three SB2 systems yielding atmospheric and fundamental stellar parameters (Sect. 4). A discussion of the results obtained (Sect. 5) is rounded off by a summary (Sect. 6).

2. Setting up the fitting function

2.1. Input physics: models and codes

Synthetic spectra are calculated by following the hybrid, non-LTE approach discussed by Przybilla et al. (2006) and Nieva & Przybilla (2006, 2007, 2008). The structure of the atmosphere, such as the stratification of temperature and density, is based on line-blanketed, plane-parallel, homogeneous, and hydrostatic LTE-model atmospheres calculated with Atlas12 (Kurucz 1996). Deviations from local thermodynamic equilibrium are then accounted for by applying updated versions of Detail and Surface (Giddings 1981; Butler & Giddings 1985). The Detail code numerically solves the coupled radiative transfer and statistical equilibrium equations to obtain population numbers in non-LTE, which in turn are input to the Surface code to compute the final synthetic spectrum that uses more detailed line-broadening data. Comprehensive, state-of-the-art model atoms (see Table 1) allow spectral lines to be modeled with high fidelity for chemical elements accessible in the optical spectra of most B- and late O-type stars. Although partly based on LTE concepts, the hybrid method outlined above is consistent with full non-LTE calculations for B-type stars (see Nieva & Przybilla 2007; Przybilla et al. 2011) but is considerably faster and, most importantly, is able to handle more realistic representations of model atoms. For the first time here, we consistently use the concept of LTE opacity sampling to treat background line opacities (employing the realization of Kurucz 1996) in all computational steps, instead of relying on pre-calculated opacity distribution functions.

2.2. Probing the parameter space

Stellar atmospheres are described by a variety of physical quantities. The implementation of those parameters in synthetic model computations varies depending on their nature. Within the scope of this paper, it is instructive to distinguish between primary atmospheric parameters, which are quantities whose effects can be modeled solely by solving the coupled radiative and structural equations describing an atmosphere, and secondary ones, which are quantities whose effects can be incorporated after an atmosphere has been built. For OB-type stars near the main sequence and their evolved progeny, BA-type supergiants, primary parameters are T_eff, log (g), ξ, Z, and {n(x)}, while ζ, ν sin(i), and ν_rad are secondary ones. By definition, the effective parameter space for which model atmospheres have to be calculated is therefore spanned by the primary parameters. Given a grid-based fitting method that interpolates a regular multi-dimensional mesh of synthetic spectra, the dimension of this parameter space is the product of the individual dimensions of the primary parameters. Due to this multiplicity, each chemical element included in the analysis drastically increases the number of required models. For instance, calculating a grid with ten elements at nine abundances each implies a factor of 9¹⁰, which has to be multiplied with the dimension given by the remaining primary parameters. Since it is infeasible to compute such a large number of models in a reasonable time, a simplification is necessary. To this end, one can exploit the unique properties of (nearby) early-type stars: due to their young ages, the metallicity of these stars is very well approximated by the cosmic standard (Nieva & Przybilla 2012), in general, so that Z = 0.014 can be assumed throughout. In massive stars, the latter implies that all metals can be considered as trace elements, which implies that moderate changes of their abundances do not significantly affect the atmospheric structure. This occurs because their respective bound-free opacities are small compared to hydrogen and helium and because line blanketing is relatively insensitive to abundance changes of a specific metal. Altering the abundance of an individual metal has therefore little to no effect on spectral lines of the other metals, except in the case of line blends. However, the density of lines is low in optical spectra of early-type stars, so that intrinsic blends are actually rare. (With some more restrictions, this even applies to the ultraviolet region.) Moreover, many of these line blends are caused by macroscopic flux redistributions, which can be, for example, via macroturbulence, stellar rotation, or the finite resolving power of the spectrograph, and are thus not an issue, since they do not affect the microscopic physics governing radiative transfer. Consequently, instead of calculating 9¹⁰ combinations of abundances, it is sufficient to compute models for nine pre-chosen abundance values for each of the ten species. All of these models contain only lines of one specific trace element, while standard abundances are assumed for the other metals to account for their background opacities. Those 10 × 9 spectra serve as a basis from which each of the 9¹⁰ combinations can be constructed, such as via multiplication of the normalized base spectra. Note that helium cannot be treated as a trace element, since corresponding abundance changes have serious impacts on the atmospheric structure, for instance, via the mean molecular weight or its contribution to the continuum opacity. Hence, the number of required synthetic models is the product of the dimensions of T_eff, log (g), ξ, n(He), and the base dimensions of the trace elements.

2.3. Numerical implementation of the model

From the considerations of the previous subsection, the fitting function is constructed as follows: for each trace element of interest, a multi-dimensional regular mesh spanned by effective temperature (12 000 to 34 000 K, step size: 1000 K), surface gravity (range depending on temperature but typically within 3.0 to 4.6 dex, 0.2 dex), microturbulence (range depending on temperature but typically within 0 to 20 km s^-1, 2 km s^-1), helium abundance (solar value minus 0.3 dex to solar value plus 0.5 dex, 0.2 dex), and abundance of the species itself (typically solar value minus 0.8 dex to solar value plus 0.8 dex, 0.2 dex) has been calculated in advance. This grid can easily be extended whenever needed. Spectra for arbitrary parameters within the mesh are approximated by linear interpolation. At this point, synthetic models contain only spectral lines of hydrogen, helium, and the element under consideration. To combine the normalized spectra of the individual trace elements, each of them is divided by a model with the same specifications but which contains only hydrogen and helium lines. Multiplying the resulting “corrected spectra” with each other and with the hydrogen-helium model leads to a model that takes all of the primary parameters into account. This simplification is well justified as long as spectral lines of different metals are not blended due to microscopic flux redistributions, such as microturbulence, natural line width, and thermal- or pressure broadening. Furthermore, even for blends of weak spectral lines, the method is a good approximation as interaction effects are tiny in that case. Secondary atmospheric parameters are incorporated afterwards: macroturbulence and stellar rotation via convolution of the synthetic spectrum with a joint profile function, which is obtained from numerical stellar disk integration (Gray 2005) with a linear limb-darkening law (using coefficients by Claret & Bloemen 2011), and radial velocity by shifting the entire wavelength scale, according to the Doppler formula.

2.4. Extension to composite spectra

Multiplicity is an important issue for massive stars, as the majority of them probably resides in systems of two (or more) components (see Sana et al. 2012; Chini et al. 2012). While visual binaries or systems with a much fainter secondary can be treated as single stars in the spectroscopic analysis, this is clearly not the case for double-lined spectroscopic binary systems where features of both stars are visible in the spectrum. Stars classified as chemically peculiar objects based on low-resolution and/or low-S/N spectra might actually be unrecognized SB2 systems where the depths of the absorption lines have been altered due to continuum emission of a companion.

The fitting method presented here is also capable of dealing with composite spectra in a simplified manner. To do so, models for the normalized single-star spectra are created separately for the primary (f_p) and secondary (f_s) component following the procedure given in Sect. 2.3. To obtain the model for their normalized composite spectrum (f_comp), these individual contributions have to be summed up while taking their weights into account. The latter are given by the components’ continuum fluxes (f_p,cont, f_s,cont) and by the ratio of their projected, effective¹ surface areas, which is parametrized by one additional free parameter A_eff,s/A_eff,p: (1)Note that this approach conflates all binary interactions, such as mutual distortions and their respective effects on the emitted spectra into the parameter A_eff,s/A_eff,p. Consequently, the geometric interpretation of A_eff,s/A_eff,p as the ratio of the effective surface areas of secondary to primary component is only valid as long as those interactions are negligible, such as for well detached systems. Otherwise, A_eff,s/A_eff,p is a combination of the ratio of effective surface areas and a fudge factor that accounts for all missing interaction effects. Of course, the ability of A_eff,s/A_eff,p to compensate these shortcomings is limited and certainly fails for very close or contact systems where a much more sophisticated method is necessary (see, for instance, Palate & Rauw 2012; Palate et al. 2013). Nevertheless, the approach outlined here offers a very fast and efficient way to derive atmospheric parameters and elemental abundances of SB2 systems from the analysis of a single epoch composite spectrum.

2.5. Comparison with observation

To compare models with observations, the synthetic spectra are finally convolved with the instrumental profile, which is a Gaussian curve whose full width at half maximum Δλ = λ/R is a function of the spectral resolving power R.

The goodness of fit of a model is then derived from its absolute χ²: (2)Here, f_i, f_model,i and δ_i are the observed flux, model flux and uncertainty in the observed flux at data point i, respectively. The sum is taken over all pixels i in the spectrum, which either excludes only those lines that have well-known shortcomings, are missing in our models, or are of non-photospheric origin, such as telluric lines and diffuse interstellar bands. The lower the value of χ², the better the quality of the fit². Other ways of measuring the fitness exist but may have different properties. One reason to choose a standard i-based over a line-based χ² criterion (as, for instance, done by Simón-Díaz et al. 2011b; Castro et al. 2012) is that it gives lines with more data points a larger weight than those with less information. Moreover, the corresponding χ² distribution is well studied, which allows statistical uncertainties of the input parameters to be deduced from the χ² statistics (see Sect. 3.3).

The analysis is carried out completely within the Interactive Spectral Interpretation System (ISIS, Houck & Denicola 2000), which simplifies the investigation tremendously. Apart from many useful functions and tools, various minimization algorithms, such as the simplex or gradient methods are available. As shown in the next section, the χ² landscape of our problem is generally very well behaved so that the absolute minimum is found after a relatively small number of steps, which allows for quick and efficient analyses.

3. Formal tests and discussion of uncertainties

3.1. Noise estimation

Calculating χ² requires proper knowledge of the measurement uncertainties δ_i, which either are systematic in nature and, for instance, caused by an incorrect continuum normalization, or, more importantly, have statistical fluctuations, such as noise n_i. The latter can be estimated from an observed spectrum, which does not need to be flux calibrated, in an easy, fast, and robust way when assuming that the noise n_i of data point i obeys a Gaussian probability distribution p(n_i) with a mean value of zero and a priori unknown standard deviation σ_i: (3)For small regions (several data points), the measured flux f_i as a function of the wavelength λ can be approximately written as the sum of a linear function a + bλ_i, which represents the first two terms in a Taylor expansion of the pure signal and a noise component n_i: (4)To estimate the noise level σ_i, consider the quantity Δ_i defined as (5)Here, w_{i − 2} and w_i+2 are weight factors and chosen such that only the noise terms in Eq. (5) remain³: (6)The reason for comparing f_i to the (weighted) average of data points i − 2 and i + 2 instead of i − 1 and i + 1 is that adjacent pixels are likely correlated. For example, this is due to detector cross-talk or actions taken during data reduction like the wavelength calibration. Assuming that there is no correlation with the next neighbor but one implies that Eq. (3) is valid for points i, i − 2, and i + 2, and the probability distribution p(Δ_i) reads: (7)Here, δ is the dirac delta function. For adjacent data points, it is well justified to assume σ_i = σ_i−2 = σ_i+2 ≡ σ, so that p(Δ_i) simplifies to (8)Consequently, the distribution of Δ_i = f_i − (w_i+2f_i−2 + w_i+2f_i+2) is a Gaussian with a standard deviation defined by Eq. (8). Extending the assumption of a constant noise level σ_i = σ to a statistically significant number of data points allows σ, which is the statistical component of the uncertainty δ_i in Eq. (2), to be derived from the measurable distribution of Δ_i. If the reduced χ² at the best fit is larger than 1, it might also be necessary to consider a systematic component of δ_i (see footnote ⁴).

3.2. Performance and reliability of the method

Before fitting real spectra, several formal tests were carried out to examine the properties of the automated method. With this aim, mock spectra were constructed from synthetic ones by adding Gaussian-distributed noise that corresponded to different S/N. A spectral range [3940 Å, 7000 Å] was chosen to match the minimum wavelength coverage of standard high-resolution spectrographs. Regions in that interval that are generally affected by telluric features were excluded. The spectral resolving power was set to R = 45 000, which is very close to the resolution of the Echelle spectra which are analyzed in this study.

Table 2 lists the results of this exercise for ten exemplary cases. The input parameters and, thus, the global minimum were recovered with excellent accuracy after only a run-time of few minutes on a standard 3.1 GHz single-core processor and independent of the choice of the starting parameters within the grid, which shows that our method is fast and reliable.

It is important to stress here that all mock spectra were constructed from complete Surface models, which, in contrast to the fitting function, treat every microscopic line blend correctly by simultaneously computing lines of all chemical species under consideration. Additionally, models off the grid points were chosen to check that our mesh is sufficiently spaced for the linear interpolation scheme applied. Because the differences of input and output values in Table 2 are often covered by the very small statistical uncertainties (see Sect. 3.3) that result from the high S/N assigned to the mock spectra, we conclude that inaccuracies introduced by simplifications in our approach are negligible.

In a second step, three mock composite spectra were created with the help of Eq. (1). The parameters chosen here are motivated by real SB2 systems and anticipate the results presented in Sect. 4. They cover a sharp-lined, well-separated and, thus, easy to analyze system and a very difficult configuration with heavily blended spectral features. Similar to the previous tests, most of the input parameters are recovered with very high precision or at least within the derived uncertainties, as seen in Table 3. In particular, the degree of accuracy in the inferred parameters of both components of the extremely blended composite spectrum is astonishing, hence, making us quite confident that our method is also highly suitable for investigations of SB2 systems.

Although our method is able to model individual abundances if necessary, we generally prefer to assume an identical chemical composition for the two components within the binary system. In this way, the number of free parameters and, consequently, the numerical complexity of the problem is significantly reduced. In cases where the parameters of the secondary component are only poorly constrained due to their little impact on the composite spectrum, it is even necessary to impose these constraints, which compensate for the lack of spectral indicators to derive reasonable atmospheric parameters. Note that the assumption of an equal chemical composition is well justified for SB2 systems containing B-type or late O-type stars. On the one hand, the components of SB2 systems are in general similar regarding to their masses (otherwise, the flux contribution of the fainter companion would not be visible in the spectrum), ages (the whole system formed at once), and pristine chemical composition (both components formed from the same building material). On the other hand, processes causing chemical anomalies are rare among B- and late O-type stars and primarily affect helium (Smith 1996). Since chemical peculiarities are possibly even less frequent in detached binary systems (Pavlovski & Southworth 2013), elemental abundances are expected to evolve in the same way in both components. To estimate the influence of this approximation on the spectral analysis, Table 3 lists the results obtained from fitting the three mock composite spectra with adjustable abundances and an equal chemical composition. Even for the system with different individual abundances, the results derived by assuming an identical chemical composition are very satisfying. In particular, this is with respect to the primary component, which dominates the spectrum and in this way also the estimates for the system abundances and their respective confidence limits. As a consequence, the actual abundances of the secondary component may sometimes lie outside of the uncertainty intervals determined for the binary system as a whole. Nevertheless, it is obvious that the decision of whether or not to use separate abundances during the fitting process depends on the individual object and has to be checked, for example, a posteriori by inspecting the final match of the model to the observation.

3.3. Discussion of statistical and systematic uncertainties

The accuracy of spectral analyses is generally limited by the quality of the obtained data and the ability of the model to reproduce the observation. As shown in this subsection, shortcomings in the model, which may be due to inaccurate atomic data or deficient line broadening theory, are the main obstacles to overcome to perform more precise investigations.

Statistical uncertainties result from the noise in the observed spectrum and can be deduced from the χ² statistics in the standard way: starting from the best fit with a reduced χ² of about one⁴, the parameter under consideration is increased/decreased, while all remaining parameters are fitted, until a certain increment Δχ² from the minimum χ² is reached (for details, see Bevington & Robinson 1992). Here, each Δχ² corresponds to a confidence level; for example, Δχ² = 6.63 is equivalent to the 99%-confidence interval (see the magenta line in Fig. 1 for an illustration). The resulting uncertainties are, of course, only trustworthy if the δ_i of Eq. (2) are reasonably estimated. The method outlined in Sect. 3.1 can do so, as shown by the tests with mock spectra with known noise level (see Tables 2 and 3). Moreover, those tests, which use the same models as the fitting routine and, thus, exclude all sources of systematic errors apart from microscopic line blends, give an estimate of the statistical uncertainties that can be expected in real data with a similar S/N.

Systematic uncertainties are much harder to cope with. Sources of systematic errors occur almost everywhere in the course of the analysis (see the discussion in Nieva & Przybilla 2010). At the same time, their effects are by no means trivial, and it is extremely difficult and sometimes even impossible to quantify them. In particular, this is true for atomic data (such as energy levels, oscillator strengths, and photo-ionization cross sections), which affect individual spectral lines and the atmospheric structure. Monte Carlo simulations in the style of Sigut (1996, 1999) offer the possibility of estimating the effects on spectra caused by variations in these input data. However, a thorough error analysis has to take all sources of systematic errors into consideration at the same time to account for correlations as well, which is an unfeasible task.

Our analysis strategy is designed to keep systematic uncertainties as small as possible. For instance, an inaccurate local continuum definition can introduce considerable uncertainties to the determination of metal abundances, especially in fast-rotating stars or low-resolution spectra where metal line blends lower the actual continuum. In our routine, these effects are allowed by re-normalizing the observed spectrum with the help of the synthetic ones. Here, the latter are used to properly locate the continuum regions, which are sufficiently frequent in optical spectra of early-type stars. For these, a correction factor is obtained by dividing the (smoothed) observed data with the model data. Interpolating this factor to the whole wavelength grid gives the local continuum correction term for all spectral lines. For this approach to work, a high degree of completeness in terms of modeled lines is necessary, which is verified by high-resolution, high S/N spectra of slow rotators, as seen in Figs. 2a–i.

Another crucial part of our strategy is that we are simultaneously fitting the maximum useful range of the optical spectrum. In this way, parameters are determined not just from one or two spectral indicators but from all available ones. As the systematic errors of the individual indicators are typically independent of each other, which can be exemplified by ionization equilibria of different metals or oscillator strengths of various multiplets, there is a good chance that their effects on the parameter determination average out because some lines systematically give higher and others lower abundances, thus reducing the impact of systematics.

To crudely estimate the systematic uncertainties, we start from the assumption that they mainly appear as inaccuracies in the determination of effective temperature and surface gravity. From our extensive experience with the applied synthetic spectra, we find it realistic but conservative to assign errors of ±2% in T_eff and ±0.1 dex in log (g). The ranges given by these errors are formally treated as 99%-confidence intervals. The precision in fixing the microturbulence and the abundances of the chemical elements is then estimated from propagating the errors in T_eff and log (g). Here, a fit of all remaining parameters is performed for each pixel (that is for each combination of temperature and surface gravity) surrounded by the solid black line in Fig. 1. In the case of a binary system, this procedure is carried out for the secondary component as well. The resulting (combined) ranges of parameter values are then taken to be 99%-confidence intervals. This approach is valid as long as uncertainties induced by variations in T_eff and log (g) dominate other sources of systematic errors. While this is likely to be true for ξ and n(x), it is clearly not the case for ν_rad, ν sin(i), and ζ. Determination of radial velocities is generally limited by the accuracy of the wavelength calibration and ranges between 0.1–2 km s^-1 for common spectrographs. Projected rotational velocity and macroturbulence are incorporated via convolution with corresponding profile functions. Because of simplifications (for example in the treatment of limb-darkening or the assumption of radial-tangential macroturbulence) during the derivation of the latter (see Gray 2005), their validity may be limited to a few km s^-1.

The comparison of statistical and systematic uncertainties, as listed in Table 2, shows that our method’s total uncertainty is dominated by systematic effects down to at least a S/N of 100. However, this by no means implies that high S/N data are an unnecessary luxury. They are indispensable to detect weak features, such as contributions from a faint companion star, which would otherwise be hidden by noise. Moreover, shortcomings in the models are much more likely to remain unrecognized in low S/N spectra. This is particularly true if they can be partly compensated by tuning some fitting parameters, which, in turn, would cause erroneous results. Instead, this comparison shows that the accuracy of the presented spectral analysis technique is currently limited by modeling and not by observation.

4. Analysis

4.1. The program stars

To illustrate the capabilities of our new method, we re-analyzed three well-studied early-type stars from the Orion region and performed a spectral analysis of three SB2 systems based on very high quality single epoch spectra. For one of the binaries, two additional spectra covering different orbital phases were also investigated. Table 4 lists the relevant information about the program stars. We focused on targets which cover a wide range of effective temperatures, have moderate projected stellar rotations, and show a variety of binary configurations.

4.2. Atmospheric parameters and abundances

Atmospheric parameters and abundances are determined by fitting synthetic to observed spectra, as outlined in detail in Sect. 2. Comparisons of final, best-fitting models with observations are shown for a large portion of the used spectral range in Figs. 2a–i for the Orion stars and in Figs. 3a–i for the SB2 systems. The overall match of metal lines is almost perfect for the cooler stars of the sample and still very good for the hotter ones where our model atoms begin to be partially incomplete because they were not optimized for this temperature regime. This is particularly true with respect to O ii lines. Consequently, more regions have to be excluded from the analysis for higher effective temperatures due to (blends with) missing spectral lines.

Note that we generally exclude several He i lines from the analysis owing to recurrent issues with their detailed spectral line shapes, which is apparent only in high-resolution spectra with high S/N. This specifically affects the diffuse He i lines that show small but perceptible systematic deficiencies in their forbidden components, which can be attributed to shortcomings in their line broadening theory. These issues are independent of T_eff and log (g) and cannot be resolved, even if these lines are fitted individually. With a sufficient number of alternative, highly trustworthy He i lines present in the optical spectral range, we generally refrain from fitting the diffuse lines to avoid a possible source of systematic error and use them instead as a consistency check. Nevertheless, given the fact that the synthetic profiles of these lines match the observed ones quite well, except for the forbidden components (see Figs. 2a–i), this decision is probably too restrictive. The diffuse He i lines could therefore be considered for spectral fitting as well with resulting changes in parameters that are well below the stated systematic uncertainties.

Furthermore, we ignore the temperature-sensitive cores of the Balmer lines during the fitting process by excluding those parts of these lines where their normalized flux is smaller than a cutoff, which is typically chosen to be 0.8. The reason for this is that they are formed in the outer stellar atmosphere, where deviations from the LTE stratification are more pronounced (see Nieva & Przybilla 2007; Przybilla et al. 2011) and where the assumption of hydrostatic equilibrium also becomes less and less valid in the accelerating (weak) stellar wind. Moreover, by simultaneously fitting the entire useful spectral range, there are enough other indicators for T_eff such as (multiple) ionization equilibria available so that it is sufficient to use the Balmer line cores as a consistency check (Figs. 2a–i and 3a–i show how well this works.). As a side product, we also reduce the otherwise overwhelming influence of these lines on the parameter determination.

Table 5 lists the atmospheric parameters and abundances of the program stars. Instead of the “classical” notation for the abundance, log (x/ H) + 12, we have chosen n(x), which is the fractional particle number of species x with respect to all elements. The motivation for this is that the helium abundance is variable from star to star, which, in turn, causes the hydrogen abundance to vary since hydrogen and helium abundances are coupled via the fixed number of total particles. As a consequence, the quantity log (x/ H) can change even if x stays constant. For a better comparison of metal abundances in stars with different helium content, we therefore prefer the notation that gives abundances relative to all elements.

4.3. Stellar parameters and mass fractions

The stellar parameters mass M, age τ, and luminosity L are derived from T_eff, log (g), and ν sin (i) by fitting single-star evolutionary tracks that account for stellar rotation (M ≤ 15 M_⊙: Georgy et al. 2013; otherwise: Ekström et al. 2012). The unknown inclination term sin (i) is replaced by its spherically averaged value π/ 4 when matching ν sin (i) to the equatorial velocity ν predicted by the evolution tracks. Using the gravitational constant G, the stellar radius R_⋆ follows then from the definition of the surface gravity . This information combined with the object’s photometry can be used to derive the spectroscopic distance d to the star (Irrgang et al. 2014). Mass fractions of hydrogen (X), helium (Y), and metals (Z) are directly computed from the deduced atmospheric abundances.

Error propagation for stellar parameters and mass fractions is analogous to the estimation of systematic uncertainties in Sect. 3.3. That is, for each pixel surrounded by the solid black line in Fig. 1, they are derived as outlined in this subsection. Minimum and maximum values of the resulting distributions are again interpreted as to define 99%-confidence intervals.

The program stars’ positions in the (T_eff,log (g)) diagram are compared to evolutionary tracks in Fig. 4 and the resulting stellar parameters and mass fractions are given in Table 6. Note that the usage of single-star tracks to deduce stellar parameters of binary stars is, of course, justified only if the two binary components have not yet interacted and is otherwise an approximation.

5. Discussion

5.1. Single B- and late O-type stars in Orion

Focusing on a wide range of effective temperatures, we have selected three slow rotators (HD 35299, HD 35912, and HD 37042) from the sample of Nieva & Simón-Díaz (2011) to check our method against previous studies.

As shown in Table 5, our atmospheric parameters have excellent agreement with those derived by Nieva & Simón-Díaz. Similarly, the results for the abundances of helium, carbon, and nitrogen are perfectly consistent with each other within the error bars, even though helium was kept fixed at the solar value in the study of Nieva & Simón-Díaz. The same applies to oxygen and silicon abundances by Simón-Díaz (2010). On the other hand, there are systematic discrepancies apparent for neon, magnesium, and iron that can be attributed either to differences in the synthetic models or in the analysis strategy. For instance, several Mg ii lines, such as λ4481 Å have shown to be very sensitive to the replacement of pre-calculated opacity distribution functions, as used by Nieva & Simón-Díaz, with the more flexible concept of opacity sampling that is coherently used here throughout all computational steps, which explains the deviations in magnesium. The disagreements in neon and iron presumably arise from the underlying analysis techniques and in particular from how the microturbulence parameter is constrained. Nevertheless, it is extremely satisfying to see that the results of the two approaches match so well despite being based on contrary conceptional designs.

5.2. Spectroscopic binaries

As a first application to SB2 systems, we have analyzed the composite spectra of three binary systems. While the lines of the two components are sharp and very well separated in our spectrum of HD 119109, the opposite is true for HD 213420 (see Figs. 3a–i). In the case of HD 75821, we have further derived parameters from spectra taken at three distinct orbital phases to investigate its influence on the results.

HD 119109 (#5): to our knowledge, there is no hint for binarity in the literature for this system so far. Nevertheless, our spectrum shows that this is doubtlessly a SB2 system owing to the many lines that appear twice in the spectrum (see Figs. 3a–i). Given the high quality of our observation and the opportune orbital phase, parameters of both components can be reliably deduced. The system turns out to be composed of two relatively unevolved (see Fig. 4), coeval (, ) late-type B-stars of masses M_p = 5.2 ± 0.2 M_⊙ and when using single-star evolutionary tracks. The corresponding squared ratio of radii, , is consistent with the surface ratio , as is the spectroscopic distance, , with the parallax, . The chemical composition resembles that of the single stars studied in Sect. 5.1.

Using published radial velocity measurements, Tetzlaff et al. (2011) have proposed that this object is a runaway star with high probability based on its peculiar space motion. This conclusion should be considered as uncertain as long as the actual system velocity of this binary is unknown.

HD 213420 (#6): this well-known binary system with a period of about 880 days and a radial velocity semi-amplitude of 9 km s^-1 (Pourbaix et al. 2004) is a clear SB2 system, given the broad absorption features superimposed to He i λ4438 Å, λ6678 Å, C ii λ4267 Å, Mg ii λ4481 Å, Si iii λ4553 Å, λ4568 Å, and S ii λ5454 Å (see Figs. 3a–i). Although the signatures of the secondary component are weak and thus only detectable in the case of a high S/N, they are apparently sufficient to determine reasonable atmospheric parameters for the companion because the resulting stellar parameters paint a consistent physical picture: In addition to the finding that the ages of both components (with masses , ) are in perfect agreement (, ), the spectroscopic parameter lies within the uncertainty interval of the squared ratio of the evolutionary-track radii, . The spectroscopic distance of the system, , finally fits to its parallax, . Apart from a slight tendency to a lower metallicity (see Table 6), the chemical composition agrees with the reference stars of Sect. 5.1.

HD 75821 (#4): this eclipsing binary has a period of about 26.3 days and a radial velocity semi-amplitude of 92 km s^-1 (Mayer et al. 1997).

The spectrum best suited for the spectral analysis is the second one (b) in Table 5, since the spectral line separation is largest in this case, which reveals several pure and unblended features of the companion (see Figs. 3a–i). Reliable atmospheric and stellar parameters for both components are, hence, determinable whereby the latter assume that single-star evolutionary tracks are appropriate. Starting from this premise, the system consists of two coeval components (, τ_s ≤ 10 Myr): a massive primary (), which is slightly evolved, and a less massive secondary (), which is almost unevolved (see Fig. 4). The spectroscopic distance lies well within the 99%-uncertainty range of the parallax, . Finally, the spectroscopically deduced effective surface ratio agrees well with the squared ratio of the evolutionary-track radii, , and is further consistent with the photometric light curve (Mayer et al. 2014). The elemental abundances of the system are in line with the single stars except for the relatively low helium, nitrogen, and oxygen content (see Table 5).

The heavily blended and, hence, almost vanishing imprints of the companion on the first (a) and third (c) spectrum are insufficient to properly constrain the secondary component’s atmospheric parameters. Instead, unphysical values and large systematic uncertainties, which are induced by variations of T_eff and log (g) of the primary, are derived for the secondary’s T_eff and log (g). These error margins are, on the one hand, a direct consequence of strong correlations among certain parameters and, on the other hand, related to the fact that contributions of the secondary component barely affect the spectrum at the corresponding orbital phases. In a simplified picture, increasing the primary’s T_eff and decreasing its log (g) at the same time causes the He ii lines to become considerably deeper than actually observed, while the He i lines still fit nicely. To compensate for this, A_eff,s/A_eff,p and, hence, the influence of the secondary component, has to be significantly increased to fill the He ii lines with the continuum which thus weakens them again. However, this makes some spectral lines of the secondary component substantially too strong, which, in turn, is corrected for by smearing them out via a larger ν sin (i) or ζ that finally leads to a more uncertain determination of ν_rad, given the high degree of line blending at these particular orbital phases.

However, the primary’s properties and the surface ratio A_eff,s/ A_eff,p are nicely recovered in all three orbital phases, which gives us confidence that the presented method is generally able to determine them from one single spectrum.

6. Summary

In this paper, a novel objective method to analyze single or composite spectra of early-type stars is presented. It is based on fitting synthetic spectra to observation by using the standard concept of χ² minimization, which requires the wavelength-dependent noise of the spectrum to be known. Therefore, a simple but precise way of estimating the local noise has been developed (see Sect. 3.1). To facilitate fast and efficient analyses, we make use of pre-calculated grids of synthetic spectra, instead of computing them on demand during the fitting procedure. To sample the entire multi-dimensional parameter space at once, we exploit the unique spectral properties of early-type stars, such as the low density of lines, which reduces the number of models required by several orders of magnitude. In this way, a simultaneous fit of all parameters is possible which has the great advantage that cumbersome iterations by hand or the risk of missing the global best solution are avoided. Moreover, parameters are not only constrained from a subset of available lines but from all useful features in the spectrum. The extension to composite spectra of double-lined binary systems proves extremely valuable in the future, given the high frequency of SB2 systems among early-type stars (see Sana et al. 2012; Chini et al. 2012). In contrast to spectral disentangling techniques like those of Simon & Sturm (1994) or Hadrava (1995), the method presented here allows for – at least – parameters of the primary and the components’ effective surface ratio to be inferred from single-epoch spectra alone.

Statistical and systematic uncertainties of our method are discussed (see Sect. 3.3). The former are based here on a clearly defined mathematical measure, namely the χ² statistics, while the latter on experience. We show that systematic effects generally dominate in the high-quality regime of our observations. The analysis of a larger sample of stars thus enables us to identify possible shortcomings in our models and to derive results with significantly reduced statistical scatter.

As a case study, we have determined parameters of three well-known stars in the Orion region that turn out to be in excellent agreement with previous studies. Additionally, three binary systems have been analyzed with all of them yielding very conclusive results. Consequently, we are now in a position to homogeneously analyze large samples of early-type stars in relatively short times. The results of a comprehensive investigation of 63 nearby mid B-type to late O-type stars will be published in a forthcoming paper.

---

Acknowledgments

A.I. acknowledges support from a research scholarship by the Elite Network of Bavaria and subsequent support by the German Research Foundation (DFG) through grant He1356/45-2. M.-F.N. acknowledges an FFL stipend from the University of Erlangen-Nuremberg. We are very grateful to Thomas Dauser, Ingo Kreykenbohm, and Fritz-Walter Schwarm for maintaining and improving the Remeis computer cluster which was used for this study, to John E. Davis for the development of the slxfig module used to prepare the figures in this paper, to Jörn Wilms for proposing the χ² strategy in combination with ISIS, and to Markus Firnstein for providing joint profiles of macroturbulence and stellar rotation obtained from numerical disk integration. We thank Horst Drechsel and Pavel Mayer for drawing our attention to HD 75821 and for discussions related to it. Mahalo nui loa to Rolf-Peter Kudritzki and the University of Hawai’i for their great hospitality during A.I.’s research stay at the Institute for Astronomy and to the German Academic Exchange Service (DAAD) as well as to the Physics Advanced study program for financial support of this stay on which parts of this paper are based on. We are grateful to Artemio Herrero for a competent and constructive referee report.

References

Online material

All Tables

All Figures