EDP Sciences
Free Access
Issue
A&A
Volume 597, January 2017
Article Number A16
Number of page(s) 10
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201629124
Published online 19 December 2016

© ESO, 2016

1. Introduction

Strong emphasis on imaging and cosmology over the last 30 yr has produced a whole generation of astronomers with a superficial understanding of stellar spectra. Recent new directions associated with the discovery of exoplanets, asteroseismology, and large-scale galactic surveys (e.g., Gaia) has rekindled interest in stars, with an emphasis on deriving accurate stellar parameters (temperature, mass, chemical composition, rotation, etc.) in a coherent fashion for a large number of targets. Spectroscopy Made Easy (SME) is a tool that we made for exactly this purpose. We emphasize the word “tool” to avoid any misconception: SME helps when used properly, but as with any tool it can be misused.

SME1 has evolved substantially since it was first released nearly 20 yr ago. In the next section we give a short introduction to SME, followed by a more detailed description of the main components that have been added or significantly changed since the original version. These are the the equation-of-state, the continuous and line opacity calculations, and the radiative transfer solver. In Sect. 4 we describe and illustrate with examples some good practices to use when analyzing stellar spectra with SME. Section 5 is dedicated to a methodology for estimating uncertainties of the resulting stellar parameters.

2. SME under the hood

SME consists of two components that are loosely connected through “input” and “output” data structures, that are usually stored in files. The first component is the graphical user interface (GUI) written entirely in IDL2. The GUI helps users create an input structure that specifies the goal of a calculation and examine the results of a calculation stored in an output structure. An SME calculation can be a straightforward spectral synthesis or a more elaborate fit of an observation. In all cases the user must provide some basic information: global stellar parameters (effective temperature, surface gravity, metallicity, and radial velocity), line data, spectral intervals, and line broadening parameters (macro- and micro-turbulence, vsini, and instrumental profile). SME reads line data in the format provided by VALD3 (Piskunov et al. 1995; Kupka et al. 1999). Fitting requires additional information: an observed spectrum and the set of parameters to vary. The GUI checks the input data for completeness and consistency, making it a convenient starting point for using SME.

Table 1

List of changes to the SME since the original paper Valenti & Piskunov (1996).

Calculations are performed by the second SME component, known as the “solver”. IDL code fits observed spectra, calling a dynamically linked external library to perform the spectral synthesis. The IDL part reads the input structure and passes all relevant information to the library. Library functions solve for molecular and ionization equilibrium, compute continuous and line opacities, and calculate intensity spectra for specified limb angles. The IDL code integrates intensities over the stellar disk and optionally solves for free parameters that yield the best fit. SME is capable of fitting global stellar parameters (effective temperature, surface gravity, metallicity, abundances of specified atoms, vsini, radial velocity, macro- and micro-turbulence) and some parameters of spectral lines. It is important to remember that SME minimizes a weighted χ2 statistic, where the usual contribution of each pixel is further weighted by the observed spectrum. This formulation gives more weight to points near the continuum at the expense of points in line cores, which helps the minimization procedure decouple the influence of continuum and line parameters. SME also calculates and returns standard χ2.

Table 1 lists most of the changes that occurred since the first SME paper. The following section describe major modifications and enhancements. Since this is more of a “progress report” paper we will not repeat the description of some important SME features, such as the observation mask, instrumental broadening, χ2 evaluation, and Marquardt-Levenberg optimization among others. More information can be found in the original paper and in documentation distributed with the package. We also emphasize SME’s compatibility with VALD as a source of atomic and molecular line data.

2.1. Equation of state

Accurate spectral synthesis requires knowledge of the number density of absorbers and perturbers at every level of the stellar atmosphere. Early versions of SME (before around 1996) use Saha-Boltzmann equations for that purpose restricting the applications to F and hotter stars. The latest version relies on a self-consistent equilibrium solver applicable across wide range of temperature and pressure. The new equation-of-state (EOS) library functions in SME solve for chemical and ionization equilibrium, returning number density and partition function for each species that affects thermodynamic state or contributes to opacities. The output of the EOS is also used for computing collisional broadening. The input parameters are temperature, total gas pressure, and atomic abundances. The EOS provides a consistent foundation for calculations of line and continuum opacities for cool stars with a variety of species competing for the same atoms. The equilibrium assumption is appropriate for collision-dominated stellar photospheres.

