1. No category

Chapter
4
Spectral analysis of the GAUDI B
star sample using a grid based fitting
method
Lefever, K.; Puls, J.; Morel, Th.; Briquet, M.; Decin. L; Aerts, C., 2007, AnalyseBstar:
an automated tool for the high-accuracy estimation of fundamental parameters for B
type stars - Application to the B-type stars in GAUDI (in preparation, to be submitted
to A&A)
In this chapter, we present the first automatic tool for the spectroscopic analysis of B-type
stars, covering the complete B-type spectral range: AnalyseBstar. We discuss our choice for
a grid method, and outline the basic properties of the algorithm. After thorough testing, we
apply the method to a sample of B-type stars, observed in preparation of the CoRoT space
mission. The wide spread in spectral types and the accuracy achieved for the fundamental
parameters, allows us to verify the currently widely used photometric calibrations in Teff and
log g.
94
GAUDI B star sample
4.1 CoRoT and GAUDI
4.1.1 The space mission CoRoT
December 27, 2006, 14:28 GMT: literally an explosive minute in the lives of every European
asteroseismologist and planetary scientist. The French-led European space mission CoRoT
(Baglin et al. 2002) has finally (and successfully!) been launched on top of a Soyuz 2.1b
carrier rocket from the ground station in Baikonour. After launch, CoRoT (COnvection,
ROtation, and planetary Transits) was placed on a circular, polar orbit, at an altitude of about
900 kilometers above the Earth atmosphere, allowing for continuous observations of two large
and opposite regions in the sky. The viewing zone is in the equatorial direction, perpendicular
to the orbital plane, thus avoiding any perturbation by the Earth straylight. The first target
field is towards Orion (the anticenter field, at right ascension 6h 50m , in the winter). As
the Earth rotates around the Sun, the telescope will repoint from time to time, in order to
keep pointing to the same field. After 150 days, when the Sun starts to interfere with the
observations, the spacecraft will turn 180 degrees towards the centre of the Milky Way (the
centre field, at right ascension 18h 50m , in the summer). Within each of these two large
regions, many different fields are selected to be monitored. CoRoT has a nominal lifetime of
2.5 years. During this period it will perform uninterrupted high-cadence (32s) high-precision
photometry (on µmag level!) per field (2.8 x 2.8◦ ) during 5 months. This never achieved
precision raises high expectations for the quality of the data and for a breakthrough in our
understanding of the stellar interiors. Indeed, one of the main objectives of the mission
is to perform asteroseismology of more than one hundred dwarfs with the goal to extract
more information from their interior. A few of the brightest among them are monitored
intensively for up to five months, whereas, simultaneously, up to nine fainter dwarfs per field
are observed. Additionally, CoRoT will provide high-accuracy photometry of thousands of
fainter objects, down to V = 16 with a cadence of 512s. In this way one hopes to achieve
information for a wide diversity of stars, covering the HRD as completely as possible.
But CoRoT has more in store. A second main objective is the search for exoplanets. Using
the transit-method, it will be attempted to detect extrasolar planets down to the size of our
own Earth.
Space-based photometry has this advantage over ground-based photometry, that it can monitor almost continuously, in this way avoiding to a very large extent any one-day alias problems. Moreover one is no longer disturbed by absorption from the Earth’s atmosphere. In
addition to these advantages, also the achieved signal-to-noise ratios are several orders of a
magnitude better than can be obtained using ground-based facilities.
4.1.2 The birth of GAUDI
The observational setup of CoRoT (observations of a small number of bright stars over a very
long period) makes target selection a crucial issue. To fully profit from the CoRoT data, it
95
4.1 CoRoT and GAUDI
Figure 4.1: Artist’s view of the satellite mission CoRoT. Courtesy: CNES
http://sci.esa.int/science-e/www/object/index.cfm?fobjectid=29381
is of major importance to assess the properties of the potential targets and thus to gather as
much a priori information as possible on them. Double-lined spectroscopic binaries (SB2s)
should be avoided, because disentangling the oscillation information of both components
from their composite frequency spectra is a very complicated process. Also, the identification of all kinds of peculiarities in both photometry and spectroscopy is necessary. Besides
this, it is important to dispose of an accurate estimate for the physical parameters of the star
(rotation, effective temperature, surface gravity, etc.), in addition to the seismic information
that will be obtained, in order to fully understand their evolution (as discussed in Chapter 1).
Due to the acute shortage of available information on the potential CoRoT targets, the need
for additional data was high. Therefore, an ambitious ground-based observing program for
more than 1500 objects was set up, under the leadership of C. Catala at Meudon, to obtain
Str¨omgren photometry as well as high-resolution spectroscopy. These data were collected in
an extended catalogue, maintained at LAEFF (Laboratorio de Astrof´ısica Espacial y F´ısica
Fundamental) and baptised GAUDI: Ground-based Asteroseismology Uniform Database Interface (Solano et al. 2005). Meanwhile, various research groups from different countries
have gathered forces to take up the huge work of analysing and interpreting all the data. Various methods have been used, each of them valid within a certain parameter range, or for a
certain type of objects. These parameters (if necessary provided with additional comments)
have been incorporated in the GAUDI database. A hierarchy between the various methods
has been adopted, to provide users with the best value available for each object. The information gathered in this database is finally fed into CoRoTSKY, the dedicated tool to prepare
and schedule the observations for CoRoT1 .
1 http://smsc.cnes.fr/COROT/A
corotsky.htm
96
GAUDI B star sample
4.1.3 GAUDI Data
Str¨omgren photometry Each summer and winter, from 2000 until 2004, two weeks of
observation time at the OSN (Observatory Sierra Nevada) were completely devoted to the
CoRoT ground-based photometry programme.
In these runs Str¨omgren-Crawford uvbyβ and Ca II H&K photometry were obtained for all
potential targets. These data have been used to derive estimates of their effective temperatures, surface gravities and metalicities (Amado et al. 2004).
Echelle spectroscopy Most spectroscopic observations were obtained by using high-resolution echelle spectrographs on different sites: FEROS (La Silla, R ∼ 48,000), ELODIE
(France, R ∼ 42,000) and SARG (La Palma, R ∼ 46,000), complemented with some additional spectra obtained by CORALIE (La Silla, R ∼ 50,000), GIRAFFE (South-Africa,
R ∼ 39,000) and the coud´e spectrograph at the 2m telescope in Tautenburg (Germany, R ∼
35,000).
4.2 Responsibilities of the Institute of Astronomy in Leuven
As a partner in the CoRoT team, ‘we’ were in charge of the detailed spectroscopic analysis
of the massive O and B-type stars. Since for these stars (except for the Be stars), no other adequate methods were available, our method (detailed spectroscopic analysis) has been ranked
first on the hierarchy list of confidence2.
The very first step to be considered, before really getting started, is to make a (rough) preselection of the targets:
– Be stars were a priori excluded. Be stars are main-sequence or slightly evolved stars, i.e.
non-supergiants. Usually they are fast rotators, surrounded by an equatorially concentrated envelope (disk) fed by discrete mass-loss events caused by yet unknown mechanisms. In the HRD, early Be stars are located at the lower border of the instability domain
of the β Cephei stars, while mid and late Be stars reside among the SPB stars. The presence of a disk causes the Balmer lines (especially Hα) to be in emission (generally, but
not necessarily, double-peaked, depending on the inclination angle). Since also the line
profiles of other atoms and the continuum level can be severely affected, they need a special approach. The FASTWIND code is not developed to treat these stars properly and
therefore we exclude them from our sample. They are the responsibility of the CoRoT Be
star working group led by A.-M. Hubert and C. Neiner: see Neiner et al. (2005), Fr´emat
et al. (2006) for the results of their investigations.
2 http://corot.oamp.fr/book/VI.3.pdf
4.2 Responsibilities of the Institute of Astronomy in Leuven
97
Figure 4.2: Breakdown by spectral type of the 190 B-type stars with spectra available in the
GAUDI database.
– Also the double-lined spectroscopic binaries (SB2) were a priori excluded from our sample, since their combined spectra make an accurate fundamental parameter estimation impossible as long as their flux ratios (and therefore their individual masses) are not known.
Also these stars need a special approach, which is the responsibility of the Binaries Thematic Team led by C. Maceroni & I. Ribas.
– The spectra with low SNR (below ±100) have to be separated from the ones with high
SNR, since bad quality data can never lead to good quality (i.e. accurate) fundamental
parameters.
– Fast rotators should be identified. Fast rotation smears out the spectral line features and
causes severe blending. Diagnostic lines become strongly affected, which makes it impossible to correctly estimate the fundamental parameters. This of course depends on
how fast the star is rotating (mainly problematic if v sin i is above 200 km/s) and on the
spectral type of the star, since the number and the strength of diagnostic lines depend on
it.
After the exclusion of the Be stars and the spectroscopic binaries, we were left with 190 Btype and 11 O-type stars. In Fig. 4.2 we show the subdivision of the B-type stars following
their spectral type (which is taken from the astronomical database SIMBAD). It is immediately clear that the number of B8 and B9 type stars is very large compared to the number of
stars at earlier spectral types. Since spectral classification relies on spectral characteristics,
these 54 B8 stars, or these 77 B9 stars, are all expected to have very similar spectra (apart
from different rotational or turbulent broadening), which means that also the resulting physical parameters are expected to be similar. Instead of analysing all stars separately, we deem it
a better approach to select a few representative cases with the best quality data to perform the
detailed spectroscopic analysis. Besides saving a lot of time, this has two more advantages:
1) it keeps our standard of high-accuracy analyses high and 2) these targets will be ideally
suited to check the validity of the currently widely used fundamental parameter calibrations.
Considering the extended sample of the GAUDI B stars to be analysed, it is no longer realistic
98
GAUDI B star sample
to perform all analyses by hand. Since no automatic tool for the spectral analysis of B-type
stars has been designed so far, we considered the time ripe to do so. In the near future, such
tool will especially prove very useful once more large-scale observing programs are being set
up, or when more and more observations of B-type stars are being gathered, and a fast and
accurate acquisition of their stellar parameters is required. Note that it is possible to apply this
procedure also for other objects (besides B stars), once suitable alterations are carried out. In
what follows we will discuss the steps considered and undertaken towards an automatic tool.
4.3 AnalyseBstar: an automatic tool for the
spectroscopic analysis of B-type stars
We present in this section the first design of an automatic method for the spectroscopic analysis of B stars with winds, which we have baptised “AnalyseBstar”. It is based on the extensive
grid of NLTE atmosphere models presented in Chapter 3 and is developed to treat large samples of stars in a homogeneous way. In contrast to a fit-by-eye method (as used for the sample
of pulsating B supergiants in Chapter 2) an automated method is not only homogeneous, but
also objective and more robust. We could not apply this to the supergiant sample of Chapter 2,
because we had insufficient spectral lines per target.
Note that for objects without winds, an automatic tool to determine simultaneously the effective temperature, gravity, microturbulent velocity and silicon abundance, has been developed
by S. Becker, using the line formation codes DETAIL, SURFACE and ANALYS of Giddings
(1981). The method is based on the isocontours of EW (and their intersection) for different
silicon and Balmer lines in the (Teff , log g)-plane (Becker & Butler 1990; Becker 1991).
4.3.1 Spectroscopic line profile fitting: methods
The detailed spectroscopic analysis of B-type stars has always been a very time consuming
job for two main reasons. The first reason is the large parameter space covered by the photospheric and wind parameters that has to be explored. The second reason is the fact that there
exists no automatic method to derive these parameters in an objective way. Fitting is still
done by using a thorough fit-by-eye method, which requires the calculation of a few tens to a
few hundreds of models per star. Fit-by-eye is a time consuming job, since new models can
only be calculated after inspection of the previous generation of models.
To cope with the large dataset of ground-based observations of CoRoT, we investigate the
possibility of automated fitting. Spectral line fitting is a clear example of an optimisation
problem: the goal is to minimise the difference (or to maximise the correspondence) between
a given observed stellar spectrum and a theoretically predicted spectrum emerging from a stellar atmosphere model. To do so, one needs to search through the parameter space spanned by
the free parameters of the stellar atmosphere model, to find the best fit (i.e. global optimum).
There are several options for the way to proceed.
4.3 AnalyseBstar
99
Fit-by-eye method
This method is the most widely spread method for performing spectral line analyses. Starting
from knowledge about spectral type and previously published values, if any, one makes a
rough guess of the initial parameters by comparing accordingly chosen grid models with the
observations. Exploiting knowledge about the specific dependency of the theoretical line
profiles on one or more of the basic parameters, one tries to obtain an optimum fit in an
iterative way through the calculation of additional models that better match the observational
data. In order to clearly see the influence of the change of a certain parameter on the line
profiles, the best strategy is to change only one or very few parameters at the same time.
However, this strategy implies a few disadvantages, such as 1) it should be avoided to solve
a multidimensional problem in a onedimensional way. A real global optimum can only be
found by allowing all parameters to vary in parallel; 2) the amount of models to be calculated
for one star can become very large. Though FASTWIND models may be calculated very
fast, this way of working remains very time consuming. Scientists unfortunately do not have
eternal time to work on one object which implies that a) analysing large numbers of objects
is out of the question, and b) the investigated part of parameter space is necessarily restricted.
A good initial guess and good knowledge of how each parameter influences the strength, the
width or the shape of a line are thus indispensable. The latter, however, can only be built
up through experience. Altogether, the risk is high that the solution that appears optimal in
the eyes of the person who performs the fit in the above way, is maybe not the global but
only a local optimum. As suggested by its name, the method severely relies on individual
eyes, which induces that the method risks of being subjective. Since there exists no standard
procedure to follow when performing a fit-by-eye method, anyone else disposing over the
same spectrum, might follow a different sequence of changing parameters and might as well
find different resulting parameters. To overcome these drawbacks, it is necessary to automate
the line profile fitting procedure. This is possible in at least two ways, which we will described
in the next two sections.
Genetic algorithm
As a first option for an automated fit procedure, one may think of using a genetic algorithm
to search for an optimal solution (e.g., Mokiem et al. 2005, see below). A genetic algorithm
(GA) is only one of the different search techniques developed within the field of artificial
intelligence (AI). The algorithm is inspired by evolutionary biology (genetics) to which it
thanks its name, and from which it inherited its terminology and techniques, such as mutation, selection, inheritance, and recombination. It is generally used to find solutions to an
optimisation and search problem. In our specific problem to find an optimum fit to an observed spectrum resulting from a set of atmosphere models, the stellar and wind parameters
play the role of genomes, which, combined, uniquely define the characterising chromosome
of each atmosphere model (so-called individual). Evolution towards a better solution usually
starts from a population of randomly generated individuals, in our case a rough grid of atmosphere models, and happens in generations. Using the DNA information of each of the
models in the current generation, new generations are created through recombination and/or
100
GAUDI B star sample
mutation, where the goodness-of-fit acts as the selection mechanism, following the principle
‘survival of the fittest’: the better a model fits, the higher the chance for it to be chosen as
a propagator. The new population is then used in the next iteration of the algorithm. Commonly, the algorithm terminates when a satisfactory fitness level has been achieved.
The main improvement compared to the fit-by-eye procedure, is the fact that this method does
no longer require human intervention, thus excluding any human bias and saving a lot of time
in the calculation of ‘bad’ models. Indeed, the selection mechanism makes sure that every
new generation consists of “fitter” solutions than the parent generation. Moreover, cleverly
written genetic algorithms are able to jump out of local optima and to find a global optimum
in a multidimensional parameter space. On the other hand the genetic algorithm approach
also has a few disadvantages. A large population of models is required for each generation,
which usually takes a huge amount of processing time. Thousands of models need to be
calculated for one single object. The information obtained through the analysis of one object,
is usually lost when treating another object, i.e. for each new object the evolution starts from
an initial guess. Parameters are drawn more or less randomly, without taking into account any
knowledge of how the parameters influence the line profiles, whereas using this experience
could considerably fasten up the fitting process.
A grid-based method
A third option for finding the best fit between theoretical and observational line profiles is
grid-based and is different from the previous automated method in this respect that the method
makes use of an existing grid without calculating additional models. This requires the grid to
be comprehensive, as dense as possible and representative for the kind of objects one wants
to analyse. The ‘subjective’ eye has to be replaced by an ‘objective’ intelligent computer programme. Through the fit-by-eye analysis of more and more different objects, one gains more
and more experience in the line profile fitting procedure. One gradually develops a ‘feeling’
for the reaction of the model and of the predicted line profiles on a change in fundamental
parameters, and slowly but surely one can even predict how large a difference is needed to
improve the fit. This human experience, as well as knowledge of the studied object need
to be translated into computer code. Similar to the fit-by-eye method, this algorithm will
follow an iterative scheme, but now in a reproducible way. Fitness is no longer judged by
visual inspection, but is fixed by a goodness-of-fit parameter. Starting from a first guess for
the fundamental parameters, based on spectral type and/or published information, improved
solutions are derived by simply comparing line profiles resulting from well-chosen existing
gridmodels to the observed line profiles. The algorithm terminates once the fit quality cannot
be improved anymore by modifying the model parameters.
4.3 AnalyseBstar
101
Justification of our choice for the grid method
The first approach, the fit-by-eye approach, is, for our purposes, out of question. The method
is a long-winded process of trial and error, and is not the appropriate choice for handling large
samples of stars. Of course, we have to keep in mind that any automatic procedure will at
some point and to some extent still require some human intervention. However, this should
be reduced as much as possible to guarantee a global and objective solution.
Mokiem et al. (2005) put effort in writing an automated tool for the spectroscopic analysis
of massive stars. Their method was developed for the quantitative analysis of O- and early
B-type stars with stellar winds. Our aim is to develop a complementary method to cover
the complete B-type spectral range. Mokiem et al. (2005) opted for the genetic algorithm
approach, which is the most accurate one, but is also very time-consuming. In first instance,
they analysed 7 O-type stars in the young open cluster Cyg OB2 and 5 other Galactic stars
with their method (Mokiem et al. 2005). They also applied their automatic method to 31 O
and early B-type stars in the Small Magellanic Cloud (Mokiem et al. 2006), and 28 O and
early B-type stars in the Large Magellanic Cloud (Mokiem et al. 2007). Essentially, Mokiem
et al. (2005) linked FASTWIND to the genetic algorithm based optimisation routine PIKAIA
from Charbonneau (1995). They determined the stellar parameters by using continuum normalised hydrogen and helium lines, assigning weights to them to account for line blends
with species that are not taken into account, or to account for potential problems from the
theoretical side. For example, the He I singlets are predicted weaker by CMFGEN than by
FASTWIND in the temperature range between 36,000 and 41,000 K for dwarfs and between
31,000 and 35,000 K for supergiants (see Najarro et al. 2006 for the origin of this discrepancy). Mokiem et al. (2005) accounted for this discrepancy by introducing a lower weighting
factor for these He I singlets for stars in these ranges. Going down in temperature from the O
and early B-type stars towards the mid and late B-type stars, requires other diagnostic lines
besides hydrogen and helium. We need to complement these with (continuum normalised)
silicon lines to determine the effective temperature. However, the inclusion of the Si lines
in the genetic algorithm approach is not straightforward. Silicon is structurally much more
complicated than the H and He atoms, with many more transitions to account for. Moreover,
it simultaneously depends on many parameters. A genetic algorithm approach would require
a sophisticated weighting scheme for the different Si lines and frequency points. For reasons
discussed in Section 3.1, we have opted for a grid-method as an alternative approach.
Of course, we are well aware of the disadvantages inherent to this method:
1. Codes are developing quickly and continuously. Since the method makes use of a predefined grid, this implies that, as soon as a new version/update of the code is released, the
grid and consequently the derived fundamental parameters are out-of-date. This can happen when, e.g., atomic data are improved. The advantage, on the other hand, is that the
line-profile fitting method itself will remain and the job is done with the computation of
a new grid. This is of course a huge work, recalling that, for the grid described in Chapter 3, seven months were needed to compute the grid, without even taking into account
the amount of time spent on the grid setup before and the grid checking after the actual
102
GAUDI B star sample
computation. Thanks to the fast performance of the FASTWIND code, however, this is
not really an insurmountable problem.
2. Whereas a genetic algorithm uses a dynamic grid, changing with every new generation and
evolving towards the final best fitting model, a grid-based method fully relies on a static
grid and no additional models are computed. Consequently, the quality of the final best fit
and the accuracy of the final physical parameters are fully determined by the density of the
grid. This underlines, again, the necessity of a grid which is as dense as possible, allowing
for interpolation. On the other hand, the grid-method has another plus-point speaking in
its advantage: it is fast! Whereas the analysis with the genetic algorithm approach can
easily take a few days for one target star, a well-designed grid-method will, on average,
be finished within half an hour, precisely because no additional model computations are
required.
4.4 Methodology of AnalyseBstar
In what follows, we will give a detailed description of AnalyseBstar: the preparation of all
input for the programme, the iteration cycle for the determination of each physical parameter,
the assumptions made, ... In Fig. 4.3 we present a flowchart of the full programme, with
references to the appropriate section, for the description of each step. The reader is advised
to follow this context diagram throughout the further description of AnalyseBstar in this
section. The AnalyseBstar programme is written in the Interactive Data Language (IDL),
which allows for interactive manipulation and visualisation of data. It fully relies on the
extensive grid of model atmospheres and the emerging line profiles, described in Chapter 3.
We use continuum normalised hydrogen, helium and silicon lines to derive the photospheric
properties of the star and the characteristics of its wind.
4.4.1 Preparation of the input
No automatic method can fully replace the ‘by-eye’ procedure, since there will always remain a few steps for which human intervention is required. This is mainly the case for the
preparation of the spectra and the input for the actual programme, including:
1. Observation of peculiarities
The first thing one should do when starting to analyse a star, is simply look at the spectrum and inspect whether any abnormal behaviour can be detected, such as line behaviour
resulting from nebulae, disks or other peculiarities. These are things that can only be
detected through inspection by eye.
4.4 Methodology of AnalyseBstar
Figure 4.3: Context diagram of AnalyseBstar. For a detailed explanation, see text.
103
104
GAUDI B star sample
2. Merging and normalisation of the spectra
Echelle spectra are characterised by the fact that the spectrum is cut into pieces (i.e. orders)
by a diffraction grating. In general, the different orders have a small overlapping region.
This can be used to merge the orders into one spectrum. The merging, however, is a
very delicate and non-trivial process. The edges of the orders are often noisy, due to the
decrease in sensitivity of the CCD near the edges. These parts are useless and should
be removed. Moreover, some spectral lines fall right on the edge of the spectral order,
resulting in a cut-off within the spectral line. Lines that are cut-off in this way, in the
overlap region of either of the two orders, are not reliable enough to derive accurate stellar
information from them. They can, in the best case, only be used as a consistency check
afterwards.
Once the merging process is finished, the spectra are ready to be normalised. The most
frequently applied method is by interactive interpolation of continuum regions. The spectrum should be normalised over a large wavelength coverage, where sufficient continuum
regions are present. If not, line depths can get corrupted, which will further propagate in
the fundamental parameters derived from these lines. In our case the normalisation will
especially prove to be very important for the determination of the gravity, which depends
sensitively on the Balmer line wings.
GAUDI sample: The standard FITS binary tables of the FEROS and ELODIE spectra in
the database contain information about both the normalised spectrum, resulting from the
pipeline reduction of the instrument, and the unnormalised spectrum. Since the quality
of both the merging and the normalisation turned out to be insufficient for our detailed
analyses, we have redone this time-consuming job for all spectral ranges around the diagnostic lines. This part of the work was done by T. Morel. To ensure a smooth merging,
first the order edges (where the SNR is really poor) were cut off. Then, the orders were
scaled relative to the counts in the overlap region, before both were actually merged (see
Fig. 4.4). The merging is generally much better for FEROS spectra than for ELODIE
spectra. The reason for this lies in the fact that the FEROS orders are well corrected for
the blaze function and are therefore ‘flat’, which makes a smooth connection rather easy.
On the contrary, the correct reconnection of the ELODIE orders is very hard due to the
incomplete removal of the blaze function during the standard reduction. They are concave, rather than flat, which causes a sawtooth-shaped continuum, which abruptly jumps
at the wavelengths where two consecutive orders connect. Also the above mentioned edge
effects are stronger for the ELODIE spectra than for the FEROS spectra, which leads to,
e.g., discontinuities in the blue part. This cannot be completely removed, not even after
cutting off the extremities. The continuum rectification is performed over a broad spectral
range using cubic splines or sometimes Legendre polynomials for the fit. The fitting is an
iterative process in which pixels that deviate by more than a certain threshold (indicated
by the user) from the fit are excluded at each iteration.
For the few available SARG data, only the normalised spectra have been inserted into the
database and we failed to get hands on the raw data. The continuum rectification was really
poor, especially Hβ suffers a lot from this (see Fig. 4.5). The drop-off in the line wings
is too sudden, and is not at all as expected from a smooth velocity distribution. We have
observed a similar line behaviour in the red wings of the Hα profiles for the B supergiants
HD 91024 and HD 106068 (see 2.5.2 and 2.6.2) and in the blue wing of HD 111904 (see
4.4 Methodology of AnalyseBstar
105
2.6.2). The refill of the wing of Hα can be explained by the wind. However, Hβ only
becomes severely affected when the wind is very strong (which is not the case as we
have concluded from the inspection of Hα), and the line shape is different from what is
observed here, and merely affects the line core. We think the rectification through a cubic
spline has been made too deep into the wing, resulting in the observed sharp decrease, too
narrow wings and a weaker line, since the inner part of the line will be pulled upwards.
Other line profiles are affected in a similar way. We have overplotted the observed Hβ
profile in Fig. 4.5 with what would be a more reasonable shape of the wings. Also the
wavelength coverage is insufficient for our purposes. The spectral range is confined to
˚ with a gap of 250 A
˚ between 6100 A
˚ and 6350 A.
˚ Moreover, due to bad
[4620, 6845] A
˚ and [5350, 5385] A
˚ cannot be used. This means that
merging, the regions [5760, 5800] A
almost all diagnostic He and Si lines are lost. Therefore, we excluded all SARG spectra
from our sample.
We want to put special emphasis on the fact that the huge and precise task of merging
and normalising all spectra has been done by the same person and in a uniform way using
IRAF3 , which is important for assuring a homogeneous treatment of the sample.
3. Line selection
The line selection is done in a automatic way. The lines that will be used in the spectral
analysis will be read directly from a file. It is the user who decides which lines will be
used, by removing or adding lines to this file. As a standard procedure we compute a
number of line profiles in the optical (see below) from the underlying model atmosphere,
complemented with a few lines from the Paschen and Brackett hydrogen line series in the
infrared. The latter infrared lines will not be used in this study, but can contain important
information for future users. These standard line profiles can directly be used. Computing
additional line profiles is, however, straightforward.
We distinguish between two types of lines:
– The first group of lines will be used for the prediction of the physical parameters
during the spectral line fitting procedure. These are:
Hα, Hβ, Hγ and Hδ
He I 4026, He I 4387, He I 4471, He I 4713, He I 4922 and He I 6678
He II 4541 and He II 4686
Si II 4128-4130 and Si II 5041-5056
Si III 4552-4567-4574, Si III 4716, Si III 4813-4819-4829 and Si III 5739
Si IV 4116 and Si IV 4212
Note that we have excluded He II 4200. This line can savely be used for O-type stars,
but not for B-type stars. At least for dwarfs, it only appears at temperatures above ±
25,000 K, where it is still overruled by N II, while at hotter temperatures, He II 4200
is for 50% blended by N III4 .
3 IRAF (Image Reduction and Analysis Facility) is distributed by the National Optical Astronomy Observatories,
operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the
National Science Foundation, USA.
4 See http://www.lsw.uni-heidelberg.de/cgi-bin/websynspec.cgi
106
GAUDI B star sample
Figure 4.4: Illustration of the merging and normalisation process for two orders of the
FEROS spectrum of HD 170580, concentrated around Hγ, which is used for deriving log g.
Upper panel: orders 9 (black) and 10 (dark grey) before merging. Bad parts of both orders are
cut off to obtain the merged spectrum (light grey), which has been moved downwards with
flux value 1 for reasons of clarity. Lower panel: After the merging, the normalisation is done
and a continuum rectified spectrum is obtained (black). We have overplotted the pipeline
normalised spectrum (grey), to illustrate why the renormalisation was necessary.
4.4 Methodology of AnalyseBstar
107
Figure 4.5: Example of the bad normalisation of the SARG spectra for Hβ. The observed
and continuum rectified spectrum of HD 48347 (black) is overplotted with an appropriate
theoretical spectrum, accounting for the presence of Stark-broadening and a smooth velocity
distribution (grey).
– The second, smaller group contains additional lines which will not be used in the
fitting procedure, for reasons of uncertainties, either in the theoretical predictions,
either in the observed spectrum (due to problems in merging or normalising). Their
goodness-of-fit is only checked to investigate systematic differences between theory
and observations. This can be useful input for considerations regarding further improvements of atomic data and/or atmospheric models.
These are: Hǫ, He I 4010, 4120, 4140 and Si IV 4950.
4. Signal-to-noise measurement
Since the signal-to-noise ratio SNR can vary a lot with wavelength, we have chosen to
calculate the local SNR for each line separately, by using a continuum region close to
the line. This is done manually through the indication of the left and right edges of the
continuum interval. The SNR is then given by the mean flux Fc divided by the standard
deviation σ(Fc ) of the flux in this interval, i.e., SNR = Fc /σ(Fc ).
5. Equivalent width determination
The equivalent widths of the He and Si lines are measured by a Gaussian non-linear least
squares fit to the observed line profiles, using the Levenberg-Marquardt algorithm. We
have thoroughly tested whether this was a justifiable approach, by comparing the equivalent widths resulting from a Gaussian fit from a proper integration of the line pixel values
themselves for different line profiles of different stars. It seemed that, indeed, the difference between both approaches is minimal. The measurement of the observed equivalent
widths is in this way reduced to a simple integral of the Gaussian profile.
The Gaussian absorption profile (emission profiles need a different treatment and are
therefore omitted so far) can be written as
"
2 #
1 λ − A2
,
f (λ) = 1 − A0 . exp − .
2
A1
108
GAUDI B star sample
where A0 is the line depth, A1 is the half-width-at-half-maximum and A2 is the centre of
the line profile.
The begin and end wavelength of the line can be identified in two ways: the manual
identification is done by simple mouse-clicks, the automatic identification is done by the
determination of the inflection points nearest to the laboratory central wavelength, through
the calculation of smoothed second order derivatives. We have developed a routine which
is very similar to the underlying technique in the ARES code of Sousa et al. (2007). The
ARES code (Automatic Routine for line Equivalent widths in stellar Spectra5 ) has recently
been developed for the automated determination of equivalent widths of the absorption
lines present in stellar spectra, and has been tested mainly for cool stars for measuring
the EW of Fe lines. The basic principles of our code are the same as ARES. However,
slight modifications have been carried out to meet our own specific requirements. The
ARES code, e.g., uses a second order polynomial fit to determine the position of the local
continuum over a small wavelength range around the, already rectified, line profile, for
which the EW is required. We considered a first order polynomial (i.e. a linear fit) to
be a better approach, since a second order polynomial might over- or underestimate the
equivalent width of the line.
An automated method for the EW determination is powerful for spectra with a very high
SNR. However, for our sample of GAUDI stars, it appeared insufficient due to the noise
level of our data. Therefore, for our sample, we kept the manual method, so that we still
have control on this process. To simplify the identification of the begin and end wavelength of the line by manual mouse-clicks and to avoid measuring the EW of other strong
lines which lie, by coincidence, very close to the desired line profile (i.e. to avoid misidentification), we have indicated all transitions of the considered ion (i.e., all transitions for
Si, when a Si line is considered, or, all transitions for He, when a He line is considered),
which have a gf -value6 above 10−3 (see upper panel Fig. 4.6). This is especially important when the considered lines are very weak.
For He I lines, a Gaussian profile is not always the best approximation due to a strong
forbidden component in the blue wing of certain lines. Therefore, in the case of He I lines,
the user can decide, by looking at the fit, whether the Gaussian approximation for these
lines is good enough or if (s)he prefers the equivalent width to be measured as the integral
of the observed profile itself.
The main source of uncertainty in the EW determination is introduced by the noise and
its influence on the exact position of the continuum within the noise. To account for this
we have shifted the continuum level up and down by a factor (1 ± 1/SNR), performed
a new Gaussian fit through both shifted profiles and measured the difference in EW (see
lower panel Fig. 4.6). Note that, by considering a constant factor over the full line profiles
(using the signal at continuum level), we implicitly give an equal weight to each wavelength point, even
p though, in reality, 1/SNR may be slightly larger in the cores, since it is
proportional to 1/S, with S the obtained signal at these wavelengths. So, the noise is
one source of errors to be taken into account. When the user deems it necessary, (s)he can
apply a correction for the continuum level by a local rerectification, e.g. in the case that
the normalisation performed over a large wavelength range is locally a bit offset. If the
5 http://www.astro.up.pt/∼sousasag/ares/
6 The
gf -value is the product of the statistical weight g and the oscillator strength f .
109
4.4 Methodology of AnalyseBstar
offset is different on the blue and red side of the line, then the factor should be wavelength
dependent. Therefore, we allow for a linear rerectification by a factor aλ + b. A third
factor is the error introduced by using a Gaussian fit to estimate the EW.
The total error on the EW can then be estimated as follows (since the individual error are
independent from each other):
p
∆EW = (∆EWSNR )2 + (∆EWcont )2 + (∆EWfit )2 .
∆EWSNR is the error in EW with respect to pure photon noise statistics. According to
Vollmann & Eversberg (2006) it equals
s
Fc (∆λ − EW )
∆EWSNR = 1 +
,
SNR
F
where Fc is the average continuum flux, measured outside the line, and F the average flux
in the line, which covers a spectral range ∆λ.
∆EWcont is the difference in the equivalent width (of the Gaussian fits) when the local
continuum is shifted by a factor (1 ± 1/SNR) (aλ + b). Due to the noise, the interactive
continuum identification is a non-unique and subjective process. By introducing the factor
(1 ± 1/SNR), we account for small shifts in the continuum level. On the other hand, local
rerectifications are accounted for by introducing the factor (aλ + b).
∆EWfit is the uncertainty in the Gaussian profile which includes the error on the line
depth A0 as well as on the half-width-at-half-maximum A1 in the integral.
Since the fitted EW can be written as
√
EW = − 2π A0 A1 ,
the error follows as
∆EWfit =
√ q
2π (∆A0 )2 A21 + (∆A1 )2 A20 .
6. Determination of v sin i
The determination of the rotational velocity is the last step in the preparation of the input
for AnalyseBstar. As mentioned in Section 2.4.1, we use the semi-automatic tool developed by Sim´on-D´ıaz et al. (2006) to derive reliable projected rotational velocities. The
method is based on the fact that only the rotational broadening function has zeros in the
Fourier transform, in contrast to, e.g., the broadening due to macroturbulence and the instrumental profiles, which are Gaussian. As the convolutions with the different broadening
profiles are transformed into multiplications in Fourier space, one can simply divide the
Fourier transform of the observed profile by the transform of the broadening functions that
have no zeroes. The frequency of the first minimum in the residual transform corresponds
to the rotational velocity. Care has to be taken at high frequencies (low v sin i ), since
also zeroes from the microturbulence are introduced there (Gray 1973; Sim´on-D´ıaz et al.
2006).
110
GAUDI B star sample
Figure 4.6: Illustration of the EW determination for He I 4026 (for a simulated spectrum).
Top: All He transitions in the considered wavelength interval, with a log gf down to −3,
are indicated as vertical lines. This makes it easier for the user to correctly identify the line,
which could sometimes be a problem in case of noisy spectra or very weak lines. Through
simple mouse-clicks, the user indicates the begin and end wavelength of the line. Bottom:
Within the indicated wavelength interval, a Gaussian fit to the observed line profile is made
to determine the EW of the line (thick grey profile). Also indicated are: the centre of the
line (vertical line) and the shift in flux, upwards and downwards (dashed line profiles), used
to account for the noise level in the determination of the EW. Also to these line profiles a
Gaussian fit is made, but these are omitted here for the sake of clarity.
4.4 Methodology of AnalyseBstar
111
Figure 4.7: Illustration of the determination of v sin i from the first minimum in the Fourier
transform of the selected line profile, as described by Sim´on-D´ıaz et al. (2006). The user can
make several attempts and finally decide which v sin i gives the best match. The difference
in the slope of the first decay gives an indication of the macroturbulence, which in this case
is absent. The Fourier transform of the observed line profile is indicated in black. The user
indicates the first minimum in the black profile, which gives the projected rotational velocity.
The Fourier transform of the rotational profile at this vsini is indicated in grey.
This Fourier Transform (FT) method seems to be a very powerful tool to separate the
effects of rotational broadening and macroturbulence. We refer to Sim´on-D´ıaz & Herrero
(2007) for a thorough discussion.
We use the following least blended, metallic lines for the determination of v sin i : Si II 5041,
Si III 4567, 4574, 4813, 4819, 4829, 5739, C II 4267, 6578, 6582, 5133, 5145, 5151 and
O II 4452.
The task of the user is twofold:
1) identify each line (if visible in the spectrum), i.e., indicate the line region, the continuum
level and the line centre. For the sake of comfort, the line region for each line is selected
automatically and all transitions of this certain species (Si, C or O) with a log gf value
above 10−3 are indicated. 2) select the first zero. The user can make up to ten trials and
then decide which v sin i -value is the most appropriate one (see Fig. 4.7). In this way,
we have obtained a reliable v sin i -value for each selected line. The projected rotational
velocity v sin i of the star is then calculated as the mean of these values, and the standard
deviation as its uncertainty.
112
GAUDI B star sample
4.4.2 Starting values for the fit parameters
Once the preparation is finished, we can start with the automatic procedure to derive an
accurate value for all physical parameters. This procedure follows an iterative scheme, in
which the fundamental parameters are derived in the following order:
– The effects of the effective temperature, the Si abundance and the microturbulence on the
Si lines are separated.
– The He abundance is fixed from the helium lines.
– The macroturbulent velocity is determined from well-chosen Si lines.
– The surface gravity is determined from the wings of Hγ, Hδ and, in the case of weak
winds, also Hβ.
– The wind parameters log Q and β are determined in parallel using Hα.
Obviously, to start the iterative procedure, we need a good initial guess for each of these
parameters. This initial value can either be user supplied or standard. If no value is set by the
user, then the following initial values are considered:
– The initial effective temperature will be determined from the spectral type of the star.
– We start from the highest gravity in the grid (i.e. log g = 4.5), but immediately improve
this value to a more realistic guess by running the subprocedure for the determination of
the gravity.
– The initial abundance for He is taken as solar, and for Si, we consider the ‘typical’ value
for B-type dwarfs, i.e. log n(Si)/n(H) = -4.79 (cf. Section 4.5.2).
– We initialise the wind parameters at the lowest values, which means negligible wind, at
log Q = -14.30 and β = 0.9. These are reasonable assumptions for dwarfs, which comprise
the largest part of our sample.
– The macroturbulent velocity is initially zero, whereas for the microturbulent velocity we
have taken a medium value of all possibilities, i.e., 10 km/s.
It is important to choose these initial parameters as well as possible, since they will fully
determine how fast the method will converge to the final solution. This could make the
difference between three minutes and a quarter of an hour of computation time.
4.4 Methodology of AnalyseBstar
113
4.4.3 The cycle of the combined determination of the effective temperature, microturbulence and abundances
The silicon lines serve multiple purposes. By using a well-defined cycle, which is based on
the procedure outlined in Urbaneja (2004), we are able to separate the effects of the effective
temperature, the Si abundance and the microturbulence on the line profiles, and derive an
accurate value for them. Our method automatically determines from the observed spectrum
which lines will be used for this purpose. The lines should be well visible (i.e. clearly distinguishable from the noise in the continuum) and the relative error on their equivalent width is
typically below 10 to 15%, dependent on the combination of the temperature and the specific
line considered. Only lines with a relative error of less than 75% on their equivalent width
will be considered in the following procedures. Such large errors are a rare exception and
only occur for the weakest lines, which almost disappear into the noise level.
Basically, we can discern two different cases, namely when multiple and consecutive ionisation stages of silicon are available, or when there is only one ionisation stage of silicon. Both
need a different approach, as we explain now.
1) Multiple ionisation stages of Si available
Effective temperature
Since the Si abundance affects all Si lines in a similar way, we can, at first instance, put aside
the effect of the Si abundance by considering the EW ratios of two different Si ionisation
stages. For each effective temperature point in the grid, the ratios of the observed equivalent
widths of Si IV to Si III and/or Si III to Si II are compared to those obtained by the model grid
of Chapter 3, given a certain log g and wind parameters (see Fig. 4.8). Whereas we consider
both the effects of microturbulence and silicon abundance simultaneously in AnalyseBstar,
we have separated them here in the figure (top versus bottom) for the sake of clarity in the
presentation. From this figure, it is immediately clear that a model with an effective temperature of 18,000 K is at this point of the iteration cycle the best solution. The effects of
both the microturbulence and the Si abundance become larger towards the hotter temperature
regime, due to the fact that there is hardly any Si II left in this region. The same effect occurs
towards the cool temperature regime, where there is less and less Si III. Consequently, the
observed effect on the line ratio is mainly due to the effect on one of both ionisation stages
alone. For each combination of Si IV/III and/or Si III/II, this gives us a range of ‘acceptable’
effective temperatures for which the observed EW ratio is reproduced, within the observed
errors (indicated as grey symbols in Fig. 4.8). From the different line ratios, slightly different
temperatures may arise. In the assumption that log g and the wind parameters are exact, the
combination of all these possibilities constitute the set of acceptable effective temperatures.
114
GAUDI B star sample
Microturbulence and silicon abundance
Once we have a list of acceptable temperatures, we can derive for each temperature the best
microturbulent velocity and the best Si abundance. To this end, we use the “curve of growth”
method, which describes how the line strength increases with the number of atoms producing
the line.
In what follows, we describe the full process, step by step:
Step 1: Deriving the abundances for each microturbulence, line per line
In this first step, we compare, for each Si line, the observed EW of the line to the theoretically
predicted EW, for each combination of Si abundance and microturbulence in the grid (see
Fig. 4.9). For each microturbulent velocity, we determine, by simple interpolation, which Si
abundance would reproduce the exact value of the observed EW. The upper and lower errors
for the Si abundance are derived from similar interpolation for EW ± ∆ EW. This is done
for each Si line separately. In Fig. 4.9, we only show a few selected lines. Several more are
used in case they are clearly present in the spectrum. In the example shown, no Si IV was
observed, so only Si II and Si III were considered.
Step 2: Deriving the mean abundance for each microturbulence and the slope of the
best fit
In a second step we compare, for each microturbulence, the abundances (with the derived
uncertainties) found in step 1 to the observed EW of each line (see Fig. 4.10). Since a star
can only have one Si abundance, we should find the same Si abundance from all the different
Si lines, i.e. the slope of the best fit to all the lines should be zero. From the different
abundances derived from each line, we can calculate the mean abundance.
Step 3: Deriving the ‘interpolated’ microturbulent velocity
In a third step, we investigate the change of the slope when jumping from one microturbulence
to the other (see upper panel Fig. 4.11). Through linear interpolation, we find the microturbulence for which the slope would be zero (i.e. for which the Si abundance derived from each
line separately would be the same). This is the estimated ‘interpolated’ microturbulence.
Step 4: Deriving the ‘interpolated’ abundance
Using now the relation between the microturbulence and the mean abundance derived in step
2, we interpolate to find the estimated ‘interpolated’ Si abundance at the estimated ‘interpolated’ microturbulence found in step 3 (see lower panel Fig. 4.11).
Step 5: Deriving the ‘interpolated’ effective temperature
Now that we have obtained the estimated ‘interpolated’ values for the microturbulence and Si
abundance, we can go back to Fig. 4.8 and interpolate in two dimensions to derive a value for
the effective temperature, which reproduces the observed EW. The mean of the effective tem-
4.4 Methodology of AnalyseBstar
115
peratures from each line ratio will then be accepted as a value for the ‘interpolated’ effective
temperature of the star.
Step 6: Determination of the ‘closest’ grid values
For the next steps in the automatic procedure, we will need the closest grid values to these
estimated ‘interpolated’ values rather than the real values themselves. These closest values
will constitute a new entry in the list of possibilities. If the real value for the Si abundance or
microturbulence falls exactly between two grid points, then both are added to the list.
We repeat these steps for each of the ‘possible’ grid temperatures, initially derived from
Fig. 4.8 (grey symbols), and obtain in this way a set of new possible solutions that optimally
reproduce the Si lines.
Helium abundance
For each (Teff , ξ, log n(Si)/n(H)) combination, found in the above described way, we can now
determine the He abundance, denoted as n(He)/n(H). As for the Si abundance, we have only
three different values to consider. We apply a similar method as for the Si abundance determination, in the sense that we derive the best-suited abundance from each line separately, and
take the mean value as the ‘interpolated’ value. Figs 4.12 and 4.13 illustrate this process. Note
that we intrinsically assume that the microturbulence is the same throughout the atmosphere,
i.e. that there is no radial stratification. In this way, the microturbulence, necessary to account
for the broadening in the He lines, will be assumed to be the same as the one derived from
the Si lines. Note also that the uncertainty in the derived He abundance will be larger when
the He abundance is lower than 0.10. Indeed, in this case we are namely forced to extrapolate
to a region where it is unclear how the dependence of the equivalent width with abundance
will change. The decrease towards abundances lower than solar may be steeper or follow a
logarithmic trend, in which case the predicted values would be underestimated. In this case
we can only state that the solar abundance is an upper limit of the true abundance.Values
lower than the primordial helium abundance of ∼0.10 can no longer be considered physical,
except when diffusion effects start to play a role and helium settles. In that case, the helium
abundance observed at the stellar surface could be lower. Such chemically peculiar stars are
referred to as He weak stars and are usually high-gravity objects. Helium settling was not
taken into account in our models.
2) Only one ionisation stage of Si available
If we do not have two consecutive stages of Si (e.g. for the late type stars, where we only
have Si II), we have to follow a different procedure. Fortunately for the late type stars, He I
is very sensitive to changes in effective temperature (see Fig. 2.3). Therefore, we can use the
joint predictive power of He I and Si II to derive, in first instance, the plausible values for
Teff . We compare the fits of all combinations of Si abundance, microturbulence and effective
temperature, and select only the best ones.
116
GAUDI B star sample
Figure 4.8: Synthetic simulation of the effect of the microturbulence (top) and the Si abundance (bottom) on the logarithmic EW ratios of silicon lines of different ionisation stages.
This example gives the (synthetic) EW ratio of Si III 5739 to Si II 4130 against a selected
range of effective temperatures. The ratios are evaluated for the different possibilities of the
microturbulence, for fixed Si abundance (upper panel: ξ = 3, 6, 10, 12, 15 and 20 km/s, respectively triangle up, square, diamond, circle, asterisk, triangle down) and for the different
possibilities of the Si abundances, for fixed microturbulence (log n(Si)/n(H) = -4.79, -4.49
and -4.19, respectively circle, triangle up and square). The horizontal lines show the observed EW ratio and its corresponding uncertainty region. Acceptable EW ratios, which fall
within these boundaries, are indicated in grey.
4.4 Methodology of AnalyseBstar
117
Figure 4.9: Step 1 (for each Teff and for each line): line per line, for each microturbulent
velocity, the abundance that reproduces the observed EW of this line is derived. The different
abundances for each microturbulent velocity are connected by the dashed-dotted lines. The
observed EW interval is indicated with horizontal grey lines. Crosses indicate, for each microturbulent velocity (from low EW to high EW: 3, 6, 10, 12, 15 km/s), the values derived for
the Si abundance (log Si/H by number, and the upper and lower limit), found through linear
interpolation. The displayed example is a synthetic simulation.
118
GAUDI B star sample
Figure 4.10: Step 2 (for each microturbulence): we plot the Si abundances and their uncertainties for all 6 lines for which we derived these values in step 1, as a function of the
observed EW. The least squares fit to the abundances are shown in grey. The microturbulent
velocity for which the slope of this fit is zero (i.e. equal abundance from each line) gives an
estimate of the ‘interpolated’ microturbulent velocity. The displayed example is a synthetic
simulation.
119
4.4 Methodology of AnalyseBstar
Figure 4.11: Step 3 & 4: the ‘interpolated’ microturbulent velocity and the accompanying
‘interpolated’ Si abundance (indicated in grey) are derived from the position where the slope
of the best fit (step 2) is zero. The silicon abundances, given in the lower panel, are the mean
abundances derived in step 2. The displayed example is a synthetic simulation.
The criterion for selection is based on the following, standard, expression of the loglikelihood
function (see Decin et al. (2007) for a broader theoretical explanation)
l≡
NX
lines
j=1
"
√
1
−ln(σ j ) − ln( 2 π) −
2
EWobs,j − EWsyn,j
σj
2 #
.
(4.1)
In essence, this formula gives the squared deviations between the observed and synthetic
equivalent widths (EWobs and EWsyn ), summed over all lines considered (Nlines ), and taking
into account the observed error σj on the EW (assuming that the theoretically predicted value
is perfect). In this way the loglikelihood function is very similar to a simple chi-square
criterion. The model for which the loglikelihood distribution reaches its maximum, will fit
the observed data best. The distribution allows for a certain deviation from this maximum,
by taking into account the observed errors on the equivalent width. The criterion determining
which models will be taken into account for the further analysis and which not, is set by the
number of free parameters p and the required significance level α. If we want to determine
only the effective temperature (i.e. one free parameter) from the equivalent widths of the lines
of He and Si, within a 95% confidence interval, then the criterion can be written as follows:
2 (lmax − l) ≤ χ21 (0.05)
(4.2)
with lmax the maximum loglikelihood, l the loglikelihood of any other model and χ21 (0.05)
the chi-square distribution with 1 degree of freedom, evaluated for the commonly used 5%
120
GAUDI B star sample
Figure 4.12: For each helium line, the most suited helium abundance is derived through interpolation. The theoretically predicted values are shown as filled circles, the interpolated
values as crosses. The horizontal lines indicate the observed equivalent width and the observed errors. The displayed example is a synthetic simulation.
significance level. Every Teff (and according model), for which the EWs of the synthetic line
profiles satisfy this criterion will be considered as a possible solution for Teff . Generally, in
this way, the temperature is very well constrained. We will further elaborate on the use and
suitability of this loglikelihood function in the next section.
Since we have only one stage of each ion, we have no means to derive any information
about the abundances of Si nor He. Moreover, it is no longer possible to follow the same
scheme as before, to derive the microturbulence. Therefore, we assume that we know at
least one of these three parameters. Although the He lines are not only dependent on the
He abundance, but also on the microturbulence, it are still the Si lines alone that are used
to fix the value for ξ (besides the Si abundance of course). This means that one can derive
only one parameter from the He lines, but two have to be fixed from the Si lines. Therefore,
the best approach is to assume we know one of those two parameters that are fixed from
the Si lines, and consequently derive the other two free parameters. We have chosen to fix
the Si abundance. We will consider all three possibilities for the Si abundance, although we
attach more credence to the lowest abundance (log n(Si)/n(H) = -4.79), since this appears
4.4 Methodology of AnalyseBstar
121
Figure 4.13: The real helium abundance is calculated as the mean over all lines. The abundance derived from each line separately is represented by the filled circles, while the arrows
show the derived errors. The best linear fit is indicated as the dotted line, while the horizontal
full line shows the mean value. The displayed example is a synthetic simulation.
to be the more ‘typical’ abundance for B-type dwarfs (cf. Section 4.5.2). For each initially
accepted Teff (from the loglikelihood function), and for each Si abundance, we search for the
best value for ξ from the Si lines (see Fig. 4.14), and consequently, using this microturbulence
allows us then to fix also the best He abundance (in the same way as before). Note that, again,
we assume that ξ is the same throughout the atmosphere, and thus for each type of ion. We do
not allow microturbulences higher than 20 km/s, which is, in view of the low temperatures, a
quite conservative and safe upper limit.
Since the surface gravity is closely related to the effective temperature, we will need to redetermine log g for every new Teff . However, before determining log g, we first determine the
macroturbulent velocity since this plays a certain role for the determination of log g as well,
as will be shown in Section 4.4.5.
4.4.4 Macroturbulence vmacro
The effect of the macroturbulence on the line profile is similar to that of the rotation, in this
sense that it only changes the shape of the line profile, but not the strength (i.e. the EW). The
macroturbulence considered is characterised by a radial-tangential model (Gray 1975). In this
model the velocity field is modelled by assuming that a fraction (Ar ) of the material is moving
radially while the complementary fraction moves tangentially (At ). In our model for vmacro ,
we use the standard assumption that the radial and the tangential components are equal in
122
GAUDI B star sample
Figure 4.14: Once a set of possible effective temperatures is fixed through the application
of the loglikelihood criterion, we can determine for each of these temperatures (in this case:
10,500 and 11,000 K) the right microturbulence, under the assumption that we know the Si
abundance. Resulting combinations are shown for log n(Si)/n(H) = -4.19 (open circles), 4.49 (grey filled circles) and -4.79 (black filled circles). The displayed example is a synthetic
simulation.
fraction. Assuming additionally that the distribution of the velocities in the radial and the
tangential directions are Gaussian and that the specific intensity is independent of the angular
coordinate, the effect of macroturbulence can be introduced in the synthetic line profiles by
means of a convolution with a function, which consists of a Gaussian and an additional term
proportional to the error function (see equation (4) in Gray 1975). Only weak lines should
be used to determine vmacro . For B-type stars, silicon lines are the best suited, since they
most obviously show the presence of vmacro in their wings. The user decides which Si lines
should be used throughout the full process. To determine the strength of vmacro , we convolve
the Si profiles of the (at that time) best fitting model with different values of vmacro 7 . The
consecutive values considered for the convolution are chosen using a ‘bisection’ method,
with initial steps of 10 km/s. The bisection continues as long as the stepsize is larger than
or equal to 0.1 km/s, which will be the final precision of the vmacro determination. For each
considered vmacro -value, we compute its loglikelihood in order to quantify the difference
between the wings of the obtained synthetic profile and the observed line profile, as follows:
"
2 #
n
X
√
1 yi − µi
,
(4.3)
−ln(σ) − ln( 2 π) −
l≡
2
σ
i=1
with n the considered number of wavelength points within the line profile, yi the flux at
wavelength point i in the wings of the observed line profile, µi the flux at point i of the
synthetic profile (convolved with the vmacro under consideration), and σ the noise, i.e. 1/SNR,
7 At
this point the profile is already convolved with the appropriate rotational and instrumental profiles.
4.4 Methodology of AnalyseBstar
123
Figure 4.15: Determination of the macroturbulent velocity of β CMa for three different Si
lines. The panels should be read from left to right as follows.
Left: The left and right edge of the line region (i.e. outer most dots at continuum level) are
indicated by the user at the start of the procedure. The blue and red edges (grey vertical lines)
are determined as the point where the flux is half of the minimal intensity. They determine
how far the wings extend towards the line centre. The black dots represent the points that are
finally used to determine the macroturbulence.
Middle: Obtained loglikelihood as a function of vmacro (see Eq. (4.3)). Note the broad maximum in the loglikelihood distribution, which gives rather large error bars.
Right: The fit with the best vmacro is shown as the grey profile.
of the considered line profile. The considered wings extend from the position where the
intensity is half as deep as the core, until the inflection points at continuum level (which were
identified during the EW determination).
The same procedure is repeated for all considered line profiles, and the mean, derived from
the different lines, gives then the final vmacro . The error is set by the standard deviation. From
Fig. 4.15, it will be clear that the shape of the loglikelihood distribution near the maximum can
be quite broad, which indicates that the errors in vmacro with respect to the loglikelihood are
significant. This might be related to the fact that, so far, we have used only the information
contained in the line wings to determine vmacro . It needs to be investigated in the future
whether the errors can be reduced when using also the line core information.
124
GAUDI B star sample
4.4.5 Surface gravity log g
The surface gravity log g is the result of fitting the Balmer lines. We will mainly use Hγ
and Hδ for this job, possibly complemented with Hβ in case of weak winds (i.e. if log Q ≤
−13.80). Since Hǫ is somewhat blended, and the merging is not always reliable for this region
(see earlier), this line will not be used as a primary gravity determinator, but will only serve
as an additional check (together with Hβ, in case of stronger winds). Since we want to define
the profile only from the extreme wings down to the strongest curvature (see discussion of
the different contributions below), we have to account for certain mechanisms, which affect
the core of the Balmer lines to some extent. Therefore, we will have to exclude the inner ∆λ
which will be the minimum of the following three ∆’s:
1. ∆1 = (λ1 - λ0 ) with λ0 the rest wavelength and λ1 the wavelength where the intensity is
half as deep as the core. We do this on both sides and use the minimum. This gives a good
guess about the position where under typical conditions the curvature is strongest.
2. ∆2 = λ0 . (v sin i + vmacro )/2c (with c the speed of light) takes into account the influence
of the rotational and macroturbulent velocity on the line core. The factor 2 enters due to
the fact that approximately only half of the width is folded into the core.
3. ∆3 = 3 vtherm λ0 /c represents the effective thermal Doppler core width of the Stark profile. vtherm combines the effect of the thermal velocity vth and the microturbulence ξ, and
can be written as
q
2
2 + ξ2
with
vth
= 2 kB Teff /mH
vtherm = vth
in which kB is the Boltzmann constant and mH the proton mass. The factor 3 in ∆3 arises
from the fact that the Doppler profile
√
2
1
e−(∆λ/∆λD ) ,
π∆λD
with
∆λD
λ0
=
c
r
2 kB Teff
+ ξ2
mH
for ∆λ/∆λD = 3 is smaller than the corresponding Stark profile, such that the Stark
wings dominate over pure Doppler broadening (as discussed, e.g., in Section 4.1 and
Fig. 3 of Repolust et al. (2005)).
Once we have removed the inner part of the line profile, we can compare the observed Balmer
lines with the synthetic line profiles, by considering all possibilities for log g which are available at this given effective temperature gridpoint. We decide which gravity is the best, by
maximising, for each log g, the loglikelihood
2 #
nj "
NX
lines X
√
1 yj,i − µj,i
−ln(σj ) − ln( 2 π) −
l≡
,
(4.4)
2
σj
j=1 i=1
4.4 Methodology of AnalyseBstar
125
with Nlines the number of Balmer lines used, nj the considered number of wavelength points
within each line profile, yj,i the observed flux at point i of line j, µj,i the expected (i.e.
theoretically predicted) flux value, and σj the noise of each separate Balmer line.
With this procedure, we came across the following problem: Hγ and Hδ (our main gravity
indicators) can sometimes be seriously affected by a large amount of blends. If one ignores
this, one will find too high a gravity, since the loglikelihood function will take all these blends
into account (see Fig. 4.16). Therefore, we have to cut away also these blends and make the
wings a bit smoother. In order to avoid a manual intervention, we have also automated this
procedure. This is done by performing a sigma-clipping algorithm which keeps only those
points that have a flux higher than both neighbouring points (the so-called ‘high points’) and
where, at the same time, the difference is less than 1/SNR. After removing all other flux
points, only the so-called ‘line continuum’ is left (see middle panel Fig. 4.16). We interpolate
this continuum to obtain the flux at all original wavelength points, and the innermost part of
the line is added again (i.e. the part between the last blend in the blue wing and the first blend
in the red wing, but obviously still without the central core). Also all original flux points
which deviate by less than 1/SNR from the interpolated ‘line continuum’ are included again.
In this way we keep only the flux points which really determine the shape of the wing to
fit the synthetic Balmer wings, and we are able to determine an accurate value for log g (see
right panel Fig. 4.16). We realize that some points that may be marked as local continuum in
this way, may still be lower than what the real local continuum level would be, because of the
transition of one blend into another. However, this will only be the case for very few points,
which will give no significant weight to the least squares fit.
The finally accepted surface gravity will be this log g which gives the best fit to all selected
Balmer lines simultaneously. Half the gridstep in log g could thus be considered a good error
estimate. However, to account for the coupling with Teff , we consider 0.1 dex as a more appropriate error. Due to stellar rotation, the gravity derived from the Balmer wings is actually
not the real gravity, but the effective gravity. If one wants to calculate e.g. the mass of a star,
one should correct this value with a factor (v sin i )2 /R∗ , as explained in Section 1.5.2.
4.4.6 Wind parameters
A change in mass loss rate will mainly affect the shape of Hα. For cool objects and weak
winds, Hα is nearly ‘photospheric’ and will, in essence, be an absorption profile, with more
or less symmetric components, since the ‘re-filling’ by the wind is low. When the mass
loss rate increases, the particles will scatter more photons in the direction of the observer.
The absorption component becomes stronger, but especially the emission component gains
significant strength. Therefore, in the case of stronger winds, Hα will take the shape of a
typical P Cygni profile. For hot objects, where Hα is dominated by recombination processes,
and high mass-loss rates, the profile may even appear in full emission.
The wind parameter β determines the velocity law, which directly influences the density. The
red wing of Hα is well suited to determine β, since it is formed in the receding part of the
126
GAUDI B star sample
Figure 4.16: Illustration of the effect of blends on the determination of log g, and the clipping
algorithm to improve the fit, for β CMa. Left: The many blends in the wings of the profile
prevent an accurate least squares fit and lead to too broad wings, implying too high a gravity.
Central: From the observed line profile (black), only the ‘high points’ (grey dots), which
deviate less than 1/SNR from their neighbours, are kept, complemented with the full inner
part of the line, i.e., the part between the last blend in the blue wing and the first blend in
the red wing (indicated as grey vertical lines). Right: After the removal of all blends through
sigma clipping (see text), we are able to obtain a good representation of the line wings (grey
profile) and, therefore, to derive a very accurate value for the gravity.
(almost) complete wind, whereas the blue wing (either in emission or absorption) forms in
the small volume just in front of the stellar core, both by absorption and emission effects.
Therefore, errors due to certain assumptions in the standard model (rotation, clumping, etc.)
will be amplified in the blue wing, whereas in the red wing, which consists of the emission
from a much larger volume, most errors are more likely to cancel out. For a fixed mass loss, a
‘slower’ velocity law (i.e. a higher β value) will result in higher densities in the lower atmosphere, close to the star. This enlarges the number of emitted photons with velocities close
to the line centre, resulting in more emission. Around the central wavelength, the absorption
component of the line profile refills and the emission component becomes stronger. Therefore
the slope of the red wing of the P Cygni profile becomes steeper. In this sense, the steepness
of the red wing is a measure for the value of β.
We have estimated the wind parameters log Q and β by comparing the observed Hα profile,
with the different synthetic Hα profiles by making combinations of log Q and β. We have
decided to make the determination of the wind parameters not too sophisticated, by simply
using the best values for log Q and β available in the grid, without interpolating and further
refining. The reason for that is that almost all stars in our sample are dwarfs (they are chosen
like this by the CoRoT-team). In this sense we do not need a very sophisticated analysis of the
wind, since they have at most a weak wind. At the start of the programme, the user is asked
to select the type of Hα profile. Choices are: absorption profile (’a’), intermediate emission
profile with only a small emission peak (’i’) or a stronger emission profile (’e’). Then, for
each log Q-value, the corresponding best β-value is the one that gives the best fit to the red
wing (see above). In the case of a thin wind (i.e. absorption profile) the red wing extends
from the lower most flux point up to the line continuum. In case of a thick wind (i.e. emission
profile), the wing goes from the top of the emission peak, until continuum level. Once we
have selected for each log Q the best β, we decide on the best combination by considering
the best fit to the entire profile.
4.4 Methodology of AnalyseBstar
127
Figure 4.17: (Synthetic simulation) For each different wind strength parameter log Q (7 values, from left to right and from top to bottom: log Q =
−14.30, −14.00, −13.80, −13.60, −13.40, −13.15, −12.70), we search for the best β (in
each panel, the 5 different values for β are indicated: 0.9, 1.2, 1.5, 2.0, 3.0) by comparing
only the red wing (grey part of the profile). Then the synthetic profiles of each best (log Q,
β) combination are compared to the entire Hα profile (bottom right).
128
GAUDI B star sample
4.4.7 Final remarks
At the end of this iteration cycle an acceptance test is performed. When the method can run
through the whole cycle without needing to update any of the fit parameters, the model is
accepted as a good model, and gets the flag ‘2’. If a better model was found (i.e. when one
or more parameters have changed), the initial model is rejected and gets flag ‘0’, whereas the
improved model is added to the list of possibilities, and gets flag ‘1’ (cf. ‘models to check’ in
the flowchart diagram on page 103). As long as the list with possibilities with flag ‘1’ is not
empty, the parameters of the next possibility are taken as new starting values. Models that
did not converge properly are automatically skipped.
Line profile fits are always performed allowing for a small shift in wavelength, or radial
velocity (5 wavelength points in either direction), to prevent flux differences to add up in case
of a small radial velocity displacement.
To the major part, this method uses only the EWs and the far wings of the profiles. Cores are
not used, except for wind features. The profile shape only plays a role for deriving v sin i
and vmacro .
4.5 Formal tests and comparison for high-resolution spectra of selected β Cephei and SPB stars
4.5.1 Formal tests of convergence on synthetic spectra
Before applying AnalyseBstar to real, observed spectra, we first test whether the method is
able to recover the input parameters of synthetic spectra. For this purpose, we have created
several synthetic datasets in various regions of parameter space. We will not dwell on discussing all of them here, but we have chosen three specific examples, each representative for
a different type of star:
- Dataset A: a B0.5 I star with a rather dense stellar wind and a strongly enhanced helium
abundance,
- Dataset B: a B3 III star with a weak stellar wind, with ‘typical’ B star values for the abundances of He and Si (cf. Section 4.5.2),
- Dataset C: a B8 V with a very thin stellar wind.
Each dataset has been convolved with both a rotational and a macroturbulent broadening profile. The projected rotational velocity adopted is 50 km/s for each dataset, characteristic for a
‘typical’ slow rotator. An additional Gaussian, instrumental convolution was carried out, chosen in such a way that it is representative for the high resolution of the FEROS and ELODIE
spectra in our sample. Artificial (normally distributed) noise was added to mimic a real spectrum with a mean local SNR of 150, which is more or less the minimal local SNR obtained
for the GAUDI sample. Finally, we have gone through the full process of the preparation of
the spectra as if it were real data, i.e., defining the EW of all available lines, measuring the
4.5 Formal tests and comparison
129
SNR to account for the errors on the EW and fixing the projected rotational velocity. Since
we are dealing with synthetic data, obviously no normalisation was required.
Table 4.1 lists the input parameters for the three synthetic datasets as well as the output parameters, obtained from the application of AnalyseBstar. First we list the derived ‘interpolated’
values and in the third column, we display the parameters of the closest (best fitting) grid
model, which in the ideal case should be exactly the same as the input parameters. This
depends, however, on how well the equivalent widths were measured, and minor deviations
occur as expected.
Besides these three cases, in which we started from a synthetic model which is one of the grid
points, we have additionally created three synthetic datasets for which the parameters lie in
between different grid points. The input parameters and the results of the analysis have been
added to Table 4.1 as datasets D, E and F.
Although the models given in this table were selected as the best fitting models, the programme additionally came up with some other models, which agree with the input model
within the errors, introduced by the artificial noise and the errors in the determination of the
equivalent widths. In all six cases, the input parameters are well recovered. We considered
numerous other test cases, which are not included in this text, but for which the input parameters were equally well recovered. Slight deviations from the input parameters are as expected.
Also in the cool temperature domain (datasets C and F), where an alternative method is required, making an assumption about the Si abundance, we are still able to deduce reliable
parameter estimates, at least for the synthetic models. Also the derived v sin i -value, which
is difficult to disentangle from the macroturbulent velocity, agrees very well with the inserted
values, irrespective of vmacro . This shows that the Fourier Transform method of Sim´on-D´ıaz
et al. (2006) indeed allows to separate both effects. All together, this gives us confidence
that our method is working fine and that the procedure will also be able to recover the true
physical parameters from real spectral data where the input is known to the star only.
Remark on the macroturbulent velocity
The macroturbulent velocity is the only fit parameter for which significant deviations arise,
if vmacro is low (< 30 km/s). Larger macroturbulences are nicely recovered, however. Especially for datasets C and E, where vmacro is only 10 km/s, the deviation is striking (of the
order of 20 to 30 km/s). To understand this discrepancy, we first had a look at the line profiles, which show that the fit is really good for this vmacro (see Fig 4.18). After verifying the
fit quality, we also verified that the discrepancies do not arise from our applied procedure,
by performing several tests on simulated spectra, in which we left vmacro as the only free
parameter. In all cases, the macroturbulent velocity was recovered very well, with deviations
within 5 km/s.
The (sometimes significant) deviations in vmacro can be understood in two ways. On the one
hand, slight deviations are expected, as the closest grid model and the observed line profiles
do not always match perfectly (due to small differences between the interpolated ‘real’ values
and the closest grid values). On the other hand, at present, we do not include the full line pro-
Input parameters for the synthetic models (IN) are compared to the actual output parameters (OUT) obtained through the application of AnalyseBstar to these synthetic data. Datasets
A, B and C are models which lie exactly on a grid point, whereas dataset D, E and F are models which lie in between the gridpoints. The ‘interpolated’ values derived for the parameters, as well as the
best fitting (i.e. closest) grid model are given. The parameters in italic are those for which no real interpolation was made, but for which the most representative grid value is chosen. We additionally
list the number of different models checked throughout the procedure and the time (in seconds) elapsed during the run of AnalyseBstar. The longer computation times for the cooler models are due to
the applied method in this range, as will be discussed in Section 4.6.3.
IN
OUT
interpolated
Dataset A
OUT
grid
IN
OUT
interpolated
Dataset B
OUT
grid
IN
OUT
interpolated
Dataset C
OUT
grid
23,000
2.7
C
2.0
0.20
-4.49
10
50
20
142
10
22,600
2.7
C
2.0
0.18
-4.49
10.2
48 ± 2
29 ± 4
23,000
2.7
C
2.0
0.20
-4.49
10
48 ± 2
29 ± 4
18,000
3.3
A
1.5
0.10
-4.79
15
50
30
494
21
18,100
3.3
b
1.2
0.08
-4.81
15
48 ± 4
33 ± 2
18,000
3.3
b
1.2
0.10
-4.79
15
48 ± 4
33 ± 2
13,000
4.2
O
0.9
0.10
-4.79
6.0
50
10
2344
128
13,000
4.2
a
0.9
0.08
-4.79
7.3
51 ± 1
30 ± 1
13,000
4.2
a
0.9
0.10
-4.79
6.0
51 ± 1
30 ± 1
fit parameter
Teff (K)
log g (cgs)
log Q (char)
β
n(He)/n(H)
log n(Si)/n(H)
ξ (km/s)
v sin i (km/s)
vmacro (km/s)
time (s)
6 checked models
=
Dataset D
21,120
3.98
b
1.42
0.14
-4.85
12
50
40
174
10
20,985
4.0
b
2.0
0.11
-4.80
15.7
52 ± 6
33 ± 5
21,000
4.0
b
1.2
0.10
-4.79
15
52 ± 6
33 ± 5
15,100
1.83
A
2.8
0.08
-4.23
11
50
10
392
26
15,030
1.80
A
3.0
0.09
-4.16
10.6
50 ± 7
39 ± 7
Dataset F
15,000
1.8
A
3.0
0.10
-4.19
10
50 ± 7
39 ± 7
11,880
2.43
a
1.02
0.18
-4.49
9
50
70
1254
95
11,500
2.3
O
0.9
0.20
-4.19
6.3
54 ± 1
67 ± 7
11,500
2.3
O
0.9
0.20
-4.19
6
54 ± 1
67 ± 7
GAUDI B star sample
Teff (K)
log g (cgs)
log Q (char)
β
n(He)/n(H)
log n(Si)/n(H)
ξ (km/s)
v sin i (km/s)
vmacro (km/s)
time (s)
6 checked models
=
Dataset E
130
Table 4.1:
4.5 Formal tests and comparison
131
Figure 4.18: Example of the fit quality for dataset E. The slight discrepancy in the cores
of the Si lines arises from the difference in Si abundance between the interpolated value
(log n(Si)/n(H) = -4.16) and the closest grid value (log n(Si)/n(H) = -4.19)
132
GAUDI B star sample
file information in the prediction of the macroturbulence. As discussed in Section 4.4.4.,
line wings alone are perhaps not sensitive enough to changes in vmacro , and core information
may have to be inserted. This has not been implemented yet, since this requires intensive
investigation and testing, beyond the scope of the present work. We will tackle this problem
in the near future though.
As a conclusion, the values derived for vmacro should be treated with caution, since they
may not be representative for the real value. This does not affect the derivation of the other
parameters, however, since a convolution with the vmacro -profile preserves the EW of the
lines. It is the EW, and not the line profile shape, which is used for the derivation of Teff ,
ξ and log n(Si)/n(H). That other parameters remain unaffected, is clear from the very good
agreement between their input and output values in Table 4.1.
4.5.2 Testing AnalyseBstar on high-resolution spectra of selected pulsators
After we convinced ourselves that AnalyseBstar converges towards the optimum solution
and that it is able to recover the input parameters of synthetic spectra, we can go one step
further and test the procedure on real spectra. Before we move on to the GAUDI dataset,
however, we have tested our method on a selected sample of high-quality, high-resolution
spectra of pulsating B stars, mainly β Cephei and SPB stars, but also of a few other candidate
pulsators which are not yet confirmed. The (mean) spectra of these stars have a very high
signal-to-noise ratio, attained through the addition of a large number of individual exposures
(see, e.g., Morel et al. 2006). The spectra were first put in the laboratory rest frame, before
the weighted (by the SNR ratio) average of the whole time series was created. Because of
the high quality of these spectral time series, the average spectra are ideally suited for testing
AnalyseBstar. The β Cephei stars and two of the SPB stars have been analysed in detail by
Morel et al. (2006)/Morel (2007) and Briquet et al. (2007), respectively. These authors have
used the latest version of the NLTE line formation codes DETAIL and SURFACE (Giddings
1981; Butler 1984), in combination with plane-parallel, fully line-blanketed LTE Kurucz atmospheric models (Kurucz 1992), to determine the atmospheric parameters (negligible wind)
and element abundances.
Figs 4.19 and 4.20 show the fit quality obtained for a representative β Cephei and SPB star.
Tables 4.2, 4.4 and 4.3 list the results of our analysis for respectively the β Cephei stars, the
SPBs and the candidate pulsators. Whenever possible, we compare with the physical parameters found by Briquet & Morel (2007) for the SPBs and Morel (2007) for the β Cephei stars
and the candidate pulsators. In Fig. 4.21 we compare the effective temperatures and the surface gravities, as derived by AnalyseBstar, with the values derived with DETAIL/SURFACE.
The shaded area around the one-to-one relation indicates the uncertainty in Teff , as well as
in log g, on the derived DETAIL/SURFACE values, i.e. 1,000 K in Teff and 0.15 in log g.
The FASTWIND nominal errors (1,000 K in Teff and 0.10 in log g) stay within this area.
This figure shows the perfect agreement between the two independent parameter determinations. FASTWIND and DETAIL/SURFACE clearly give consistent results, which strongly
4.5 Formal tests and comparison
133
indicates that AnalyseBstar is indeed converging properly to the correct values, also for real
stars. In the lower panel of Fig. 4.23, we show the position of the pulsators within the HRD.
All β Cephei stars nicely fall in the corresponding instability strip, whereas there are a few
SPB stars which seem hotter than what is expected from the SPB instability strip.
The v sin i values derived by Morel (2007) are typically slightly larger than those derived
by us. This is readily understood as their value includes not only the rotational, but also
the macroturbulent velocity (as macroturbulence is not included in DETAIL/SURFACE). In
those cases where we derive vmacro = 0 km/s, the v sin i values are in perfect agreement.
As a conclusion, 13 out of 31 sample stars seem to display no or only negligible macroturbulence. Even though we stated earlier that we should be careful with the derived macroturbulences (see Section 4.5.1), we believe that these ‘zero’ values are quite representative,
and that indeed only negligible macroturbulence is present in these stars. Indeed, this ‘zero’
macroturbulence results as the average value from all considered Si lines of the star.
Niemczura & Daszy´nska-Daszkiewicz (2005) derived basic stellar parameters (besides metalicities and angular diameters, also effective temperatures and surface gravities) for all β
Cephei stars observed by the International Ultraviolet Explorer (IUE) satellite mission. The
parameters are derived by means of an algorithmic procedure of fitting theoretical flux distributions to the low-resolution IUE spectra and ground-based spectrophotometric observations.
In this way they provide us with a means to compare the basic stellar properties derived from
an optical analysis (this work) with the ones derived from UV studies, at least for those β
Cephei stars in common. Once again, we can be delighted: Fig. 4.22 shows that most temperatures and gravities agree within the errors (grey zone). Only for ξ 1 CMa, the temperature
from the optical analysis is at least 2000 K higher than the one derived from the low-resolution
UV spectra. The gravity, however, agrees really well (3.89 compared to our value of 3.8). In
view of the spectral type, and the effective temperature and gravity derived by Morel (2007),
we can be quite convinced that, in this case, the UV value is offset.
Depletion in Si: real or artefact?
For the majority of our sample objects we find a depletion in silicon, with abundances generally between the solar and depleted grid point. We find a mean abundance of log n(Si)/n(H) =
−4.75 ± 0.06 for the β Cephei stars in our sample. The upper panel of Fig. 4.23 shows the
good agreement with the values derived by Morel (2007) for the sample of β Cephei stars.
Indeed, Morel (2007) also finds abundances which are typically lower than solar (on the average, log n(Si)/n(H) = −4.81 ± 0.08), which gives us confidence that the depletion in these
stars might be real. The abundance of the non-confirmed pulsators, however, leads to a mean
value of log n(Si)/n(H) = −4.65 ± 0.16, which is intermediate, and thus consistent with both
the solar and depleted Si abundance. We have verified that no significant trend exists between
the Si abundance on the one hand and the effective temperature or the microturbulence on the
other hand, since these parameters are derived from the same lines.
134
GAUDI B star sample
Figure 4.19: Example of the obtained fit quality for a hot giant (β Cephei) star: 15 CMa
(B1 III).
4.5 Formal tests and comparison
135
Figure 4.20: Example of the obtained fit quality for a cool main sequence (SPB) star:
HD 215573 (B6 IV). We have left out all line profiles which had zero EW. Note that the
Si II line prediction appears slightly stronger and broader than observed. This is due to the
fact that the real microturbulence is actually lower than the displayed 3 km/s, which is the
lowest grid value available at this grid point.
136
GAUDI B star sample
Table 4.2: Basic fundamental parameters of the sample of β Cephei stars, as derived from
AnalyseBstar. We list the values of the closest grid model, except for the Si abundance, where
the interpolated value is given. In italic, we display the results from the DETAIL/SURFACE
analysis by Morel (2007). Column ’He’ indicates the He-abundance n(He)/n(H), and column ’Si’ indicates the Si abundances log n(Si)/n(H). For a summary of the errors, see Section 4.5.3.
name SpT
Teff
(K)
log g
(cgs)
He
3.80
3.75
3.40
3.60
3.80
3.70
3.50
3.50
3.60
3.65
3.90
3.80
3.80
3.75
3.70
3.75
0.10
ξ 1 CMa B0.5-B1IV 27000
27500
15 CMa B1III
24000
26000
β Cep B1IV
25000
26000
β CMa B1.5III
24000
24000
12 Lac B1.5IV
24000
24500
δ Ceti B1.5-B2IV 23000
23000
γ Peg B1.5-B2IV 23000
22500
ν Eri B1.5-B2IV 23000
23500
0.10
0.10
0.10
0.10
0.10
0.10
0.10
ξ log Q
(km/s) (char)
Si
-4.69
6
-4.87 ± 0.21 6 ± 2
-4.79
12
-4.69 ± 0.30 7 ± 3
-4.70
6
-4.89 ± 0.23 6 ± 3
-4.76
15
-4.83 ± 0.23 14 ± 3
-4.70
10
-4.89 ± 0.27 10 ± 4
-4.80
6
-4.72 ± 0.29
1+3
−1
-4.84
<3
-4.81 ± 0.29
1+2
−1
-4.73
10
-4.79 ± 0.26 10 ± 4
A
O
A
O
a
a
O
O
-
β
v sin i vmacro
(km/s) (km/s)
2.0
9±2
10 ± 2
3.0 34 ± 4
45 ± 3
3.0 26 ± 3
29 ± 2
2.0 19 ± 4
23 ± 2
2.0 39 ± 12
42 ± 4
2.0 14 ± 2
14 ± 1
2.0 10 ± 1
10 ± 1
2.0 21 ± 3
36 ± 3
11
38
24
20
46
0
0
39
-
Table 4.3: Same as for Table 4.2, but now for the hot B pulsators which are not yet confirmed.
For τ Sco, we additionally compare with the values found by Mokiem et al. (2006), using his
genetic algorithm approach. These are indicated with the superscript M .
name
SpT
θ Car
B0Vp
Teff
(K)
31000
31000
τ Sco
B0.2V
32000
31500
31900M
HD 36591 B1V
26000
27000
ǫ CMa
B1.5-B2II/III 22000
23000
γ Ori
B2II-III
22000
22000
ι CMa
B2.5Ib-II
17000
17500
log g
(cgs)
He
4.20
4.20
4.20
4.05
4.15M
3.90
4.00
3.20
3.30
3.60
3.50
2.60
2.75
0.15
0.10
0.12M
0.10
0.10
0.10
0.10
-
Si
-4.38
-4.57 ± 0.23
-4.55
-4.76 ± 0.14
-4.73
-4.75 ± 0.29
-4.69
-4.77 ± 0.24
-4.79
-5.00 ± 0.19
-4.78
-4.82 ± 0.31
ξ log Q β
(km/s) (char)
10
12 ± 4
6
2±2
M
10.8
6
3±2
15
16 ± 4
10
13 ± 5
15
15 ± 5
a
B
O
b
O
b
-
v sin i vmacro
(km/s) (km/s)
1.5
108 ± 3
113 ± 8
1.2
10 ± 2
8±2
M
0.8
5M
3.0
13 ± 1
16 ± 2
2.0
32 ± 2
28 ± 2
3.0
46 ± 8
51 ± 4
1.2
27 ± 4
32 ± 3
0
0
0
20
37
39
-
137
4.5 Formal tests and comparison
Table 4.4: Same as for Table 4.2, but for the SPBs. For two SPBs, we also show (in italic)
the results of the DETAIL/SURFACE analysis by Briquet & Morel (2007).
HD number SpT
3360 B2IV
Teff
(K)
22000
22000
85953 B2IV
20000
21000
3379 B2.5IV 20000
74195 B3IV
16500
160762 B3IV
19500
25558 B3V
17500
181558 B5III
15000
26326 B5IV
15500
206540 B5IV
13500
24587 B5V
14500
28114 B6IV
14000
138764 B6IV
14500
215573 B6IV
13500
39844 B6V
14500
191295 B7III
13000
21071 B7V
13500
37151 B8V
12500
log g
(cgs)
He
Si
3.80
3.70
3.80
3.80
4.30
3.70
4.10
4.00
4.00
3.60
3.80
4.00
3.50
3.90
3.80
3.70
3.70
3.70
3.80
0.10
-4.76
-4.72 ± 0.30
-4.84
-4.75 ± 0.30
-4.73
-4.79
-4.86
-4.72
-4.79
-4.49
-4.79
-4.79
-4.79
-4.51
-4.49
-4.39
-4.79
-4.79
-4.79
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
ξ log Q
(km/s) (char)
β
v sin i vmacro
(km/s) (km/s)
<3
1±1
6
1±1
<3
<3
3
6
6
<3
3
<3
<3
3
<3
<3
3
3
<3
3.0
3.0
2.0
1.2
3.0
1.2
0.9
1.2
2.0
3.0
3.0
1.5
0.9
0.9
2.0
1.2
1.5
18 ± 2
19 ± 1
30 ± 1
29 ± 2
48 ± 8
18 ± 2
8±2
28 ± 2
17 ± 2
17 ± 1
15 ± 2
25 ± 4
21 ± 4
21 ± 2
8±1
16 ± 1
16 ± 1
22 ± 1
20 ± 2
O
O
O
a
O
O
a
a
O
O
O
O
A
O
O
O
O
13
20
43
18
0
31
0
17
0
21
17
0
0
0
15
0
0
Although this is not the first time that such low abundances have been observed, nobody
seems to have a clear answer to the question why Si should be underabundant. Without
attempting to be complete, we list here a few literature studies, which confirm the depletion in
Si. Kilian (1992) found a mean abundance of log n(Si)/n(H) ∼ −4.75 ± 0.2 for 21 unevolved
B-type stars in the local field and the nearby associations Ori OB1, Sco Cen, and Sgr OB1.
Daflon & Cunha (2004) performed a comprehensive and homogeneous abundance analysis
for 90 non-evolved late O to early B Galactic stars. They found a mean abundance of ∼ -4.80.
Morel et al. (2006) find an unweighted mean value for the Si abundance of −4.83 ± 0.05 for
their sample of nine β Cephei stars, which are prime targets for seismic theoretical modelling.
All nine targets are early B type stars on (or close to) the main sequence. All these studies
find similar depletions for the elements C, O, Mg, Al, S and Fe for the B-type dwarfs. From
these studies, it seems clear that the Sun is indeed more metal-rich than young nearby B
dwarfs.
On the other hand, there are also studies that show the contrary, namely that there exist Btype stars which do have values consistent with the solar abundance, as shown by, e.g., Gies
& Lambert (1992) and Rolleston et al. (2000). Crowther et al. (2006), in their study of early
B-type supergiants, fixed the silicon abundance to solar values and did not report on any
discrepancies which might have resulted, so we presume that their results did not show any
evidence to reject this assumption. Also Markova & Puls (2007) found values which agree
within 0.1 dex with the solar values for seven out of eight B supergiants. For one star they
find an increase in Si with 0.44 dex.
138
GAUDI B star sample
The scatter in derived abundances is large, and each determination depends both on the specific method applied and on the assumed atomic data for silicon. N. Przybilla is currently
updating the Si model. Some preliminary studies indicate that Si III 4813-4819-4829 becomes stronger, whereas Si III 4552-4567-4574 remains more or less unaffected. We believe
that the change in the Si III 4800 triplet will not affect the Si abundance to such an extent
that this might cancel the depletion completely, certainly not because these lines are used in
parallel with several other Si lines. Moreover, only for 3 or 4 stars, lines from the Si III 4800
triplet was considered due to the fact that they are often severely affected by noise and/or
blends. Since the Si III 4500 triplet does not change significantly, we do not expect much
difference with the updated atomic data.
There might, however, exist a difference between supergiants and dwarfs, since several studies find values consistent with solar Si abundance (or, no discrepancies arose when assuming
it) for B-type supergiants, whereas there is a vast amount of literature pointing out a clear
depletion in the Si abundance for B-type dwarfs. However, a physical explanation has not
been established so far. We suggest the following one. All stars in our sample are solar
neighbourhood stars. As explained in the introduction (cf. Section 1.1), B-type stars evolve
so quickly that they generally stay near their birthplaces in the galactic plane. This means that
our sample stars must have been born in the vicinity of the Sun. Once a star is born, Si will
not be nuclearly produced until close to the supernova stage, as Si is a product of advanced
nucleosynthesis, which means that the sample stars should have Si abundances which do not
differ significantly from the solar one. The lower silicon abundance must thus be due to diffusion. Gravitational settling due to diffusion effects can indeed imply that the Si sinks towards
the stellar interior, such that the abundance observed at the surface is lower than the actual
one. This would be similar to the accumulation of iron in the driving zone of a 10 M⊙ star, as
suggested by Bourge et al. (2006), as well as an explanation of particular surface abundances
in terms of diffusive effects (Bourge et al. 2007). On the other hand, there apparently exists
a difference between main sequence and supergiant B stars regarding their surface Si abundance. We can explain this with the following scenario. During the main sequence phase,
silicon sinks due to diffusion effects. It accumulates in a region, located supposedly near the
iron opacity bump, where the radiative force overcomes the gravitational force acting on Si.
This means that Si is stuck there, which leads to the observed depletion at the surface. When
the B main sequence star leaves the main sequence, two convection zones are created. The
first one is located near the iron opacity bump, and, as Teff decreases, another one appears
near the He ionisation zone. These two zones are likely to merge and the result of this will
be a dredge-up of Si from its accumulation zone to the surface. In this way, the observed Si
abundance at the surface of B supergiants is restored to values similar to the one of the Sun.
This scenario explains both the observed depletion in main sequence dwarfs and the normal
Si abundance of B supergiants (Arlette Noels, private communication). Theoretical diffusion
computations of this scenario are not yet available.
If this scenario is true, we can be quite confident that this does not change the results presented in Chapter 2, where we adopted solar Si abundances for the analysis of the sample of
periodically variable B-type supergiants. For the majority of the SPBs, on the other hand, we
gave preference to a model with a depleted Si abundance, when a choice needed to be made,
i.e. when only Si II was available.
4.5 Formal tests and comparison
139
Figure 4.21: Comparison of the effective temperatures (top) and surface gravities (bottom)
derived with FASTWIND (AnalyseBstar) with those found by Morel (2007) and Briquet &
Morel (2007), using DETAIL/SURFACE, for the β Cephei stars (black), SPBs (grey) and
pulsators which are not yet confirmed to belong to one of these two classes (open symbols).
Supergiants are indicated as diamonds, giants as triangles and dwarfs as circles. The shaded
area around the one-to-one relation (straight line) denotes the uncertainty on the derived fundamental parameters of Morel (2007) and Briquet & Morel (2007) (i.e. 1,000 K in Teff and
0.15 in log g).
140
GAUDI B star sample
Figure 4.22: Comparison of the effective temperatures (top) and surface gravities (bottom)
derived from detailed spectroscopic analysis in the optical region (this work) with those derived from UV flux distribution fitting by Niemczura & Daszy´nska-Daszkiewicz (2005), for
the β Cephei stars in common. The shaded area around the one-to-one relation (straight line)
denotes the uncertainty region of 1,000 K in Teff and of 0.15 in log g.
4.5 Formal tests and comparison
141
Figure 4.23: Top: Position of the considered B-type pulsators in the HRD. Theoretical instability domains are the same as in Fig. 2.14. Symbols have the same meaning as in Fig. 4.21.
Bottom: Si abundances derived from AnalyseBstar, compared to those derived from DETAIL/SURFACE, for the β Cephei stars and the non-confirmed pulsator. The black dotted
lines represent the standard solar abundance of log n(Si)/n(H) = -4.49, the grey dotted lines
represent a depletion and enhancement with 0.3 dex (i.e. -4.79 and -4.19 respectively). Symbols have the same meaning as in Fig. 4.21.
142
GAUDI B star sample
4.5.3 Error estimates
Since our method is fully based on finding the best fit within a predefined grid, the precision and the accuracy of the fundamental parameters is reflected by the grid mesh. Indeed,
under the assumption that all parameters would be independent, and could be determined irrespective of all other parameters, the ‘typical’ errors would be half the grid steps. In reality,
however, the determination of most parameters is to some extent linked to the value for the
others, and their accuracy will thus depend on how well each of these can be determined.
For example, as log g and Teff are closely coupled, an error in the determination of Teff will
automatically propagate to the derived log g-value. This is why we have tried to search for
as many independent diagnostics as possible. We are confident that we have found a good
way to reduce the error maximally, by using the described scheme to decouple Teff , ξ and
log n(Si)/n(H), using two ionisation stages and a sufficient number of lines. When we have
lines of only one ion present, the errors will be much larger and will depend on the assumed
abundance.
In general, a thorough statistically valid error propagation, including the interdependency of
all fundamental parameters in this multidimensional space, is extremely complicated (especially because the errors do not propagate in a linear way), and is certainly beyond the scope
of this thesis. The best way to get an idea of the errors is to use an empirical error propagation, by varying each parameter one by one, as has been done in Chapter 2. On the other
hand, it should be noted that the major uncertainty will come from the models themselves
in the sense that their input physics is an approximation of reality. If, e.g., there would be a
stratified ξ, then the results might be rather different.
For this reason, we avoided a time-consuming error propagation. In practice, we have adopted
the following ‘typical’ errors throughout the remainder of this investigation:
– ∆Teff : The errors on the effective temperature will be different above and below 20,000 K,
due to the different grid step, which is 1,000 K for temperatures above 20,000 K and 500 K
below 20,000 K. Since the derived Teff will to some extent depend on the accuracy of the
determined EW, we consider the grid steps as a good estimate for the errors on Teff .
– ∆log g: Even though half the grid step would give us an accuracy of 0.05 dex, we will,
also for log g consider the grid step of 0.10 dex as the typical error, in view of the strong
coupling with the effective temperature.
– ∆ξ: The grid steps can be considered as representative for the errors on ξ. Their values
depend somewhat on the actual derived values, since the grid microturbulences are not
equidistant. The errors will not exceed 6 km/s.
– ∆log n(Si)/n(H): In view of the crude grid steps in log n(Si)/n(H), we consider half the
grid step, i.e. 0.15, as a reasonable estimate for the error on the derived Si abundance.
– ∆n(He)/n(H): In analogy to the Si abundance, also for the He abundance, half the grid
step, i.e., 0.025 will be taken as a conservative error estimate.
4.6 Alternative strategy for a grid-based method
143
– ∆log Q: Since we restrict to the grid points to find log Q, without further refining, the
absolute errors are about half the grid steps, i.e. generally about 0.10 dex, at least for
strong winds. For weak winds (M˙ . 10−7 M⊙ /yr or more or less for log Q . −13.80),
the errors soon become significantly larger. Indeed, in case of low mass loss rates, Hα is
nearly ‘photospheric’ and (almost) only the core will be affected.
– ∆β: Likewise the errors in log Q, the error in β is set by approximately half the grid step,
i.e. of the order of 0.20 dex for the three lowest values (0.9, 1.2 and 1.5) and of the order
of 0.5 dex for the larger values (2.0 and 3.0).
– ∆vmacro and ∆v sin i : The errors on the macroturbulent and rotational profiles are derived from the standard deviation of the values resulting from different line profiles, and
will be different for every object.
4.6 Alternative strategy for a grid-based method
The method described in Section 4.4 does not make full use of the grid, since it will only
check a few tens or hundreds well chosen models within the grid of almost 265,000 models.
The following questions may thus arise: “Is the optimal fit we find with the above described
method, really a global minimum? Or is it just a local optimum and are there other solutions
which might fit equally well, or even better, in another subspace of the grid?”. To answer
these questions, we have additionally considered an alternative approach for an automated
procedure for the line profile fitting procedure.
4.6.1 Description of the alternative method
Instead of working with a method which iteratively updates the fundamental parameters by
selecting only the best values at each iteration, another way to proceed would be the following. One could ‘simply’ compare the EWs of all models in the grid and the accompanying
line profiles with the observed ones. To select then, from this grid, only the best fitting models, we need to introduce a criterion which gives the goodness-of-fit of each model. This
criterion should help to decide which model is most likely the best representation of the true
stellar model.
For this purpose, we have introduced the following loglikelihood function, in analogy with
the one introduced in Eq. (4.1) for the EW:
2 #
nj "
NX
lines X
√
1 yj,i − µj,i
−ln(σj ) − ln( 2 π) −
l≡
(4.5)
2
σj
j=1 i=1
with Nlines the total number of (H, He or Si) lines used, nj the considered number of wavelength points within the line profile, yj,i the observed quantity at point i of line j, µj,i the
144
GAUDI B star sample
Table 4.5: (1-α) 100% quantile of the χ2p -distribution, with p degrees of freedom and a
significance level α = 0.05, taken from Decin et al. (2007). We have repeated here only the
values for up to eight free parameters.
p
1
2
3
4
5
6
7
8
χ2p
3.8414588
5.9914645
7.8147279
9.4877290
11.070498
12.591587
14.067140
15.507313
expected (i.e. theoretically predicted) value for this quantity, and σj the error on the observed
quantity. For the basic statistics behind this formula, and the assumptions made, we refer to
Decin et al. (2007). The maximum of this loglikelihood distribution will give us the model
which fits the observations best. However, models that are sufficiently close to this maximum
loglikelihood can still be accepted as good models. Similar to Eq. (4.2), sufficiently close is
quantified as
2 (lmax − lmodel) ≤ χ2p ,
(4.6)
where lmodel is the evaluation of the loglikelihood function for the model, or the free parameters, under consideration, p is the number of free parameters that one wants to fix by
applying this criterion, and χ2p determines the 95% quantile of the chi-square distribution
with p degrees of freedom as listed in Table 4.5 and taken from Decin et al. (2007).
We will use Eq. (4.5) with the corresponding criterion in Eq. (4.6) in two different ways.
First, we will use it to fit the equivalent widths of the lines, to determine, e.g., the effective
temperature. In this case, we use Eq. (4.1) and (4.2). We have also developed an alternative
likelihood function by using line-ratios instead of absolute equivalent widths (with an automated detection of which ionisation stages can be used). Second, we will use it to fit the line
profile shapes. In this case, yj,i represents the flux at wavelength point i in the j-th line profile, and µj,i the flux corresponding to the same wavelength point. Note that, for this reason,
the number and position of the wavelength points should be adjusted to the observations if
necessary, before comparing.
The preparation (e.g. EW determination) remains as it is for the AnalyseBstar approach.
Further, for the alternative method, we work in two big steps:
The first step consists of the simultaneous determination of the effective temperature, the
surface gravity, the microturbulent velocity and the abundances of silicon as well as helium,
by exploitation of our knowledge of the EWs of both ions. We have decided to include both
He and Si lines at the same time. In this way we hope to be independent of the temperature
domain (i.e. early versus late spectral type). For the late type stars the He lines will have
a higher weight where the Si lines almost all disappear (except for Si II), as helium reacts
4.6 Alternative strategy for a grid-based method
145
sensitively to changes in the effective temperature. On the other hand, in the domain where
the He lines become independent from Teff , the Si lines will be given a higher weight. We still
make the additional assumption in the cool regime that the Si abundance is known. Therefore,
we should obtain for all three abundance grids solutions with different Teff /ξ, at comparable
likelihoods. The surface gravity log g has also been included in this first iteration, since it
also has a certain influence on the strength of the Si lines, and on the He lines in the hot
temperature regime. In this way. we will fix in first instance 5 parameters from the equivalent
widths. Therefore the models should satisfy the following criterion:
2 (lmax − lmodel ) ≤ 11.070498.
To be able to apply the loglikelihood formula for the EWs and this criterion, we first need
the synthetic EWs. For this purpose we have created extremely large (IDL) structures which
contain the theoretically predicted EWsyn,j for all lines of all models. Once we have this
information, it is straightforward to compare the observed EWs with all theoretical EWs in
the grid.
After a first selection of models that agree with the observed EW (which was done irrespective
the wind parameters), in a second step, we aim at selecting the best line profile fits, not only
for the Balmer lines, which give us information about the three parameters log g, β and log Q,
but also for the He and Si lines, which give a second constraint for the remaining parameters.
First, we derive the best vmacro (in the same way as described in Section 4.4.4) for each of
the models selected in the first step. Then, we use the Balmer lines to constrain log g and the
wind parameters.
Whereas the two wind parameters are determined from only one line (i.e. Hα) and the gravity
from only two or three lines (depending on the wind strength), the He and Si lines are much
more numerous. Therefore, we approximated the overall loglikelihood by weighting the
separate loglikelihoods with the number of lines, N , used to derive the parameters:
l = llogg .
Nlines
Nlines
+ lwind . Nlines + lSi,He .
,
NH
NSi + NHe
where NH is the amount of hydrogen lines used to derive log g (usually 3), NSi and NHe
are the amount of He and Si lines used to derive Teff , vmicro and the abundances of Si and
He (for the second constraint from the line profiles, and not from the equivalent width), and
Nlines is the total amount of lines used (= NH + NHe + NSi ). llogg , lwind and lSi,He are the
corresponding loglikelihoods, calculated respectively from the Balmer lines for log g, from
Hα for the wind and from all He and Si lines for the remaining parameters. This approach
is motivated by two facts: 1) the l values are additive (by construction), and 2) since the
individual l values are a sum over the considered lines, by dividing by the number of lines
used, we construct in principle for each process an l per line, and add those. Multiplication by
Nlines does not change the shape of the likelihood function and is done as final normalisation.
By assigning these weights, the loglikelihoods will get lower weight, when more lines are
used, and higher weight when only few diagnostic lines are used. On the other hand, all
wavelength points have been given equal weight. In a more sophisticated implementation,
one could think about giving higher weight to, e.g., the forbidden lines in the blue wings of
He I 4471 and He I 4922.
146
GAUDI B star sample
The criterion for the models to be selected as well-fitting models is then
2 (lmax − lmodel ) ≤ 15.507313,
where lmax is now the maximum of the combined loglikelihood function of all line profiles,
as mentioned above. In this way we are, in most cases, left with only one solution. It would
obviously be very interesting to compare this unique solution to the one we found for some
of the stars with AnalyseBstar. In the next section, we show the results of such a comparison
for the B2.5 Ib/II supergiant ι CMa.
For the sake of clarity, in the remainder, we will refer to this alternative method as the full
grid search or loglikelihood method (shortly, the l-method). We stress that we have only
implemented a basic version of this method, and it should be further finetuned if one really
wants to rely on the outcoming results. We have refrained from spending time on this, because
it soon became clear that the method has a few disadvantages which make it an insufficient
method for our purposes. We refer to Section 4.6.3 for a discussion of the usability and
reliability of this method.
4.6.2 Comparison with AnalyseBstar for ι CMa
For ι CMa we applied both AnalyseBstar and the full grid search, coupled to the likelihood
functions in Eqs (4.1) and (4.5), to determine the physical parameters. We compare the outcome in Table 4.6. Most parameters agree fairly well, except for the Si abundance. This could
be the result of forcing an optimal fit to one of the grid points when applying the l-method,
whereas in AnalyseBstar we leave the possibility for the Si abundance to lie in between two
grid points. The value for the Si-abundance, derived from our AnalyseBstar analysis, is
log n(Si)/n(H) = −4.78, which makes ‘-4.79’ the closest grid value, whereas the l-method
indicates log n(Si)/n(H) = −4.49 as the best Si abundance.
In Fig. 4.25, we compare the Si line profiles of both models (we verified that the difference in
H and He is minimal). For the Si III lines, there is only little difference between both models.
However, it appears that we are dealing with the same problem as in Chapter 2, namely that
the different components of the Si II doublets cannot be equally well fitted. Where AnalyseBstar chose Si II 4130 and 5056, the l-method seemed to have given preference to Si II 4128.
The Si II 5041 line lies in the middle between the two grid models, as would be expected
from the interpolated Si abundance.
In the first step of the l-method, 7652 models throughout the whole grid have been selected
as well-fitting models on the basis of the equivalent widths of He and Si. Then, finally, all
synthetic line profiles of these almost 7700 models had to be compared with the observed line
profiles. It took about three hours to finally decide on the best model, whereas AnalyseBstar
only needed half of the time to get the same results. Only 66 different possibilities had to be
checked throughout the procedure.
Note that we find an extremely low He abundance. This could be the result of an extrapolation
147
4.6 Alternative strategy for a grid-based method
Table 4.6: Comparison between the resulting parameters for ι CMa obtained by, on the one
hand, the automatic tool ’AnalyseBstar’ (Section 4.4) and, on the other hand, the alternative
method which makes use of a loglikelihood function, and a full grid search.
fit parameters
Teff (K)
log g (cgs)
log Q (char)
β
n(He)/n(H)
log Si/H
ξ (km/s)
vmacro (km/s)
interpolated
(AnalyseBstar)
17200
0.05
-4.78
16
39 ± 16
best fit
(AnalyseBstar)
17000
2.6
b
1.2
0.10
-4.79
15
39 ± 16
loglikelihood
method
16500
2.6
b
1.5
0.10
-4.49
10
34
Figure 4.24: Contour plot of the best fitting models for ι CMa through the application of
loglikelihood functions. In the first step, the grid is constrained to a subspace of the grid
by comparing the equivalent widths of He and Si lines (black dots). In a second step, we
compare the line profile shapes and we are left with only one model, the one indicated in
white. The rainbow colours represent the third dimension of the loglikelihood distribution
over the subgrid (low l (black) to high l (red)).
towards abundances lower than solar, but, on the other hand, Morel (2007) obtained a similar,
suspiciously low He abundance of 0.045 ± 0.02, which may indicate that it might be real.
148
GAUDI B star sample
Figure 4.25: Comparison of the two different models chosen as the best fitting model for the
observed Si lines for ι CMa (grey) by the l-method (full) and AnalyseBstar (dashed).
4.6.3 Discussion points
Several points can be raised against the actual application of this method. In this section we
address a few points of discussion, which should be kept in mind when applying either of the
two methods.
1. A first important point, which becomes especially of crucial importance when aiming at
the analysis of large samples, is the speed of performance. Comparison of equivalent
widths is straightforward and can be done within a minimum of time. However, the large
number of line profiles that need to be checked, slow down the l-method considerably. For
high-quality data with a high SNR (and thus small errors), the l-method requires at least
one hour of computation time to select the best model, rising exponentially with the error
4.6 Alternative strategy for a grid-based method
149
on the equivalent widths. Indeed, for low SNR spectra, where the errors on the equivalent
widths are much larger, the amount of models that were accepted in the first step is very
large, and, therefore, it can take up to more than 20 hours8 before the final solution is
found. Moreover, in most cases, this final solution is very similar to the one found by
AnalyseBstar in less than a quarter of an hour. For AnalyseBstar, the time needed to
converge to the optimal solution will not only depend on the quality of the data, but also
on the temperature region and on the starting values. The closer the initial values are to
the real values, the faster convergence is reached (typically within a few minutes). When
more than one Si ionisation stage is available (generally this is no problem for effective
temperatures above 15,000 K), the final solution is, on average, found within 15 minutes.
When only one ionisation stage of Si is available, also AnalyseBstar is slower. This can
easily be understood, when one remembers that, in the cool regime, we use a procedure
based on the same principles as this alternative method (i.e. using loglikelihood functions).
Therefore, for cool stars, the computation time can go up to one hour and a half, which is
still much faster than the alternative method.
2. A second disadvantage of the alternative method, is the fact that it is bound by the grid.
Whereas our method is able to search in between the grid points (thanks to the parallel
determination of Teff , ξ and abundance), the alternative method searches the best solution
only among the grid models. This means that it forces the observed line profiles to be
exactly represented by the best model. We know that, e.g., with respect to Si abundance,
the grid is too coarse. A lot of B stars seem to display a depletion in Si (cf. Section 4.5.2),
which means that, likely, almost all abundances lie somewhere between -4.49 and -4.79
(or even lower). The alternative method will force the abundance to either of the two
edges, and, consequently, also the temperature and microturbulence will be affected. For
the microturbulence, a same reasoning can be made, although the effect might be a bit less
worrisome because of the smaller grid steps.
3. The l-method assumes that, apart from noise, all observed profiles can be represented
perfectly by the model. However, as far as FASTWIND (or any alternative code) concerns, we know that this is not the case and that there may still be systematic deviations
between observed and theoretical profiles (see, e.g., the discrepancy between the two lines
of Si II 4128-4130, as discussed in Section 2.7). As long as these systematic deviations are
rather similar over the whole grid, their presence might not influence the likelihood ratios
too much. In our case, however, deviations between observed and theoretically predicted
profiles seem to be rather different across the B-type spectral range. Though different,
the deviations are typically small.If we could quantify this ”theoretical” bias as a function
of stellar parameters, we could account for this by introducing carefully chosen weights.
This is exactly what the eye is doing in a fit-by-eye procedure (and, therefore, we were
able to correct for this discrepancy in the analysis of the supergiants in Chapter 2). However, it is not evident to translate this into a code. For O-type stars, the reliability of the
individual lines as a function of parameters is well known, thus allowing Mokiem et al.
(2005) to account for it in their genetic algorithm approach. However, for B-type stars,
8 Even though the efficiency might be improved upon, due to the large number of models to be checked the
method will still remain sufficiently slower than the iterative method (see, e.g., the ana- lysis of the late B stars,
which is partly based on the same approach and already takes sufficiently longer computation time), and therefore
(while keeping the other disadvantages in mind), we have not invested more time in it.
150
GAUDI B star sample
this knowledge can only be built up through the analysis of large samples covering a wide
variety of B-type stars. Hopefully we can provide part of this information after the analysis of the large CoRoT sample. Since almost all GAUDI stars are main sequence stars,
we will only be able to give more information for the dwarfs. Of course this problem of
imperfect line predictions is not inherently connected to the l-method. Also our present
approach (AnalyseBstar) suffers from it.
4. One advantage of the l-method is the fact that it is able to point to all relevant optima
in the possible parameter space, and, in this way, it avoids to end up in a local optimum
instead of the global optimum. However, also with our iterative method, we are sure to
end up in the correct domain of the HRD, since there should be only one maximum, due
to the monotonic behaviour of the ionisation fractions.
For both methods the accurate measurement of the observed EWs is of major importance,
since these will fully determine which models are accepted or rejected.
4.7 Results for the GAUDI stars
In Table 4.7 we list the results of the application of AnalyseBstar to a selected sample of
the B-type stars in the GAUDI database9 . As mentioned in Section 4.2 we have immediately
excluded the known spectroscopic binaries and the Be stars, since our methods are not appropriate to handle these stars, and they have been analysed by other members of the CoRoT
team. The quality of both the SARG and CATANIA spectra was insufficient to perform an
accurate detailed analysis, so these have not been considered either.
For HD 174069, we had both an ELODIE and a FEROS spectrum available. From inspection
of the spectra, it was immediately clear that there are systematic differences in the line profiles, in particular the wings of the Balmer lines are much less pronounced in the ELODIE
spectra (see Fig. 4.26). It is not clear what causes this discrepancy. Also other stars, for which
we had both a FEROS and an ELODIE spectrum available, showed the same discrepancy, so
it seems to concern a systematic, maybe instrumental, effect. We fitted both spectra and came
up with a different set of parameters (see Table 4.7). The difference in log g is large and certainly worrisome. When accounting for the difference in Teff , this is a discrepancy of at least
0.6 dex. Our first thought was that the discrepancy was caused by the shorter wavelength
range of each order in the ELODIE spectra, which affects the normalisation. We have downscaled the FEROS normalisation to exactly the same wavelength range as for the ELODIE
spectra. The problem remained, however, and we did not succeed in getting a hold on it. Due
to the fact that the merging of the ELODIE spectra is in any case less reliable than the one
for the FEROS spectra (cf. Section 4.4.1), we have decided to use the FEROS spectra, and to
leave the ELODIE spectra out of our sample, in order to keep our study as homogeneous as
possible. Of course, this significantly reduces the extent of our sample. With our final aims
9 Based on GAUDI, the data archive and access system of the ground-based asteroseismology programme of the
CoRoT mission. The GAUDI system is maintained at LAEFF, which is part of the Space Science Division of INTA.
151
4.7 Results for the GAUDI stars
Figure 4.26: HD 174069: Significant (and apparently systematic) discrepancies are observed
between the line profiles of the FEROS (black) and the ELODIE (grey) spectrum. This leads
to considerable discrepancies in the derived stellar parameters.
in mind, namely to derive accurate fundamental parameter calibrations, we preferred to base
our analysis on thrustworthy stars, with qualitative spectra.
For the same reason, we have also excluded stars with an very bad SNR and the fast rotators,
because an accurate equivalent width determination is impossible in these cases. In this way,
we were left with 35 reliable targets. In Table 4.7, we give an overview of their derived
properties. Their position in the HRD can be found in Fig. 4.27.
To enlarge the number of stars in the hot temperature side of the spectral range and to allow
for a continuous flow towards the O-type regime, R. Mokiem has been so kind to additionally
analyse 2 of the 11 O-type stars in the GAUDI stars, using his genetic algorithm approach
(Mokiem et al. 2005). Since this approach is very accurate and is based on the same atmosphere code (i.e. FASTWIND), they can be added to our sample, without any objection. We
list the resulting parameters in Table 4.8.
4.7.1 Comparison between observed and synthetic equivalent widths
For all the analysed GAUDI stars, we have compared the observed equivalent widths with
the theoretically predicted equivalent widths of the best fitting model, in order to see if there
exist any systematic differences. This simple exercise provides very useful information, to
improve the theoretical line prediction or to confirm the currently existing predictions. We
have investigated the distribution of
EWobs − EWtheo
EWtheo
152
GAUDI B star sample
Table 4.7: Basic fundamental parameters of a selected sample of GAUDI stars, resulting from the
application of AnalyseBstar. ‘Interpolated values’ as well as the derived closest grid points are given.
Spectral types are taken from SIMBAD. For a summary of the errors, we refer to Section 4.5.3. The
errors on v sin i and vmacro are the standard deviations of the values derived from different Si lines. If
only one line was used for deriving vmacro , we fixed the error at 5 km/s. For most of the latest spectral
types (up to B5), there is no Si III anymore and we had to use a different approach, adopting a value for
the Si abundance. In this case the derivation of the interpolated Teff was impossible, and we could only
derive the interpolated microturbulence.
HD SpT
48434
172488
52382
174069F
174069E
45911
45546
170580
44700
45726
51507
181074
43157
48215
58973
177880
178744
48807
51360
42677
44948
47964
48497
49935
51150
53083
44321
44354
48808
48957
53004
54761
173370
181440
182198
B0III
B0.5V
B1I
B1.5V
B1.5V
B2IV-V
B2V
B2V
B3V
B3
B3V
B3
B5V
B5V
B5
B5V
B5
B7Iab
B7III
B8
B8Vp
B8III
B8
B8
B8
B8
B9
B9
B9
B9
B9
B9
B9V
B9III
B9
interpolated
Teff
(K)
grid
Teff
(K)
log geff
(cgs)
28100
22200
23100
22100
19900
20800
20100
19700
17100
21000
17500
19800
16500
15000
15000
15000
14500
13000
13800
10000
11500
12000
14000
12500
13500
11500
11500
13000
12000
12000
11000
11500
11000
11000
11000
28000
22000
23000
22000
20000
21000
20000
19500
17000
21000
17500
20000
16500
15000
15000
15000
14500
13000
14000
10000
11500
12000
14000
12500
13500
11500
11500
13000
12000
12000
11000
11500
11000
11000
11000
3.10
3.10
2.70
4.20
3.40
4.00
3.80
4.00
3.90
4.00
3.50
3.60
3.60
3.90
3.30
3.90
3.50
2.10
3.30
3.10
3.50
3.10
3.90
3.30
2.80
2.90
3.80
3.90
3.20
3.10
3.90
3.10
3.20
3.50
3.30
interpolated
He
Si
0.10
0.10
0.15
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.20
0.10
0.10
0.10
0.10
0.10
0.10
0.10
0.15
0.10
0.10
0.10
0.10
0.10
0.10
0.10
-4.50
-4.36
-4.71
-4.50
-4.81
-4.81
-4.81
-4.86
-4.62
-4.44
-4.51
-4.35
-4.79
-4.49
-4.79
-4.79
-4.49
-4.55
-4.04
-4.79
-4.79
-4.79
-4.49
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
grid
Si
-4.49
-4.49
-4.79
-4.49
-4.79
-4.79
-4.79
-4.79
-4.49
-4.49
-4.49
-4.49
-4.79
-4.49
-4.79
-4.79
-4.49
-4.49
-4.19
-4.79
-4.79
-4.79
-4.49
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
-4.79
interpolated grid
ξ
ξ
(km/s) (km/s)
17.40
13.00
19
3.80
8.50
6.60
4.50
5.90
1.40
5.10
4.40
12.6
1.25
3.94
7.38
9.57
11.80
9.70
2.00
7.59
16.30
3.59
1.63
5.78
4.38
3.39
7.10
9.75
4.00
1.63
7.90
4.81
5.84
2.75
6.83
15
12
20
6
10
6
6
6
3
6
3
12
3
3
6
10
12
10
3
6
15
3
3
6
3
3
6
10
3
3
6
6
6
6
3
v sin i
(km/s)
vmacro
(km/s)
63 ± 5
101 ± 7
61 ± 6
35 ± 5
45 ± 4
23 ± 4
66 ± 10
9±1
12 ± 2
113 ± 10
147 ± 6
133 ± 8
26 ± 3
103 ± 1
95 ± 4
45 ± 4
236 ± 32
25 ± 2
68 ± 8
118 ± 17
74 ± 10
50 ± 4
9±2
71 ± 4
12 ± 2
76 ± 20
98 ± 8
121 ± 0
47 ± 15
24 ± 3
51 ± 7
54 ± 1
262 ± 53
53 ± 7
25 ± 1
49.2 ± 5.0
80.0 ± 9.2
61.7 ± 10.3
34.9 ± 7.8
37.6 ± 5.0
0.0 ± 5.0
57.6 ± 3.6
12.7 ± 5.0
0.0 ± 5.0
55.6 ± 1.0
0.0 ± 5.0
0.0 ± 2.8
20.4 ± 3.2
100.2 ± 5.0
82.9 ± 5.0
52.1 ± 1.9
0.0 ± 5.0
26.0 ± 5.0
53.8 ± 9.5
81.2 ± 5.0
106.9 ± 0.0
52.8 ± 29.1
18.0 ± 5.0
77.4 ± 5.0
39.6 ± 5.2
90.0 ± 5.0
0.0 ± 24.6
72.9 ± 5.0
63.3 ± 10.5
23.1 ± 3.9
54.2 ± 9.4
71.8 ± 5.0
0.0 ± 5.0
0.0 ± 5.0
0.0 ± 5.0
as a function of the effective temperature, surface gravity and microturbulent velocity (Figs
4.28 till 4.30). These diagrams not only teach us about systematic differences between the
observed and synthetic equivalent widths, but also tell us which lines are observable in which
regions, and which lines are reliable to use as indicator for Teff , log g or ξ. For each line,
we have also indicated the linear least squares fit. Positive values in the diagram mean that
the observed EWs are larger than the theoretically predicted ones, negative values mean that
153
4.7 Results for the GAUDI stars
Figure 4.27: Position of the GAUDI stars in the HRD, based on the detailed spectroscopic
analysis performed with AnalyseBstar.
Table 4.8: Fundamental parameters of two GAUDI O-type stars, as derived by R. Mokiem,
using the genetic algorithm of Mokiem et al. (2005), and FASTWIND model atmospheres.
Teff (K)
log g (cm/s2 )
log gc (cm/s2 )
n(He)/n(H)
R∗ /R⊙
L∗ /L⊙
ξ (km/s)
M˙ (M⊙ /yr)
β
v∞ (km/s)
v sin i (km/s)
M∗ /M⊙
HD 47432 (O9.5 III)
29300
3.14
3.18
0.10
17
1.9e+05
18
6.975e-07
1.64
1600
126
16
HD 52266 (O9 V)
31400
3.64
3.73
0.10
11
1.0e+05
19
6.864e-07
0.66
2100
274
22
FASTWIND predicts too large EWs. As a first overall conclusion, we can say that, generally,
the lines are all very well predicted with only a few exceptions. Note that we have not used
the three He I lines, He I 4010, 4120 and 4140, in our spectral line fitting procedure, because
we knew on beforehand that the broadening functions are very approximative and need to
be updated. We see that all three of them are systematically underestimated. The He I lines
can be observed all over the B-type spectral range, whereas we only have one star for which
He II 4541 could be measured, and only two stars for which He II 4686 is detectable. The
Si II lines are visible up to 22,000 K, the Si III triplet at 4552 from 13,000 up to 28,000 K,
whereas the Si III 4716 and the Si III 4813 triplet are both only visible in a very limited range
154
GAUDI B star sample
Figure 4.28: Behaviour of the relative differences between the observed EW and the theoretically predicted EW of the best fitting model EWobs /EWtheo − 1 as a function of the
effective temperature, surface gravity and microturbulent velocity, for eight He I lines.
4.7 Results for the GAUDI stars
Figure 4.29: Same as for Fig. 4.28, but for different lines.
155
156
GAUDI B star sample
Figure 4.30: Same as for Fig. 4.28, but for Si lines.
4.8 Photometric calibrations
157
between 18,000 and 23,000 K. Si IV is detectable only in three cases, which all seem to be
supergiants with large microturbulent velocities. Most other stars have low microturbulent
velocities, as was expected.
We have to be very careful in drawing additional conclusions from these diagrams. Differences between observational and synthetic EWs are not necessarily due to errors in the
theoretical predictions, but can also be due to a misfit of the Gaussian profile to the observed
line, in measuring the observed EW. As mentioned before, for lines such as He I 4471 and
He I 4922, it is difficult to get an accurate value for the EW because of the forbidden component in the left wing and the non-Gaussian distribution as well as due to the many blends in
the right wings. We have looked into the line profile fits to see if there are indeed systematic
differences, or if this is just an artificial effect due to an inaccurate EW measurement. As a
means of illustration, we display in Fig. 4.31 until 4.36, a few examples for different types of
stars: a cool/hot dwarf, a cool/hot giant and a hot supergiant (the only one in our sample).
They show the quality of the GAUDI spectra and the difficulties with which we have to deal.
The figures are discussed in their captions.
4.8 Photometric calibrations
Direct means to determine the physical parameters of stars are available, such as, e.g., cluster main sequence fitting, parallax measurement and analysis of binary motion to determine
respectively the age, luminosity, and mass and radius of a star. The unfortunate thing is that
these methods can only be applied to a small number of targets, and they rapidly loose power
when going to faint objects. One way to overcome this problem, apart from the spectroscopic
analysis, is brought to us through multicolour photometric measurements. Photometric colour
indices reflect, by definition, the difference in emitted energy at different wavelength ranges.
Therefore, they are sensitive to stellar properties, such as gravity and effective temperature,
which alter the energy distribution over the spectral range. To derive information about the
star (e.g., its effective temperature), one has to find the most suited colour indices of a given
photometric system, i.e. those colour indices which are most sensitive to changes in this parameter.
For early type stars, the depth of the Balmer discontinuity is very sensitive to temperature
changes, and, therefore, a good colour index is one that combines filters on either side of the
discontinuity. The gravity, on the other hand, can be determined from the equivalent width
of the Hβ line (in case there is no wind contamination). Thus, a good colour index is one
that combines a broad and a narrow filter, both centred on the Hβ rest wavelength. The broad
filter measures the EW of the line and the narrow filter the line depth, which, combined, give
the FWHM of the line. Once we dispose of the right colour indices, the ‘only’ thing required to derive the fundamental parameters, is an accurate calibration of these parameters as
a function of the observed colour indices. Though this may sound easy, much effort has been
done in order to derive such calibrations. Yet, for the hottest stars, it seems hard to establish. The main reason for this is that, at these hot temperatures, we are on the Rayleigh-Jeans
tail throughout the optical spectral range. Therefore, optical colours are rather insensitive to
changes in the effective temperature for O and early B-type stars.
158
GAUDI B star sample
Figure 4.31: Example of the line profile fits for a cool dwarf : HD 44321. The preparation
was not easy for this spectrum: the EW determination of He I 4387 and He I 4922 is hampered
by strong blends, the normalisation for Si II 4128/4130 is uncertain due to their position in
the broad Stark wings of Hδ and the normalisation of He I 4026 seems to bit a bit offset in
the blue wing. At this temperature, there is also no Si III left anymore, so we need to make
an assumption for the Si abundance. Despite these difficulties, we were still able to obtain a
satisfying fit quality. From the Si II profiles it is immediately clear that this is a fast rotator,
where the macroturbulence is absent. Omitted lines were not present in the spectrum.
4.8 Photometric calibrations
159
Figure 4.32: Example of the line profile fits for a cool giant: HD 181440. Despite the low
quality of the data, and the fact that we have no more Si III left in this temperature region, we
are still able to obtain a satisfying fit quality. Omitted lines were not present in the spectrum.
160
GAUDI B star sample
Figure 4.33: Example of the line profile fits for the only cool supergiant in our sample:
HD 48807. The Balmer lines show some bumps in the wings. It is unclear whether these
bumps are real or not. They also appear in the other supergiant in the GAUDI sample. The
He and Si abundances are both lower than solar, though lying closest to the solar grid point.
This clarifies the small discrepancies in the cores of He as well as Si. Note the excellent fit of
the He I 4471 and 4922 forbidden components.
4.8 Photometric calibrations
161
Figure 4.34: Example of the line profile fits for a hot dwarf : HD 44700. The strong lines
on the blue side of Si III 4716 and Si III 4819 are easily mistaken for the Si lines themselves.
This is the reason why we always indicate the exact position of the transition during the
full preparation process. The misfit of He I 4010 and 4140 is due to incomplete broadening
functions of these lines. Remember that these lines were not used during the fitting procedure,
but only as a double check afterwards. We observe a similar behaviour in several other stars.
All other lines, even the weakest, fit very well. The shape of He I 4471 is even perfect.
162
GAUDI B star sample
Figure 4.35: Example of the line profile fits to the ELODIE spectrum for a hot giant:
HD 48434. Both the He and Si abundance are a bit lower than solar, but for the remainder, this model gives a very good fit quality.
4.8 Photometric calibrations
163
Figure 4.36: Example of the line profile fits for the only hot supergiant in our sample: HD 52382. The peak of Hα is reasonably well represented, despite the fact that
the fitting of the wind is restricted to the grid combinations of log Q and β. The Balmer
lines show bumps in the wings. It is not clear whether this is real and due to the strong
wind, or if this is a spectral artefact. It also arises in the cool supergiant. The interpolated Si abundance (log n(Si)/n(H) = −4.71) is higher than the closest grid abundance
(log n(Si)/n(H) = −4.79), which clarifies the discrepancies in the Si line cores.
164
GAUDI B star sample
We aim at verifying and improving (if necessary) the commonly used Teff and log g calibrations for both the Geneva and the Str¨omgren photometric system. We attempt to do this
by providing accurate, fundamental parameter estimates (as derived from the application of
AnalyseBstar to high-resolution spectroscopic data) for targets for which also photometric
colour indices are available. These will then be used to check and/or to recalibrate both
photometric systems. The GAUDI sample is ideally suited for this job. As mentioned in
Section 4.1.2, in the framework of CoRoT, an intensive ground-based observing programme
was set up to gather Str¨omgren colour indices of more than 1500 stars, among which some
320 B-type stars. For 250 of these targets, also high-resolution spectroscopy was available.
However, 60 targets were excluded, as they are spectroscopic binaries or Be stars. On the
other hand, we have 183 B-type stars for which we have Geneva photometry at our disposal,
and for 153 among them we have also spectroscopic information. At present, 25 of these
targets have been spectroscopically analysed. In what follows, we will mainly concentrate on
the Str¨omgren photometric system, but we still thought it useful to discuss also the Geneva
photometric system, as it allows a comparison between the physical parameters derived from
the fundamental parameter calibrations of both systems.
4.8.1 Geneva photometric system
The Geneva system (Golay 1974, 1980) consists of seven filters (U, B, V, B1, B2, V1 and
G). In this multicolour parameter space, Cramer & Maeder (1979) defined an orthogonal
coordinate system (X, Y, Z):
X = 1.3764 [U − B] − 1.2162 [B1 − B] − 0.8498 [B2 − B]
−0.1554 [V 1 − B] + 0.8450 [G − B] + 0.3788,
Y = 0.3235 [U − B] − 2.3228 [B1 − B] + 2.3363 [B2 − B]
+0.7495 [V 1 − B] − 1.0865 [G − B] − 0.8288,
Z = 0.0255 [U − B] − 0.1740 [B1 − B] + 0.4696 [B2 − B]
−1.1205 [V 1 − B] + 0.7994 [G − B] − 0.4572.
This (X, Y , Z)-representation has brought considerable improvement in the separation of
the effects of gravity and temperature, compared to the original colours. Indeed, in the HRD,
the X-axis is oriented in the direction of the main sequence O and B stars, and is therefore
a good temperature indicator (for log g ≥ 3.5). The Y -axis is perpendicular to X, and lies
in the direction of the luminosity variations. Therefore it is a good gravity indicator. Z is
orthogonal to the (X, Y ) plane, and is a measure of stellar peculiarities (North & Nicolet
1990).
Based on this (X, Y , Z)-representation, Cramer & Maeder (1979) have defined an effective
temperature calibration for dwarfs and giants, based on the temperature estimates of a number
4.8 Photometric calibrations
165
of stars by Code et al. (1976), which were derived from measurements of the angular diameter
and the total absolute flux, integrated over the entire spectrum10:
log Teff = 4.496 − 0.453X + 0.086X 2.
At that time, three stars had not been sufficiently well measured in the Geneva system and
the relation performed bad for Teff ≥ 25, 000 K. Cramer (1984b) improved this effective
temperature calibration by using the latest Geneva colour indices for these stars. Moreover,
he excluded the peculiar stars in the list, and corrected the colours of the binaries to the values
of the primary component. In such way, he came up with the following calibration for O-, Band early A-type stars:
log Teff = 4.586 − 1.038X + 1.094X 2 − 0.646X 3 + 0.139X 4.
However, this calibration did not take into account the slight dependence of X on log g.
Therefore, North & Nicolet (1990) proposed a new calibration for B-type stars, allowing
to predict accurate Teff and log g values, for the measured X and Y parameters, for main
sequence stars with Teff ≥ 10, 000 K, based on Kurucz (1979) atmosphere models and a set
of ‘standard’ stars of Code et al. (1976) and Lanz (1987).
With improved Kurucz models (Kurucz 1993), also this calibration needed revision. This was
done by Kunzli et al. (1997). Instead of comparing the observed colour indices with those interpolated in the “direct” grids of synthetic colours from the known fundamental parameters
(as e.g. Lester et al. (1986) for Str¨omgren), they preferred to compare the fundamental physical parameters with those interpolated in the inverted grid from the observed colours. With
the aim to obtain the position of the CoRoT target stars in the HRD, C. Aerts interpolated in
the (X, Y ) grids of Kunzli et al. (1997), to obtain the fundamental parameters Teff and log g
of all stars in the eyes of CoRoT for which Geneva photometry is available. Fig. 4.37 shows
the derived Teff and log g estimates.
4.8.2 Str¨omgren photometric system
The Str¨omgren or ubvyβ system, is the photometric system for which the most attempts have
been made to establish an accurate calibration in terms of effective temperature and gravity
for early type stars. Generally the respective colour indices c0 and β are used for this. The
index c0 is the dereddened version of c1 = (u - v) - (v - b), which is a measure of the depth
of the Balmer discontinuity, and thus of the effective temperature. As mentioned above, the
equivalent width of Hβ (and therefore the Str¨omgren parameter β) is a measure of the gravity.
Balona (1984) used Kurucz (1979) model atmospheres to establish a theoretical functional
relationship between effective temperature and gravity on the one hand, and the Str¨omgren
indices c0 and β on the other hand. He fixed the zero-point of this relationship by using the
effective temperatures of the 32 stars as derived by Code et al. (1976), and found that
log Teff = 3.9036 − 0.4816[c] − 0.5290[β] − 0.1260[c]2 + 0.0924[β][c] − 0.4013[β]2,
10 The angular diameters have been measured with the Sydney University Stellar Interferometer (SUSI) at the
Narrabri Observatory (Brown 1968). The entire spectrum covers the range of UV and visual up to infrared.
166
GAUDI B star sample
Figure 4.37: Effective temperatures and gravities of all stars in the eyes of CoRoT for which
Geneva photometry is available, as a function of the Geneva colours X and Y , as derived
from interpolation in the (X, Y ) grids of Kunzli et al. (1997).
4.8 Photometric calibrations
167
with [c] = log (c0 + 0.20) and [β] = log (β - 2.5). With the update of the theoretical ubvyβ
indices from the published and unpublished grids of Kurucz (1979) by Lester et al. (1986),
and an adjustment of the temperature scale of Code et al. (1976) by Beeckmans (1977) (which
was initially not taken into account, see below), also the calibration of Balona (1984) needed
revision. So, ten years after the original paper, Balona 1994 (hereafter BA94) recalibrated the
effective temperature in terms of Str¨omgren indices, from the latest data available at that time.
He used the theoretical indices of Lester et al. (1986) to determine the functional behaviour
θ = 0.7515 + 0.2915[c] + 0.4875[c]2 + 0.1216[c]3 + 1.0306[β] + 0.7958[β]2 − 0.4392[β][c],
with θ = 5040/Teff , [c] = log (c0 + 0.25) and [β] = log (β - 2.5). This theoretical relation was
then compared to empirical determinations of Teff from Beeckmans (1977) and Malagnini
et al. (1986). Beeckmans (1977) revised the effective temperatures of Code et al. (1976),
by applying a correction to the absolute UV fluxes (resulting from different instrumental
calibrations). Malagnini et al. (1986) applied a different approach and derived the effective
temperature, gravity, bolometric correction and angular diameter in a self-consistent way by
comparing the observed energy distribution to the flux taken from model atmospheres. For
temperatures hotter than 20,000 K, there seemed to be a significant discrepancy between
these ‘observed’ values for Teff and the theoretically predicted ones, which BA94 suggested
to be probably due to a problem with the synthetic colours of Lester et al. (1986). On the
other hand, Malagnini et al. (1986) already pointed out larger uncertainties in the hotter temperature range, due to both the inadequacy of the LTE assumption and the incompleteness of
the adopted blanketing. Balona corrected for this discrepancy by comparing the theoretical
calibration with empirically determined temperatures of Malagnini et al. (1986) to fix the
zero-point. In this way he finally came up with the following improved relation:
θ = 0.7951 + 0.3226[c] + 0.5395[c]2 + 0.1346[c]3 + 1.1406[β] + 0.8808[β]2 − 0.4861[β][c].
(4.7)
Since, for early type main sequence and giant stars, the photometric error in the observed β
index outweighs the change in temperature as a function of luminosity class, he suggested to
use the mean β value for these stars as a function of c0 , as soon as c0 < −0.18 or β < 2.60:
β = 2.620 + 0.2517 c0 − 0.1400 c20 + 0.1704 c30 .
(4.8)
For the surface gravity, Balona followed a similar procedure as for the effective temperature.
First, he used the theoretical ubvyβ indices of Lester et al. (1986) to define the functional
relationship for a wide variety of effective temperatures and luminosities. Subsequently, he
fixed the zero-point using data from double-lined eclipsing binaries (from Andersen 1991),
for which the surface gravities can directly be determined from the masses and radii of both
components:
log g = 7.492 − 5.7432[c] − 1.9660[c]2 − 0.8074[c]3 + 5.9449[β] − 1.6611[β][c]. (4.9)
Using the calibrations in Eq. (4.7) and (4.9), we have calculated the effective temperatures
and gravities for all GAUDI B stars, for which the c0 and β indices were obtained during the
CoRoT preparatory campaigns. For the stars with a spectral type B3 or earlier, and for stars
which have c0 < −0.18 or β < 2.60, we have used Eq. (4.8), following the suggestion by
168
GAUDI B star sample
BA94 (and private communication). Fig. 4.38 shows the derived values for both the effective
temperature (upper panel) and the surface gravity (lower panel) as a function of the diagnostic
colours. The curved trend visible in the lower panel refers to those objects for which Eq. (4.8)
was applied.
4.8.3 Comparison between the two photometric systems
Effective temperatures
We can now compare in how far the effective temperatures derived from the two different
photometric systems agree. In the upper panel of Fig. 4.39, we compare the colour indices,
used to derive the effective temperatures, for the Str¨omgren (c0 ) and Geneva system (X). For
the range X > 0.5 or c0 > 0.2, there exists a tight relation between both colour indices.
Below these values, the relation is rather chaotic. In the lower panel, the resulting effective
temperatures are shown. From this figure, we see that for log Teff < 4.10, the Geneva
temperatures are generally a bit higher than the Str¨omgren temperatures, while from log
Teff = 4.10 on, the opposite is true. We can compare this result with what Heynderickx
(1991) found, since he performed a similar study, but restricted to a sample of β Cephei
stars. For comparison, the ranges of his research are indicated as boxes in Fig. 4.39. He
concluded that the effective temperatures obtained by the Geneva system are systematically
higher than in the Str¨omgren system, which was attributed to a discontinuity in the Str¨omgren
Teff -calibration around 20,000 K (log Teff = 4.3). This is exactly the temperature where the
scatter around the one-to-one relation starts to grow (see Fig. 4.39). Also Heynderickx (1991)
noted that the scatter at these temperatures is rather large, which we confirm.
Surface gravities
Cramer (1984a) has linked the Geneva colours to the Str¨omgren β-index for O, B and early Atype stars, through the definition of a similar β, which is a polynomial of the Geneva colours
X and Y :
β = 2.5909 + 0.0667X + 0.1748X 2 − 0.0612X 3 − 0.6801Y − 2.4676Y 2
+0.4418Y 3 − 0.2559XY + 0.1448XY 2 + 0.2582X 2Y.
Although a linear correlation can be observed between the two photometric β indices, the derived gravities differ appreciably (see Fig. 4.40). For the majority of the stars, the Str¨omgren
gravities are higher than the Geneva gravities. The Str¨omgren gravities of the GAUDI stars
also occupy a wider range in log g (extending from 2.9 to 5.5) whereas the Geneva gravities
lie in a range between 3.2 and 4.5. Heynderickx et al. (1994) found exactly the opposite for
their sample of β Cephei stars, namely that the Geneva gravities generally occupy a wider
range, which they ascribed to two factors. Either the uncertainties in the Geneva gravity determination are larger, or the Str¨omgren index β may not be sensitive enough in the β Cephei
4.8 Photometric calibrations
169
Figure 4.38: Effective temperatures and gravities of the GAUDI B stars as derived from the
observed Str¨omgren colour indices c0 and β, using the Teff and log g calibrations of Balona
(1994). For early spectral types, Eq. (4.8) was used to determine β. Therefore, at the lowest
β values a curved trend is visible.
170
GAUDI B star sample
instability strip, which would lead to a rather constant log g in this range. However, considering a wide variety of B-type stars of all spectral classes, we find exactly the opposite, namely
that the uncertainty in the Str¨omgren gravity is larger. Gravities higher than 4.5 for normal
B-type main sequence stars, seem rather implausible.
The main conclusions we can draw from this comparative study is that different photometric
systems do not agree, as far as the prediction of the surface gravities concerns. The effective
temperatures do compare well.
4.8.4 Calibration based on standard stars
Both the Str¨omgren and Geneva photometric calibrations are based on the fundamental parameter determination of the stars studied in Code et al. (1976). Therefore, it is of crucial
importance that these parameters are accurate. The most recent revision of the fundamental
parameters of these stars was performed by Smalley & Dworetsky (1995). They have reevaluated the values for Teff , as derived by Code et al. (1976), to see if modern observational
data and newer model atmosphere fluxes would have any significant effect.
The result was negative: only small changes in Teff , in agreement with the Code et al. (1976)
values within the derived errors, were found. When scanning through the subsample of B
stars in Code et al. (1976), it is striking how ‘non-standard’ the stars are that were used to
derive the Teff -calibration. At first sight, these targets do not seem the most suited to derive
a reliable calibration. From the 16 B stars, there are 3 Be stars (α Eri, β Ori and ǫ Ori), 3
fast rotators (α Leo at 350 km/s, ǫ Sgr at 240 km/s and α Gru at 250 km/s), 3 spectroscopic
binaries (β Cru, α Vir and α Pav) and 1 with a peculiar spectrum (γ Crv, spectral type B8
IIIp HgMn). From the remaining 6, there are 3 supergiants (ǫ CMa, κ Ori and η CMa) and 3
giants (β CMa, γ Ori and ǫ Cen). Moreover, at least half of these stars undergo pulsations,
which can lead to variations of 1,000 K in Teff during the pulsation cycle. If we exclude all
‘special cases’, only a handful of early B-type stars remains, and thus both the Str¨omgren and
Geneva calibrations are rather fragile. For three of them, we had a spectrum available and
we derived the fundamental parameters using FASTWIND (see Section 4.5.2). We compare
them to the values derived by Code et al. (1976) and Smalley & Dworetsky (1995). Our
values agree very well with both the original Code et al. (1976), and the revised Smalley &
Dworetsky (1995) determination.
Table 4.9: Comparison of Teff and/or log g between Code et al. (1976), Smalley & Dworetsky
(1995) and this study for the three B stars in common.
Teff
Teff
Teff
(Code)
(Smalley)
(Lefever)
γ Ori 21580 ± 790 20930 ± 950 22000 ± 1000
β CMa 25180 ± 1130 24020 ± 1150 24000 ± 1000
ǫ CMa 20990 ± 760 20210 ± 950 22000 ± 1000
Star
log g
(Smalley)
3.6 ± 0.1
3.6 ± 0.1
3.1 ± 0.1
log g
(Lefever)
3.6 ± 0.10
3.4 ± 0.10
3.2 ± 0.10
4.8 Photometric calibrations
171
Figure 4.39: Comparison between the temperature indicators for both the Str¨omgren and
Geneva photometric systems: the colour indices c0 and X, respectively (upper panel), and
their resulting effective temperatures (lower panel). For both comparisons we show the linear regression curves in grey. The boxes show the range studied by Heynderickx (1991).
Open symbols in the lower panel refer to objects for which the Geneva effective temperatures
become unreliable due to extrapolation outside the (X, Y ) grid. The one-to-one relation between both derived effective temperatures is indicated as a dotted line. The indicated error on
the Str¨omgren photometry is an upper limit of 5% for a typical error of 0.01 mag in c0 and β.
172
GAUDI B star sample
Figure 4.40: Comparison between the β indices of the Geneva and the Str¨omgren system,
as well as the derived gravities within each system. Line styles have the same meaning as in
Fig. 4.39.
4.8 Photometric calibrations
173
4.8.5 Spectroscopy versus photometry
Now that we realise that large discrepancies occur among parameters derived from the different photometric systems, especially as far as the surface gravity concerns, we can have a
look if spectroscopy can bring an improvement.
Spectroscopy versus Geneva photometry
In the upper panel of Fig. 4.41, we compare the effective temperatures derived from AnalyseBstar, with those derived from interpolation in the (X, Y ) grid of the Geneva colours.
Generally, the effective temperatures agree fairly very well, within an uncertainty of 1,000 K,
at least for temperatures below 20,000 K. The hotter stars (among which a giant and a supergiant), however, lie considerably below the one-to-one relation, leaving a large discrepancy.
Unfortunately, we have too few objects in the hot temperature regime to decide whether this
trend is real. In the upper panel of Fig. 4.42, we show the same comparison, but now as a
function of the Geneva colours. Obviously, the same trends can be observed. Especially for
X > 1 (i.e. the coolest stars), the agreement between the effective temperatures is perfect.
Altogether, there seems no need to improve the Teff -calibration at this point. Even though
there is a large discrepancy for the hottest stars, we need more hot star data to be able to improve the calibration in this range. The lower panels of both Figs 4.41 and 4.42 show similar
comparisons for the surface gravities. They seem to undergo a large systematic shift of, on
the average, about 0.5 dex in log g, between photometric and spectroscopic values, the latter
being systematically lower.
Spectroscopy versus Str¨omgren photometry
In the upper panel of Fig. 4.43, we compare the spectroscopic effective temperatures (for all
the analysed GAUDI B-type stars for which Str¨omgren-indices have been measured) with the
photometric values, as derived from Eq. (4.7). We have applied Eq. (4.8) in the appropriate
regions, i.e. when c0 < −0.18 or β < 2.60, or for spectral types B3 or earlier. In the lower
panel, we compare their values as a function of c0 , which is the main temperature indicator11.
We additionally show the Be stars from the GAUDI sample, analysed by Fr´emat et al. (2006),
excluding 4 objects for which they state possibly inaccurate parameter determination. For
reasons of consistency, we have also omitted those stars, for which the analysis is based
on ELODIE spectra, as we have done for our own analysis (see discussion in Section 4.7
and Fig. 4.26). Fr´emat et al. (2006) applied a careful and detailed modeling of the stellar
spectra, taking into account the veiling caused by the envelope, as well as the gravitational
darkening and stellar flattening due to rapid rotation. The plane-parallel atmosphere models
they used for effective temperatures ranging from 15,000 K to 27,000 K were computed in
two consecutive steps.
11 Note that also β will still have some impact on the effective temperature (through Eq. (4.7)), as the effective
temperature and gravity are closely coupled.
174
GAUDI B star sample
Figure 4.41: Comparison between the spectroscopically determined Teff and log g, and the
Geneva photometric values (derived from interpolation in the (X, Y ) grid) for the few B-type
stars in GAUDI for which Geneva photometry was available. Diamonds denote supergiants,
with luminosity class (l.c.) I, triangles denote ‘giants’, with l.c. II, II-III, III or III-IV, filled
circles represent the ‘dwarfs’, with l.c. IV, IV-V or V. Open symbols are used for those objects
for which the l.c. is unknown.
4.8 Photometric calibrations
175
Figure 4.42: Position of the GAUDI stars in the calibration window of the Geneva indices.
The photometric (black) and spectroscopic (grey) determinations are connected by a vertical
line segment. Symbols have the same meaning as in Fig. 4.41.
176
GAUDI B star sample
To account in the most effective way for line-blanketing, the temperature structure of the atmospheres was computed by using Kurucz (1993) and the ATLAS9 FORTRAN program.
Non-LTE level populations were then calculated for each of the atoms considered using
TLUSTY (Hubeny & Lanz 1995) and keeping the temperature and density distributions fixed.
Model atmospheres with effective temperatures lower than 15,000 K were treated assuming
full LTE, while those hotter than 27,000 K were taken from the OSTAR2002 NLTE grid
(Lanz & Hubeny 2003).
Generally, the spectroscopically derived effective temperatures agree very well with those
derived from the photometry (see Fig. 4.43), except for some of the Be stars for which the
Str¨omgren indices, and particularly β, are probably not so reliable, as narrow-band hydrogenline photometry is strongly affected by emission. Most of the stars for which the temperatures
disagree, lie in the hot temperature regime, above roughly 20,000 K, as of which the scatter
seems to grow without any clear trend. In the same temperature range, we find an O star (open
circle), a giant for which we analysed the ELODIE spectrum and which we consider less
reliable (grey triangle), and an additional star (HD 172488) for which the spectroscopically
derived temperature seems very low with respect to its spectral type (B0.5V). In view of the
very good agreement between the spectroscopic and the photometric values, we can confirm
the Teff calibration and there is no need for improvement.
We have made a similar comparison for the surface gravities. The result can be found in
Fig. 4.44. Whereas the effective temperatures agreed very well, the situation is completely
different for the surface gravity. Only four non-Be stars agree within the ‘allowed’ error
of 0.15 dex. Three objects have a higher spectroscopic gravity, whereas, for the remaining
objects, the gravities derived from the BA94 relation are systematically higher. The transition
point seems to lie around log g ≈ 3.5. For spectroscopic gravities larger than roughly 3.5,
the photometric value is too large, whereas for lower spectroscopic gravities, the photometric
value is too low. Note that this is fully consistent with what Balona (1994) also found on the
basis of a confrontation of his calibration with evolutionary models.
Interestingly, most of the stars for which the photometric value is significantly underestimated, are Be stars. The reason for this probably lies in the fact that these stars have Hβ in
emission, which affects the β colour index. Indeed, when the core of Hβ is slightly refilled
by emission, the line suddenly appears much weaker, implying a lower log g value.
Stars which are classified as giants (4 targets) do have lower gravities (lower panel Fig. 4.44),
which gives confidence in our derived values. However, one has to keep in mind that spectral
classification is very heterogeneous, and for a lot of these stars, spectral types are uncertain.
It seems that, among the sample, there are many giants, which were classified as dwarfs or
had no luminosity class assigned.
The top panel of Fig. 4.45 illustrates the position of our sample stars (black filled circles)
with the position of the eclipsing binaries of Andersen (1991) on which the currently applied
calibration is based. The gravity from these eclipsing binaries is based on their masses and
radii, which can usually be obtained in a very precise and reliable way. However, the colour
indices come from the combined light, using the assumption that most systems have (nearly)
4.8 Photometric calibrations
177
identical components. This implies that the uncertainties in the colours may be large. On
the other hand, the colours for our GAUDI sample should be well determined, if they are
single stars, whereas the uncertainties for the derived gravities may be a bit larger than for
the binaries.
We have also investigated the position of a handful of targets for which an asteroseismic
determination of log g was available (from the derived radius and mass). Such an estimate
results from the construction of a structure model for the stellar interior that fits oscillation
frequencies with an extremely high precision. At present there are only 6 such targets, namely
HD 129929 (Dupret et al. 2004), β CMa (Mazumdar et al. 2006), δ Ceti (Aerts et al. 2006),
12 Lac (Handler et al. 2006), θ Oph (Briquet et al. 2007) and β Cep (C. Aerts, in preparation). For each of these six targets, we obtained the Str¨omgren β-index from the Galactic
β Cephei catalogue of Stankov & Handler (2005). The position of these stars in the top panel
of Fig. 4.45 (dark grey circles) seem to prolonge the curved trend of the eclipsing binaries
towards lower β-values. Unfortunately, seismic estimates of log g are not yet available for
cooler stars and/or further evolved objects.
A new calibration for log g does seem appropriate, but is hampered by the large scatter in
the data. At first sight, no clear correlation between log g and β is visible in the top panel of
Fig. 4.45. Note, however, that this figure includes a dependency of the surface gravity on the
effective temperature, via the colour index c0 in the calibration Eq. (4.9), as log g is dependent
on both β and c0 : log g(β, c0 ). The 3D representation of this dependence is given in the lower
panel of Fig. 4.45. The height of the bars is representative for log g, while the width reflects
the errors on the observed colour indices. Projecting this figure onto the β-plane (fixing c0 )
brings us back to Fig. 4.44. In Fig. 4.46, we make a similar projection, but now onto the
c0 -plane (fixing β). These figures illustrate that fixing a connection between log g, β and
c0 into one formula, is far from obvious. We have made several attempts to calibrate log g
as a function of the Str¨omgren photometric indices, going from simple linear calibrations to
complicated higher order polynomials (also in the combined terms), but none of them gave a
satisfying and convincing fit. This can be for two reasons: either there simply does not exist a
single relation which is able to cover all luminosity classes at the same time (from dwarfs up
to supergiants), or we were unable to detect it, due to the fact that we have too few datapoints
for which both the gravity and the observed colour indices are very strictly defined. In the
second case, it would be interesting to see how the trend which is now initiated by the 6
asteroseismic dwarfs, would be continued when giants and cooler dwarfs would be added, as
these provide us with very accurate values of log g. On the other hand, it is worthwhile to
investigate if the colours would change significantly when decoupling the combined light of
the eclipsing binaries, as also for this method, log g is very well defined.
In any case, the message from this investigation should be that one should refrain from using
the calibration outlined by Balona (1994), as long as the discrepancies with the spectroscopic
values are not clarified, as the photometry generally leads to an overestimation of the surface
gravity. We have clearly shown that the photometric colours are too rough a measure to
provide an accurate estimate of the stellar gravity and that they cannot compete with accurate
line fitting through high-resolution spectroscopy.
178
GAUDI B star sample
Figure 4.43: Comparison between the spectroscopically determined effective temperatures
and the Str¨omgren photometric values, derived from the application of the BA94 relation.
Crosses indicate the Be stars from Fr´emat et al. (2006). The grey symbols denote two stars
for which we have used an ELODIE spectrum, and for which the parameter determination is
less reliable. The dark grey symbol represents one of the two O-type dwarfs analysed by R.
Mokiem, the only one for which the photometric indices were available. Circles, triangles
and diamonds have the same meaning as in Fig. 4.41. The filled symbols in the lower panel
refer to the spectroscopic values for Teff , the open symbols to the photometric values.
4.8 Photometric calibrations
179
Figure 4.44: Comparison between the spectroscopically determined surface gravities and the
Str¨omgren photometric values, derived from the application of the BA94 relation. Symbols
and colours have the same meaning as in Fig. 4.43. Bottom: stars with unknown luminosity
class are now indicated by squares, to distinguish them from the ‘known’ spectral types (from
SIMBAD): supergiants (diamonds), giants (triangles) and dwarfs (filled circles).
180
GAUDI B star sample
Figure 4.45: Top: Position of the 35 analysed GAUDI B stars (filled circles) as compared
with the eclipsing binaries (open circles) of Andersen (1991) on which the log g calibration
of Balona (1994) is based. We additionally show the 6 B-type dwarfs for which we have a
very accurate asteroseismic determination of log g.
Bottom: 3D representation, showing the variation of the spectroscopically derived surface
gravity as a function of the colour indices β and c0 . The height of the shaded area’s gives
the value of the surface gravity and the width is representative for the error on the observed
colour indices.
4.8 Photometric calibrations
181
Figure 4.46: Comparison between the photometric and spectroscopic determinations of the
surface gravity as a function of the colour index c0 . Symbols and colours have the same
meaning as in the bottom panel of Fig. 4.44 and the top panel of Fig. 4.45.