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/00046361/201629124  
Published online  19 December 2016 
Spectroscopy Made Easy: Evolution
^{1} Department of Physics and Astronomy, Uppsala University, 75120 Uppsala, Sweden
email: piskunov@astro.uu.se
^{2} Space Telescope Science Institute, Baltimore, MD 21218, USA
email: valenti@stsci.edu
Received: 15 June 2016
Accepted: 2 August 2016
Context. The Spectroscopy Made Easy (SME) package has become a popular tool for analyzing stellar spectra, often in connection with large surveys or exoplanet research. SME has evolved significantly since it was first described in 1996, but many of the original caveats and potholes still haunt users. The main drivers for this paper are complexity of the modeling task, the large user community, and the massive effort that has gone into SME.
Aims. We do not intend to give a comprehensive introduction to stellar atmospheres, but will describe changes to key components of SME: the equation of state, opacities, and radiative transfer. We will describe the analysis and fitting procedure and investigate various error sources that affect inferred parameters.
Methods. We review the current status of SME, emphasizing new algorithms and methods. We describe some best practices for using the package, based on lessons learned over two decades of SME usage. We present a new way to assess uncertainties in derived stellar parameters.
Results. Improvements made to SME, better line data, and new model atmospheres yield more realistic stellar spectra, but in many cases systematic errors still dominate over measurement uncertainty. Future enhancements are outlined.
Key words: stars: abundances / radiative transfer / stars: fundamental parameters / stars: atmospheres / techniques: spectroscopic
© 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 largescale 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.
SME^{1} 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 equationofstate, 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 IDL^{2}. 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 microturbulence, vsini, and instrumental profile). SME reads line data in the format provided by VALD^{3} (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.
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 microturbulence) 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 MarquardtLevenberg 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 SahaBoltzmann equations for that purpose restricting the applications to F and hotter stars. The latest version relies on a selfconsistent equilibrium solver applicable across wide range of temperature and pressure. The new equationofstate (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 collisiondominated stellar photospheres.
The EOS solver is based on work by Phil Bennet (1991). Let Z_{elem} = n_{elem}/n_{total} be the abundance of a particular chemical element, where n_{elem} is the number density of nuclei of that element in any form (atomic or molecular, neutral or ion) and n_{total} is the number density of all elements in any form. Note that the reference is all elements, not just hydrogen. Let n_{species} be the number density of a particular atomic or molecular species with charge q_{species}. Let X^{species} 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 n_{species} 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 n_{neutral} is the number density of neutral counterpart of a qtimes 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.
List of molecular species and negative ions included the EOS.
Equations (4), (5) can be combined to provide explicit expressions for n_{species}(6)After eliminating n_{species} 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 (n_{atom}) 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 D_{species} is a dissociation energy.
The ionization equilibrium constant I for two consecutive ionization stages is given by the wellknown Saha equation: (8)where E_{ion} 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. Kurucz^{4}. 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 NewtonRaphson 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 NewtonRaphson 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 timeconsuming 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 precompute 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 nonhydrogen 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 welldeveloped Lorentzian wings.
Barklem et al. (2000a) computed twoparameter van der Waals approximations for several thousand strong transitions of elements from Li to Ni, better reproducing line shapes in the lowtemperature, highpressure regime. Their approach takes advantage of the quantummechanical description of collisions with neutral hydrogen developed by Anstee, Barklem and O’Mara (ABO, see references in Barklem et al. 2000a). The collisional crosssection σ can be converted to the broadening of the Voigt profile for a given collisional velocity v as: (10)where v_{ref} is the reference velocity value and N_{H} 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 database^{5}. 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 selfbroadening 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 reevaluated every time the atmospheric structure, abundances or line parameters are changed although in the last two cases only the opacities for affected transitions are recomputed.
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 RungeKutta, 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 planeparallel 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 rewritten 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 higherorder 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/underestimating 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.
Fig. 1 Geometric representation of the sphericallysymmetric 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θ. 
The planeparallel 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 T_{eff}, 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 T_{eff} in the subgrid that bracket the requested T_{eff}.
Next SME pairwise 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 T_{eff} 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 T_{eff}, producing a single output atmosphere at the specified T_{eff}, 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} goodnessoffit 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 T_{eff}, 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 solartype atmospheres in an ATLAS9 grid, the maximum errors are 1.4% for interpolation along the T_{eff} 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 linebyline 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.
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 T_{eff}, 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. 
3. NLTE calculations with SME
SME can account for nonlocal 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 precalculated 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 precomputed 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 b_{i} 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 b_{l} and b_{u} 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 wellcharacterized 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., T_{eff}, 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 (T_{eff}), 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 signaltonoise 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.
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 T_{eff}, 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 crossdependence of parameters and deficiency of numerical derivatives. 
Stellar parameters can be grouped into a subset that mainly affects spectral line strength (T_{eff}, 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 T_{eff}, 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 mainsequence stars cooler than T_{eff} = 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 ironpeak 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 ironpeak 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 T_{eff}. For each star in the sample, they used this relationship to tie macroturbulence to T_{eff} 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 semimajor 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 bestfitted 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 curveofgrowth 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 LevenbergMarquardt 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 T_{eff}, 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 T_{eff} 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 T_{eff} 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.
Fig. 4 Examples of normalized cumulative distribution functions of parameter offsets Δp_{i} for an SME fit. The xaxis is parameter offset Δp relative to the value returned by SME. The horizontal dash and dashdotted lines show the median and the ±1σ levels. 
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 (∂F_{i}/∂p ≠ 0) and the observed minus model residual (R_{i}) is less than 5 times the measurement uncertainty for each observed spectrum pixel. For each sensitive pixel, we estimate the parameter change (Δp_{i}) needed to make R_{i} equal to zero, assuming that R_{i} is a linear function of p. Mathematically, (16)If model errors dominate, then the width of the distribution of Δp_{i} 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 C(Δp), the normalized cumulative distribution function of Δp_{i} values. Then σ_{lo} = Δp where C(Δp) = 0.16, σ_{hi} = Δp where C(Δp) = 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 crosstalk between parameters. This is well illustrated in Fig. 5 showing the slice of the χ^{2} space along the T_{eff}log g plane. Doubleline 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.
Fig. 5 Slice of the χ^{2} space along the T_{eff}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. 
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 GaussSeidel iterations to account for scattering. Proper inclusion of scattering is critically important for lowabsorption environments such as the atmospheres of metalpoor 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 precomputed 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.
Available at http://www.stsci.edu/~valenti/sme.html
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
 Barklem, P. S., & Collet, R. 2016, A&A, 588, 96 [Google Scholar]
 Barklem, P. S., Piskunov, N., & O’Mara, B. J. 2000a, A&AS, 142, 467 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Barklem, P. S., Piskunov, N., & O’Mara, B. J. 2000b, A&A, 363, 1091 [NASA ADS] [Google Scholar]
 Bennett, Ph. 1991, Ph.D. Thesis, University of British Columbia, Vancouver, Canada [Google Scholar]
 Bergemann, M., Ruchti, G. R., Serenelli, A., et al. 2014, A&A, 565, 89 [Google Scholar]
 Brewer, J. M., Fischer, D. A., Basu, S., Valenti, J. A., & Piskunov, N. 2015, ApJ, 805, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Bruntt, H., Basu, S., Smalley, B., et al. 2012, MNRAS, 423, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Carlsson, M. 1986, Uppsala Astronomical Observatory Report 33: A Computer Program for Solving Multilevel NonLTE Problems in Moving or Static Atmospheres [Google Scholar]
 Castelli, F., & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, Uppsala, Sweden, eds. N. Piskunov, W.W. Weiss, & D.F. Gray. (ASP), Proc. IAU Symp., 210, A20 [Google Scholar]
 de la Cruz Rodríguez, J., & Piskunov, N. 2013, ApJ, 764, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Gray, D. F. 2008, The Observation and Analysis of Stellar Photospheres (Cambridge: Cambridge University Press) [Google Scholar]
 Gruyters, P., Lind, K., Richard, O., et al. 2016, A&A, 589, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kupka, F., Piskunov, N., Ryabchikova, T. A., Stempels, H. C., & Weiss, W. W. 1999, A&AS, 138, 119 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Kurucz, R. L. 1993, CDROM 13 (Cambridge, MA: Smithsonian Astrophysical Observatory) [Google Scholar]
 Kurucz, R. L., Furenlid, I., Brault, J., & Testerman, L. 1984, NSO Atlas No. 1: Solar Flux Atlas from 296 to 1300 nm (Sunspot, NSO) [Google Scholar]
 Lind, K., Asplund, M., Barklem, P. S. 2009, A&A, 503, 541 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lind, K., Asplund, M., Barklem, P. S., & Belyaev, A. K. 2011, A&A, 528, A103 [Google Scholar]
 Lind, K., Bergemann, M., & Asplund, M. 2012, MNRAS, 427, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Mashonkina, L., Zhao, G., Gehren, T., et al. 2008, A&A, 478, 529 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Osorio, Y., Barklem, P. S., Lind, K., et al. 2015, A&A, 579, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piskunov, N., & Kochukhov, O. 2002, A&A, 381, 736 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piskunov, N., Kupka, F., Ryabchikova, T. A., Weiss, W. W., & Jeffery, C. S. 1995, A&AS, 112, 525 [NASA ADS] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P., 2002, Numerical recipes in C++: the art of scientific computing (CUP) [Google Scholar]
 Sauval, A. J., & Tatum, J. B. 1984, ApJS, 56, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Shulyak, D., Tsymbal, V., Ryabchikova, T., Stütz, Ch., & Weiss, W. W. 2004, A&A, 428, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sitnova, T. M., Mashonkina, L. I., & Ryabchikova, T. A. 2013, Astron. Lett., 39, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Sozzetti, A., Torres, G., & Charbonneau, D. 2007, ApJ, 664, 1190 [NASA ADS] [CrossRef] [Google Scholar]
 Stehlé, C. 1994, A&AS, 104, 509 [NASA ADS] [Google Scholar]
 Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Valenti, J. A., & Piskunov, N. 1996, A&AS, 118, 595 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valenti, J. A., Fischer, D. A., Marcy, G. W., et al. 2009, ApJ, 702, 989 [NASA ADS] [CrossRef] [Google Scholar]
 Vidal, C. R., Cooper, J., & Smith, E. W. 1970, J. Quant. Spec. Radiat. Transf., 10, 1011 [Google Scholar]
All Tables
All Figures
Fig. 1 Geometric representation of the sphericallysymmetric 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θ. 

In the text 
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 T_{eff}, 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. 

In the text 
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 T_{eff}, 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 crossdependence of parameters and deficiency of numerical derivatives. 

In the text 
Fig. 4 Examples of normalized cumulative distribution functions of parameter offsets Δp_{i} for an SME fit. The xaxis is parameter offset Δp relative to the value returned by SME. The horizontal dash and dashdotted lines show the median and the ±1σ levels. 

In the text 
Fig. 5 Slice of the χ^{2} space along the T_{eff}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. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.