The EOS solver is based on work by Phil Bennet (1991). Let Zelem = nelem/ntotal be the abundance of a particular chemical element, where nelem is the number density of nuclei of that element in any form (atomic or molecular, neutral or ion) and ntotal is the number density of all elements in any form. Note that the reference is all elements, not just hydrogen. Let nspecies be the number density of a particular atomic or molecular species with charge qspecies. Let Xspecies be the count of nuclei in a particular species, that is, 0 for electrons, 1 for atoms or ions, 2 for diatomic molecules, etc. Let be the count of nuclei of the specified element in a particular species, e.g., , , and . With these definitions, we can write N equations for N elements, expressing our abundance definitions: (1)These abundance equations are supplemented by the particle and charge conservation equations: Eliminating nspecies from Eqs. (1)–(3) requires a relationship between the number density of a species and its constituents. In “equilibrium” the ratio of a reaction product to the reactants depends only on temperature. This ratio is denoted by K(T) for chemical equilibrium and I(T) for ionization equilibrium: where nneutral is the number density of neutral counterpart of a q-times ionized species. For example, these equations: connect the number of water molecules with its constituent atoms and the numbers of neutral and ionized titanium monoxide molecule with the concentration of free electrons.

Table 2

List of molecular species and negative ions included the EOS.

Equations (4), (5) can be combined to provide explicit expressions for nspecies(6)After eliminating nspecies using Expression (6) above that combines chemical and ionization equilibrium Eqs. (4), (5), we are left with only N + 1 independent variables: number densities of neutral atoms (natom) for N elements and the number density of free electrons (). This seems to be smaller than the number of equations (N + 2), but the sum of all N equations in Eq. (1) is unity by definition, reducing by one the number of constraints restoring the match with the number of unknowns.

The chemical equilibrium constant K can be written as combination of masses, partition functions, and temperature: (7)where k is the Boltzmann constant, m is a mass, Q is a partition function, and Dspecies is a dissociation energy.

The ionization equilibrium constant I for two consecutive ionization stages is given by the well-known Saha equation: (8)where Eion is the ionization energy from ion to ion+1. Combining Saha equations permits writing the expression of I connecting higher ionization stages with the neutral form.

Currently, the SME library includes partition function for 99 atoms from Hydrogen to Einsteinium with six or more ionization stages for all abundant elements. These partition functions were generously provided by R. L. Kurucz4. For molecules we collected partition functions and molecular equilibrium constants for 197 molecules and six negative ions, mostly from Sauval & Tatum (1984). Barklem & Collet (2016) recomputed molecular partition functions for diatomic molecules and a few negative ions, using updated energy level data from the literature and the NIST database. Their calculations cover temperatures from 10 K to 10 000 K, suitable even for interstellar medium conditions. We fit their finely tabulated partition function values with a function that minimizes the maximum error over the full temperature range. Table 2 lists the molecular species currently in SME and the source of each partition function.

SME uses a Newton-Raphson approach with Ng acceleration to solve Eqs. (1)–(8) iteratively. Convergence of this method depends critically on having a good initial guess. SME uses the following elegant and remarkably robust algorithm to obtain an initial guess. Apportion the input total pressure to atomic species according to their relative abundance. Use the equilibrium constants to calculate partial pressures (Eq. (6)) of most abundant species. Scale the atomic partial pressures by the ratio of the input pressure divided by the sum of current atomic and molecular partial pressures. Iterate the cycle above a few times until the scale factor deviates from unity by less than 1%. The result is close to the true solution for all physically possible conditions, even though the procedure ignores coupling between molecular species. With this initial guess, the Newton-Raphson solver typically achieves a precision of 10-5 in 3–5 iterations for every partial pressure.

2.2. Line opacities in SME

SME has a new and improved approximation for the atomic hydrogen line opacity (see below), while Voigt function: (9)is used for spectral lines of other species. Here a and v are the Lorentz broadening and detuning from the central wavelength expressed in units of Doppler width. Line opacity calculations are the most time-consuming part of spectral synthesis, so optimization is important. SME computes line opacity in two steps. First, we compute the central line opacity at each atmospheric layer using the number of absorbers provided by the EOS, line oscillator strength, and excitation energy of the lower level. At this step we identify the lines that will have no detectable contribution to the final spectrum and ignore them in subsequent opacity calculation. For lines other than hydrogen, we also pre-compute the Voigt function parameter a. In the second step, we multiply the central opacity by the opacity profile and determine the line contribution interval. The latter is done by comparing line opacity with continuous opacity at different wavelength offsets from line center, We assume a Voigt profile for all non-hydrogen lines. The Doppler broadening parameter of the Voigt function is a combination of thermal and microturbulent parts. The microturbulent velocity can be specified explicitly and adjusted by the SME. Three mechanisms contribute to the Lorentzian broadening of the Voigt profile: natural or radiative damping, interaction with charged perturbers (quadratic Stark effect), and with neutral perturbers (van der Waals effect). The latter effect is critical for determining the pressure and thus surface gravity and typical approximation using power low extrapolation from a broadening value at one temperature is not good enough for strong transitions with well-developed Lorentzian wings.

Barklem et al. (2000a) computed two-parameter van der Waals approximations for several thousand strong transitions of elements from Li to Ni, better reproducing line shapes in the low-temperature, high-pressure regime. Their approach takes advantage of the quantum-mechanical description of collisions with neutral hydrogen developed by Anstee, Barklem and O’Mara (ABO, see references in Barklem et al. 2000a). The collisional cross-section σ can be converted to the broadening of the Voigt profile for a given collisional velocity v as: (10)where vref is the reference velocity value and NH is number density of neutral hydrogen. A power law fit to the dependence on the mean collisional velocity adds the second parameter α to the approximation. These parameters, computed for many intermediate to strong lines in the solar spectrum, are available from the VALD database5. SME readily accepts line data in the VALD format.

Note that the Δγ6 fudge factor does not apply to these “extended” van der Waals parameters. You can still use this factor to enhance conventional γ6 values or results of the Unsöld approximation when no γ6 data is available. As theoretical values improve, enhancement becomes less necessary.

The new hydrogen line absorption code written by Barklem and Piskunov combines the self-broadening theory developed by Barklem et al. (2000b), the unified theory to model the interaction of hydrogen atoms with charged particles by Vidal et al. (1970) implemented in a code by Kurucz (1993), and the dynamic effects of ions using the model microfield method by Stehlé (1994).

Line opacities are re-evaluated every time the atmospheric structure, abundances or line parameters are changed although in the last two cases only the opacities for affected transitions are re-computed.

2.3. Radiative transfer

The radiative transfer (RT) solution is performed using an attenuation operator method with Bézier spline approximation for the source function. This algorithm was selected after many performance and precision tests for several options including the original SME Runge-Kutta, Feautrier and Hermitian solvers. A detailed description of the new algorithm can be found in de la Cruz Rodríguez & Piskunov (2013). Here we give a very short summary of the method. In a plane-parallel case the direction of light propagation is conveniently characterized by μ, which is the cosine of the angle between the local vertical and the propagation direction. The specific intensity Iλ is then given by the radiative transfer equation: (11)where τλ is the monochromatic optical depth measured from the surface into the star and Sλ is the source function. The boundary condition is set deep in the atmosphere where the intensity can be assumed equal to the local black body radiation. The formal solution of RT Eq. (11)from point τ0 to τ is given by: (12)where Δτ = ττ0. Note, that this equation is invariant with respect to the direction of τ and can be easily re-written in terms of stellar radius or column mass as independent variables. Assuming that the source function is changing smoothly along the integration path we can approximate it with an analytical function and evaluate the integral analytically. The simplest linear approximation requires the knowledge of the source function at the two ends of the interval. Numerical analysis shows that higher-order approximations significantly improve precision without refining step size. A quadratic fit to the source function makes this solver comparable with or superior to the Feautrier method (see Piskunov & Kochukhov 2002, Fig. 3). A simple parabolic fit will require the value of the source function at three grid points and it may result in “overshooting” – creating unphysical values inside the interval. To avoid this we used a second order Bézier fit to the source function. It achieves high precision on a sparse grid without over/under-estimating the source function. Algorithmic details can be found in the paper by de la Cruz Rodríguez & Piskunov (2013) that also describes the cubic version and even the application to the polarized RT.

thumbnail Fig. 1

Geometric representation of the spherically-symmetric radiative transfer problem. The full radius mesh used for representing the variables as functions of radius is indicated along the radius vector by a set of concentric circles. The distance of rays from the stellar center is given by the impact parameter p. The set of rays that do not hit the core are represented by five parallel rays. The mesh of points used for solving the equation of radiative transfer are the crossing points between the circles and the horizontal lines. Distances along the rays are measured by z = μr where μ = cosθ.

Open with DEXTER

The plane-parallel case can be generalised to the case of spherical geometry, as presented in Fig. 1. In this case, a ray crosses each spherical shell at different μ angles. Grazing rays actually cross the atmosphere from surface to surface. The change in geometry affects the calculation of the optical path from one layer to the next and the solver. To facilitate solution of RT in a spherical geometry, SME converts height in the model to impact parameter p, relative to the reference radius. SME now includes and interpolates in a grid of spherical models.

The attenuation operator algorithm is fast and accurate. Further acceleration comes from an adaptive wavelength grid described in the original paper. The same algorithm is used to compute the continuum intensities. SME evaluates continuum at the same wavelength where the line intensity is computed avoiding continuum interpolation. When constructing the wavelength grid, SME now takes advantage of the contribution intervals described in Sect. 2.2, resulting in a substantial performance boost.

2.4. Model atmosphere grid interpolation

SME can use a single atmosphere (e.g., for the Sun), but more commonly SME generates a custom atmosphere by interpolating in a grid of model atmospheres, based on specified values of Teff, log g, and [m/H]. SME uses the following logic to identify eight “corner” models that bracket the desired model. Find values of [m/H] in the grid that bracket the requested [m/H]. Then in this subset of models, find values of log g in the subgrid that bracket the requested log g. Finally in this subset of models, find values of Teff in the subgrid that bracket the requested Teff.

Next SME pair-wise interpolates (or extrapolates) the eight corner models to produce an output atmosphere. Pairs of models at each of the four combinations of log g and Teff are interpolated to the desired value of [m/H]. These four new models are then interpolated to the desired value of log g, yielding two models at the requested [m/H] and log g. Finally, this pair of models is interpolated to the desired Teff, producing a single output atmosphere at the specified Teff, log g, and [m/H]. The logic is designed for a regular grid, but handles occasional missing models or changes in grid spacing, both of which occur in practice.

SME linearly interpolates the logarithm of atmospheric quantities (temperature, electron number density, atomic number density, mass density) versus the logarithm of the depth scale. The depth scale can be mass column or continuum optical depth at a reference wavelength (e.g., 5000 Å). By default SME interpolates atmospheres versus continuum optical depth and solves radiative transfer on a mass column depth scale because our tests suggest that this often yields better precision.

For each atmospheric quantity, SME fits the first atmosphere in a pair with the second, allowing a scalar offset in the logarithm of the depth scale and in the logarithm of each atmospheric quantity. Valenti & Fischer (2005, Sect. 4.1) illustrate why the depth scale needs to be adjusted before interpolation. A penalty function is added to the χ2 goodness-of-fit metric to discourage large offsets in the depth scale. SME constructs a common output depth scale by adopting the mean of the offsets determined for each atmospheric quantity. Both atmospheres are then shifted onto the output depth scale and combined in proportion to the offset from the grid point to the specified value of Teff, log g, or [m/H].

We evaluated the precision of SME atmosphere interpolation by comparing each atmosphere in the grid with the atmosphere obtained by interpolating the two adjacent atmospheres along one parameter axis. In these tests, the grid spacing is effectively doubled along one axis, yielding larger errors than during normal operation. We calculated the maximum fractional error (“maximum error”) in temperature at any depth in the atmosphere and the mean absolute value of the fractional error (“mean error”) in temperature at all depths. For solar-type atmospheres in an ATLAS9 grid, the maximum errors are 1.4% for interpolation along the Teff axis with a grid spacing of 250 K, 0.1% along the log g axis for 0.5 dex spacing, and 0.2% along the [m/H] axis with 0.1 dex spacing.

The SME distribution includes atmosphere grids generated with MARCS (Gustafsson et al. 2008), ATLAS9 (Castelli & Kurucz 2003), and a line-by-line code derived from ATLAS9 (Shulyak et al. 2004). Some of the grids are segregated into separate files based on microturbulence, grid spacing, or atmosphere geometry (plane parallel vs. spherical). The suite of available atmosphere grids is described in the SME Manual and associated documents distributed with the software.

thumbnail Fig. 2

Example of NLTE spectrum synthesis for the Sun. The SME calculations used departure coefficients for sodium precomputed for the MARCS atmosphere grid and interpolated to the solar Teff, log g, and [m/H]. The MULTI calculation was done separately for each sodium line using a MARCS model for the Sun. For illustration we also include the observed solar flux spectrum and the SME LTE synthesis. The insert shows a magnified view of the core of the 5890 Å line.

Open with DEXTER

3. NLTE calculations with SME

SME can account for non-local thermodynamic equilibrium (NLTE) effects on stellar parameters and individual abundances by applying departure coefficients obtained by interpolating in externally generated grids of departure coefficients. Departure coefficients are pre-calculated for all relevant energy levels of the specified element, for every layer of every atmosphere in the specified atmosphere grid, and for a range of abundances for the species in question. SME needs departure coefficients for a range of abundances because elemental abundances may deviate from the solar pattern.

SME has been used with departure coefficients for Li, C, O, Na, Mg, Al, Si, Ca, Fe, Ba, and Eu (e.g., Bergemann et al. 2014; Osorio et al. 2015; Gruyters et al. 2016; Nordlander et al., in prep.). The current distribution includes grids for Li (Lind et al. 2009), O (Sitnova et al. 2013), Na (Lind et al. 2011; Mashonkina et al. 2008), Ca (Mashonkina et al. 2008), Fe (Lind et al. 2012), and Ba (Mashonkina et al. 2008). SME documentation describes the format of the grids and the distribution includes a tool for generating grids from user departure coefficient data.

The use of pre-computed departure coefficients in SME requires associations between each transition and specific energy levels. To make these associations, we use species names, total angular momentum quantum number J, and term designation. These data are all available from VALD using long format extraction. SME still supports the VALD short format for LTE calculations. SME currently uses a file naming convention to associate each departure coefficient grid with the corresponding model atmosphere grid. Once the user selects a model atmosphere grid, the GUI identifies which transitions can be treated in NLTE and allows the user to activate NLTE treatment for the corresponding species.

For completeness, we describe here how the radiative transfer solver incorporates departure coefficients into line opacities and the source function. SME computes the continuum in LTE, which seems to be a good approximation for cool stars, but may be a limiting factor for hot stars. Departure coefficients bi are defined as the ratio of NLTE level populations to LTE level populations, . The line opacity is then computed as: (13)where is the LTE line opacity and bl and bu are the departure coefficients of the lower and the upper energy levels. The corresponding change to the line source function is: (14)and the total source function (ignoring scattering) is: (15)Figure 2 shows the comparison between NLTE SME spectral synthesis based on interpolation of departure coefficients and the synthetic spectra computed with the MULTI code (Carlsson 1986).

4. Using SME to derive stellar parameters

The use of SME to derive stellar parameters involves several steps.

  • 1.

    Define spectral intervals that can adequately constrain themodels.

  • 2.

    For each spectral interval, obtain atomic and molecular line data from, e.g., VALD. If the stellar sample spans a wide range of temperatures, use VALD extract stellar mode to obtain line data at multiple characteristic temperatures.

  • 3.

    Verify atomic or molecular line parameters, using the spectrum of a well-characterized star (e.g., the Sun) as a constraint.

  • 4.

    Compute an initial synthesis, solving for radial velocity and, perhaps, continuum normalization. The current version of SME applies continuum adjustments to the observation, rather than the model.

  • 5.

    Use the GUI to adjust the mask, marking continuum regions that will be used to renormalize observations, line regions, which will be used to constrain model parameters, and bad regions, which will be ignored in the fit.

  • 6.

    Solve for global parameters (e.g., Teff, log g, [m/H]), keeping individual abundances fixed.

  • 7.

    Solve for individual abundances, keeping global parameters fixed, unless certain that the spectral intervals can adequately constrain global parameters and individual abundances simultaneously.

The details of how to perform these steps are described in the SME documentation.

SME can solve for various model parameters by fitting an observed spectrum. The stellar parameters are effective temperature (Teff), surface gravity (log g), metallicity ([m/H]), individual elemental abundances, microturbulence, macroturbulence, projected equatorial rotation velocity (vsini), and radial velocity. The atomic and molecular parameters are individual oscillator strengths (log gf), individual van der Waals damping parameters (γ6), and a global scale factor for van der Waals damping (Δγ6).

Finally, SME can fit some properties of the observations, such as continuum level and spectral resolution.

Individual model parameters can be well or poorly constrained, depending on the information content in the observed spectrum and physical properties of the star. Subsets of model parameters may be partially or completely degenerate with each other. Spectral line diversity is a key factor, but spectral resolution and signal-to-noise ratio are also relevant. Effective use of SME requires careful assessment of constraints in the observed spectrum, when deciding which SME parameters to determine. The sensitivity to specific free parameters (or the lack of) can be verified by trying different initial guesses, preferably bracketing the final solution.

thumbnail Fig. 3

Thousand guesses test, using Vesta observations with the Keck HIRES spectrometer. The axes are effective temperature, surface gravity and metallicity. On the left panel crosses show the initial guesses covering 1200 K in Teff, 0.5 in log g and 1 dex in metallicity. A tight clump of circles slightly above the panel center shows the results illustrating robustness of SME optimization. The right panel shows a close up of the results. The spread is very small but has a structure reflecting cross-dependence of parameters and deficiency of numerical derivatives.

Open with DEXTER

Stellar parameters can be grouped into a subset that mainly affects spectral line strength (Teff, log g, [m/H], abundances, microturbulence) and a subset that mainly affects line shape (macroturbulence, vsini). Instrumental resolution also affects line shape. Degeneracies tend to be stronger within these subsets and weaker across the two subsets.

Standard spectroscopy texts (e.g., Gray 2008) discuss how different types of spectral lines respond to changes in stellar Teff, log g, and abundance. SME can only disentangle these physical parameters if the observed spectrum includes adequate constraints. For example, Valenti & Fischer (2005) use the damping wings of the Mg 1 b triplet to constrain log g in main-sequence stars cooler than Teff = 6200 K. The damping wings disappear at higher temperature (due to Mg ionization) or lower gravity (due to lower photospheric density), limiting the effectiveness of this gravity diagnostic outside the domain where it was originally used. Ionization balance provides another constraint on photospheric density and hence gravity. Including spectral lines from neutral and ionized species helps SME constrain gravity in warmer or more evolved stars (e.g., Brewer et al. 2015).

SME uses [m/H] to interpolate in an atmosphere grid and also to scale all elemental abundances. When SME determines [m/H] by fitting a spectrum, the result is a compromise between the value that yields the best atmosphere and the value that yields the best scaled abundance pattern. In principle there should be no difference, but in practice SME and the atmosphere grids likely have slightly different assumptions and errors.

Care is required when solving for [m/H] and individual abundances of selected elements. Effectively, [m/H] becomes the abundance scale factor for any remaining elements constrained (perhaps poorly) by the observed spectrum. For example, solving for individual abundances of some iron-peak elements (V, Cr, Mn, Fe, Co, Ni) while using [m/H] to determine scaled solar abundances for some α elements (O, Ne, Mg, Si, S, Ar, Ca, Ti) would yield an [m/H] that is formally inconsistent with metallicity as defined in α-enriched atmosphere models. A more consistent approach would be to solve for individual abundances of α elements, while using [m/H] to determine scaled solar abundances for iron-peak elements.

For slowly rotating stars, distinguishing macroturbulence from vsini is difficult, especially when observed spectral lines are asymmetric due to surface granulation. In this case, it may be necessary to assume a value for one of these parameters. For example, Valenti & Fischer (2005) set vsini = 0 and determined “macroturbulence” for many stars. Most stars had some contribution from rotation, but a few had negligible rotation. This lower envelope defined a relationship between macroturbulence and Teff. For each star in the sample, they used this relationship to tie macroturbulence to Teff and then solve only for vsini.

External constraints can help breaking spectroscopic degeneracies. For example, Sozzetti et al. (2007) forced log g to be consistent with stellar structure models, constrained by the ratio of semi-major axis to stellar radius inferred from the light curve of a transiting planet. Valenti et al. (2009) forced log g to be consistent with stellar models, constrained by spectroscopic parameters from SME, apparent visual magnitude, and a precise parallax. External constraints on log g are particularly useful because gravity has a subtle effect on spectra (see, for example, the comparison of various methods by Bruntt et al. 2012).

Several applications of SME focused on strictly differential analysis of similar stars. In the case of solar twins, such a comparison can be simplified by selecting the best-fitted spectral regions and tuning the corresponding line data (oscillator strength and van der Waals broadening) using solar spectra. This hides inconsistencies in the data and the model shifting the derived parameters into the reference frame of the Sun. Of course, the downside of such approach is that accuracy of the results will progressively decrease for objects less similar to the Sun.

Since SME fits the synthesis to the observations in individual pixels there is no obvious way to derive microturbulent velocity as in the curve-of-growth analysis: by forcing the abundance to be independent of of reduced equivalent width. This condition can be checked after the SME fit, though. Using a number of lines with good parameters covering a significant range in equivalent width in SME fit will usually come close to this condition.

5. Deriving realistic uncertainties

When assessing uncertainties in derived parameters, one must distinguish between numerical errors in the solver, measurement errors in the data, and physical errors in the model. For observed spectra of good quality, physical errors in the model typically dominate other sources of error, yet formal uncertainties returned by the Levenberg-Marquardt assume measurement errors in the data dominate. In particular, our formal uncertainties are the square root of elements on the main diagonal of the covariance matrix, which is the inverse of the curvature matrix (Press et al. 2002). These formal uncertainties typically underestimate the true error (χ2 ≫ 1), so SME also returns a heuristic uncertainty estimate that considers model errors.

Before describing our algorithm for estimating model errors, we first demonstrated that numerical errors due to the SME solver are negligible. Components of the solver that can produce numerical errors include the minimization algorithm, calculation of numerical derivatives, interpolation algorithms, and ultimately numerical precision. As a straightforward test, we used SME to recover stellar parameters for 1000 simulated solar observations, all based on a single spectrum generated by SME. We set instrumental resolution to 100 000 and added 0.5% noise to each simulated observation. In this case, the SME model is perfect by construction. For each realization we adopted a random (uniformly distributed) initial guess for the free parameters. The initial guesses spanned 1200 K in Teff, 0.5 dex in log g, and 1 dex in [m/H]. The distribution of final values had a standard deviation of 1 K for Teff and 0.001 dex for the other two free parameters. This precise convergence from a wide range of initial values demonstrates that the SME solver is robust. However, the final parameter errors are much smaller than errors encountered when analyzing actual spectra.

Next we used SME to fit an actual solar spectrum obtained by observing the asteroid Vesta with the Keck HIRES spectrometer. We used an FTS solar atlas (Kurucz et al. 1984) to normalize the stellar continuum. Again we used SME to recover stellar parameters, this time starting from a uniform grid of initial guesses for the free parameters. The distribution of final values had a standard deviation of 12 K for Teff and 0.01 dex or less for log g and [m/H]. This is an order of magnitude larger than the previous test. If the spectral intervals poorly constrained the free parameters, then degeneracies and parameter errors would be even larger. Figure 3 illustrates the initial and final values of the free parameters. The final values are mostly coalesced into a primary sequence, but two alternate sequences indicate adjacent local minima. The extent of the parameter sequences indicate partial degeneracy between the free parameters, which is apparently exacerbated by model errors. These model errors include imperfect atomic line data, poor treatment of granulation, neglect of NLTE effects, the exact shape of the instrumental profile, etc.

thumbnail Fig. 4

Examples of normalized cumulative distribution functions of parameter offsets Δpi for an SME fit. The x-axis is parameter offset Δp relative to the value returned by SME. The horizontal dash and dash-dotted lines show the median and the ±1σ levels.

Open with DEXTER

From version 433 onwards, SME provides an alternative parameter uncertainty estimate that assumes model errors dominate, rather than measurement error. There are many ways to assess model error, none perfect. We begin by identifying pixels that are sensitive to small changes in a free parameter. In practice, we select all unmasked pixels (with index i) where the partial derivative of the synthetic spectrum with respect to the free parameter is nonzero (Fi/∂p ≠ 0) and the observed minus model residual (Ri) is less than 5 times the measurement uncertainty for each observed spectrum pixel. For each sensitive pixel, we estimate the parameter change (Δpi) needed to make Ri equal to zero, assuming that Ri is a linear function of p. Mathematically, (16)If model errors dominate, then the width of the distribution of Δpi values is a measure of parameter uncertainty. For normally distributed errors, the distribution is Gaussian, but model errors need not have a Gaussian distribution. Instead we construct Cp), the normalized cumulative distribution function of Δpi values. Then σlo = Δp where Cp) = 0.16, σhi = Δp where Cp) = 0.84, and the mean is σ = (σlo + σhi)/2. This cumulative distribution function approach is robust against large residuals, as long as they comprise fewer than 15% of the points at either extreme. The extent of the central part of the cumulative distribution function provides a measure of model uncertainties. In many cases the actual parameter uncertainty is approximately bounded by the formal uncertainty (from the covariance matrix) and the model uncertainty.

Figure 4 shows examples of normalized cumulative distributions. Note that the estimated median values do not exactly match the SME result. The reason is that the error estimate procedure presented here is not taking into account the cross-talk between parameters. This is well illustrated in Fig. 5 showing the slice of the χ2 space along the Teff-log g plane. Double-line contours show the uncertainties of the numerical procedure (small circle in the center) and the new uncertainty estimates for the effective temperature and surface gravity that attempt to account for the deficiencies of the model. The thick lines show the constant value contours of the χ2 surface. While the new estimate is very close (in size) to the true 1σ contour it does not account for interdependence of the two parameters, as demonstrated by the difference in orientation. When using SME to derive stellar parameters, it is imperative to understand the underlying physics and provide balanced and (as much as possible) orthogonal constraints for the free parameters.

thumbnail Fig. 5

Slice of the χ2 space along the Teff-log g plane for the example presented in Fig. 4. Contours show 1, 2 and 3σ levels above the minimum determined with procedure described in Sect. 5. The double circle in the center shows the uncertainties given by the main diagonal elements of the covariance matrix constructed by SME at the end of iterations. A larger double ellipse comes from the new uncertainty estimate.

Open with DEXTER

6. Miscellaneous and coming features

One of “undocumented” capabilities of SME is the direct calling of the library functions. Such calls could be useful for getting more information about the final model that is not saved in the output structure (e.g., partial number densities of various species, continuous and line opacities, etc.). Alternatively, direct calls can be used for other purposes, such as computing monochromatic optical depth. The SME documentation explains now how to make direct library calls and the distribution includes examples of such calls from IDL for solving the equation of state and continuous opacity at specified wavelengths.

Another useful tool included in the latest distribution simplifies porting a mask from one SME structure to another. Such need often comes up when analysing a group of similar targets using the same spectral intervals. Both structures should include observations and radial velocities. The tool called port_mask takes care of resampling and interpolating the mask. It can work on structures loaded in IDL session or stored in files.

We have successfully tested a new radiative transfer solver, that combines the current attenuation operator algorithm with Gauss-Seidel iterations to account for scattering. Proper inclusion of scattering is critically important for low-absorption environments such as the atmospheres of metal-poor giants. At moment of writing, the new solver is still in the testing phase, but we expect it to be included in SME in the near future.

7. Concluding remarks

The current version of the SME incorporates all the capabilities needed to determine stellar parameters and chemical composition for a broad range of stars, including cool dwarfs and giants. The package includes an equation of state that solves for chemical equilibrium, a faster and more accurate radiative transfer algorithm, new grids of model atmospheres, and grids of pre-computed departure coefficients for key elements. While useful, these tools do not replace the insight of a trained researcher. The SME user must still obtain and reduce the observations, select spectral regions, adjust the mask, collect atomic and molecular data, and interpret the final results.


2

Interactive Data Language, Excelis Visual Information Solutions.

Acknowledgments

The authors are very grateful to Ulrike Heiter for writing and maintaining SME documentation, to Thomas Nordlander for contributing NLTE tools and testing the model atmosphere grid interpolation, and to Paul Barklem and Remo Collet for contributing the partition functions for atoms, molecules and negative ions as part of our collaboration on the EOS package. We would also like to thank all SME users for continuous and useful feedback on the package.

References

All Tables

Table 1

List of changes to the SME since the original paper Valenti & Piskunov (1996).

Table 2

List of molecular species and negative ions included the EOS.

All Figures

thumbnail Fig. 1

Geometric representation of the spherically-symmetric radiative transfer problem. The full radius mesh used for representing the variables as functions of radius is indicated along the radius vector by a set of concentric circles. The distance of rays from the stellar center is given by the impact parameter p. The set of rays that do not hit the core are represented by five parallel rays. The mesh of points used for solving the equation of radiative transfer are the crossing points between the circles and the horizontal lines. Distances along the rays are measured by z = μr where μ = cosθ.

Open with DEXTER
In the text
thumbnail Fig. 2

Example of NLTE spectrum synthesis for the Sun. The SME calculations used departure coefficients for sodium precomputed for the MARCS atmosphere grid and interpolated to the solar Teff, log g, and [m/H]. The MULTI calculation was done separately for each sodium line using a MARCS model for the Sun. For illustration we also include the observed solar flux spectrum and the SME LTE synthesis. The insert shows a magnified view of the core of the 5890 Å line.

Open with DEXTER
In the text
thumbnail Fig. 3

Thousand guesses test, using Vesta observations with the Keck HIRES spectrometer. The axes are effective temperature, surface gravity and metallicity. On the left panel crosses show the initial guesses covering 1200 K in Teff, 0.5 in log g and 1 dex in metallicity. A tight clump of circles slightly above the panel center shows the results illustrating robustness of SME optimization. The right panel shows a close up of the results. The spread is very small but has a structure reflecting cross-dependence of parameters and deficiency of numerical derivatives.

Open with DEXTER
In the text
thumbnail Fig. 4

Examples of normalized cumulative distribution functions of parameter offsets Δpi for an SME fit. The x-axis is parameter offset Δp relative to the value returned by SME. The horizontal dash and dash-dotted lines show the median and the ±1σ levels.

Open with DEXTER
In the text
thumbnail Fig. 5

Slice of the χ2 space along the Teff-log g plane for the example presented in Fig. 4. Contours show 1, 2 and 3σ levels above the minimum determined with procedure described in Sect. 5. The double circle in the center shows the uncertainties given by the main diagonal elements of the covariance matrix constructed by SME at the end of iterations. A larger double ellipse comes from the new uncertainty estimate.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.