© ESO, 2011

1. Introduction

Comparing the solution of the fundamental equations that describe the interior of the stars to stellar surface observations has produced impressive improvements in the knowledge of the internal structure and the evolution of stars (Cassisi 2009). A considerable part of this study has been focused on the solar model (e.g., Basu & Antia 2008) but this model appears to be applicable to other stellar regimes, in particular to those of mass and chemical composition close to the solar values. Indeed, the physics of the solar model seems well suited to reproduce the one of FGK stars of Population I and capable as well of accurately recovering their global parameters (such as luminosity and effective temperature), see Lebreton et al. (1999, 2008). On the other hand, since the discovery in 1995 of the first planet around a solar-type star (Mayor & Queloz 1995), the analysis of the F, G, and K dwarfs has received considerable attention, and different research groups started to perform detailed spectroscopic studies of large samples of stars (e.g., Santos et al. 2004; Laws et al. 2003). Moreover, one can now find data for a large number of FGK stars; see The N2K Consortium (Ammons et al. 2006).

It is important to recall now that the typical output of those spectroscopic analyses are the metalicity ([Fe/H]), the effective temperature (T_eff), and the gravity (currently log    g), although the luminosity (L) can also be determined. Note that because a large number of these dwarfs are not in binary systems, their masses cannot be obtained directly by astrometry and/or spectroscopy. Their ages can be estimated with an empirical relation involving stellar activity, but this process leads to inaccuracies (e.g., Saffe et al. 2005). The determination of unknown stellar parameters such as the mass and the age can be done by fitting the output of stellar evolution models (typically direct interpolations in grids) to the known observations in the Hertzsprung-Russel Diagram (HRD) (e.g., Laws et al. 2003), or using more sophisticated methods such as the Bayesian approaches (Jørgensen & Lindegren 2005; Nordström et al. 2004) – we suggest an alternative approach here. We stress that the knowledge of the mass and the age has a considerable impact on the theories that try to explain the formation of a star/planetary-system in the context of exoplanets (e.g., Santos et al. 2004; Johnson et al. 2007).

Let us recall also that the knowledge of the stellar HRD position is a function of the stellar mass (M), of the initial individual abundance of helium (Y) and metals (Z), and of the age (t _⋆ ). However, the physical inputs chosen to describe the stellar interior also constrain the evolution of the star in the HRD. In particular, some mechanisms that are as yet insufficiently known, such as convection, rotation, and diffusion, are dependent on free parameters. For instance, in the framework of the mixing length theory (MLT) that is currently used to model the stellar convection, two more (unknown) variables have to be considered: the mixing length parameter (α) and the overshooting (ov). Thus, to appropriately model a star, we must determine these six parameters, using the fact that the solution of the equations of the internal structure provides values for T_eff, L,  [Fe/H] , and g.

For simplicity, the analysis on the HRD is frequently performed using solar-scaled stellar models, where solar values are assumed for those parameters that lack strong observational constraints. For instance, one knows the minimum possible value for the helium abundance, which is equal to the primordial value, and is about 0.23, but the individual value is not available from observation for the FGK dwarfs (and a similar argument applies to the mixing length parameter). Currently, the helium abundance value is determined as a function of the metals, assuming the same proportion Y − Z as in the Sun. This proportion is materialized by the helium-to-metal chemical enrichment parameter  (e.g., Jimenez et al. 2003; Casagrande et al. 2007). The first grid of isochrones with non-solar scaled helium values was published by Bertelli et al. (2008). For the mixing length parameter, most parts of the grids use the same value, equal to the solar one, in all stellar models, independently of the mass, the chemical composition, and the age (e.g., Pietrinferni et al. 2007). Concerning the overshooting, for stars of masses equal to or smaller than the solar mass, it is common to consider this parameter to be zero because those stars do not exhibit a permanent convective core, and for stars of higher mass a fixed value of about 0.20H_p or 0.25H_p is used.

A question that naturally arises is whether one can consider the universality and uniqueness of the mixing length parameter, the overshooting, and the helium-to-metal chemical enrichment parameter. Many authors showed that solar-scaled values for the helium and the mixing length parameter are not adapted to the nearby visual binary star α   Centauri and for the Hyades (Gennaro et al. 2010; Casagrande et al. 2007; Eggenberger et al. 2004; Miglio & Montalbán 2005; Lebreton et al. 2001). There is also no clear answer regarding the overshooting parameter (see the studies in Ribas et al. 2000a; Claret 2007, for double-line eclipsing binary systems).

Thus, the universality of the mixing length parameter, the overshooting, and the helium-to-metal chemical enrichment parameter should be investigated in great detail. The main goal of this paper is to propose the application of a global derivative-free optimization method to estimate simultaneously the stellar mass, the initial individual abundances of helium and metals, the age, and the two convection parameters (mixing length and overshooting), taking into account the observed metalicity, the effective temperature, the gravity, and the luminosity for each star. This optimization methodology will be adapted and applied to a sample of about 115 FGK stars issued from detailed spectroscopic analysis.

We also provide an online computational tool for other researchers to use and test our methodology: www.norg.uminho.pt/aivaz/astro.

This paper is organized as follows. In Sect. 2 we describe the input physics chosen for our work and introduce our optimization methodology. We then report in Sect. 3 results for the Sun and for a sample of fictitious stars as well as for true FGK stars selected among those observed by Santos and co-workers. Finally, in Sect. 4 we discuss the results, draw conclusions, and suggest future work.

2. The estimation methodology

2.1. Input physics in stellar evolutionary models

We adjusted known estimation models for the stars studied in this work. The stellar evolution calculations were made with the CESAM code, version 3 (Morel 1997), selecting the following physical features. The nuclear reactions rates were taken from Caughlan & Fowler (1988). The OPAL opacities were the ones from Iglesias & Rogers (1996) combined at low temperatures with data from Alexander & Ferguson (1994), following the approach in Houdek & Rogl (1996). The atmosphere was described by an Eddington T(τ)-law. The convection was treated according to the MLT from Böhm-Vitense (1958), leaving the mixing length parameter (proportional to the pressure scale height H_p) as an unknown, and thus letting α and ov be free parameters. The equation of state was chosen as the so-called EFF (Eggleton et al. 1973), which is valid for solar-type stars where the departure from the ideal gas is not a relevant one (note that owing to the large amount of computation required in our work, an analytical equation of state is clearly more suitable than a tabulated one). Finally, the solar mixture was taken from Grevesse & Noels (1993). Using this input physics, we were able to reproduce to within five digits of accuracy the observed solar luminosity, radius, and metalicity (see Sect. 3.1).

We point out that there are other mechanisms that we did not consider here that could affect the HRD position of a model, such as the rotation (Maeder & Meynet 2000) and the helium and metal gravitational settling (microscopic diffusion), but one should keep in mind that the stars addressed in our work are slow rotators. On the other hand, diffusion has only a marginal effect in the HRD position for Population I solar-mass stars (e.g., Lastennet et al. 2003). It is therefore reasonable to expect that the exclusion of these mechanisms does not considerably change our results.

2.2. The estimation model

We considered the following parameters for estimation: the stellar mass M (M_⊙), the abundances of hydrogen X, helium Y, and metals Z, the stellar age t_★ in Myr, the stellar surface convection α, and the stellar nucleus overshooting ov. Because X, Y, and Z are related by Z = 1 − X − Y, only two of these variables need to be considered, and we chose to work with X and Y.

For a given star, the evolution simulation is performed for specific values of these six parameters (M, t_★, X, Y, α, and ov). Among other features obtained from the simulation are the effective surface temperature T_eff, the luminosity L, and the radius R (from which we can then obtain the stellar gravity by g = 27397M/R²).

From the observational point of view, we have an estimate for the metalicity Z/X, the luminosity, the effective temperature, and the gravity. Absolute errors are available for all these observed quantities.

We are therefore left with four observed quantities and six unknown stellar parameters. The proposed optimization methodology consists in the computation of the six unknown stellar parameters by solving a minimization problem whose objective function reflects the fit between the simulated and observed stellar features.

2.3. The optimization formulation

The minimization problem is thus defined by subject to (3)where Z = 1 − X − Y and the subscript obs and the prefix δ denote the observed data and the corresponding absolute errors, respectively. The vectors and represent lower and upper bounds on the variables, respectively.

Note that the objective function (2) is nonlinear in the optimization variables (M, t_★, X, Y, α, and ov), because its evaluation depends on the solution of the underlying stellar evolution equations. We face an inverse or parameter estimation problem of the simulation-based type, for which the objective function is expensive to evaluate and its derivatives are unavailable. Our numerical experience has also shown that this problem has non-unique global minimizers. Thus, to properly solve it, one must select a solver capable of global optimization without the use of derivatives (and, preferably, in a parallel environment, given the cost of the numerical simulations).

2.4. The optimization solver

The algorithm selected (PSwarm) to solve (2), (3) is a direct search method enhanced by particle swarm, which has been recently shown to perform well in a large variety of bound-constrained optimization problems, compared to other global derivative-free solvers (Vaz & Vicente 2007).

Direct search methods are a popular class of iterative optimization methods that try to minimize a function by comparing objective function values at a finite number of points (a point in this context represents a star configuration). Direct search methods do not use derivative information of the objective function nor do they try to approximate it (see Conn et al. 2008; Kolda et al. 2003) and typically organize their calculations at each iteration around two main steps, a poll step and a search step.

The poll step evaluates the objective function at a carefully chosen sample set (using positive spanning sets), looking for points relatively close to the current iterate that have a lower objective function value. The poll step obeys specific rules that ensure the overall algorithmic convergence to local stationary points.

The search step is essentially free from any rules, except that it must evaluate the objective function a finite number of times. The search step is optional but, when tried and incapable of finding a point with lower objective function value, must necessarily be followed by a poll step. The search step may accommodate heuristic methods for a global optimization search. The scheme used in the search step by PSwarm is particle swarm, a well-known population-based heuristic (where, again, in this context, particle or points represent star configurations).

To handle the simple bounds on the variables, one initializes the algorithm with a feasible initial guess (a star configuration that satisfies the bound constraints (3)) and replaces the objective function f(z) in the comparisons of the search and the poll steps by the extreme barrier or indicator function  (given by if z is feasible and otherwise).

PSwarm takes advantage of using a population in the search step to then poll at the best particle (the point in the population with best objective function so far), which improves the overall robustness of the population scheme. One is able to prove (see Vaz & Vicente 2007) that the algorithm is convergent to first-order stationary points from arbitrary points and, under some additional conditions, that it can eventually meet the stopping criteria used in the search and poll steps.

2.5. Implementation details

We point out that PSwarm explores the variables domain regardless of their physical meaning. Thus, it is not uncommon for the simulation code (CESAM) to fail or return an unexpected result. In these cases, one assigns  + ∞ to the objective function value (exactly as we do for infeasible points).

PSwarm was initialized with a randomly generated population (a set of star configurations satisfying (3)) and ran with the following parameters: 42 particles in the population, a precision of ϵ = 10^-5 in the stopping criteria, and a maximum of 2000 function evaluations. The remaining parameters were left as default, and we refer the reader to Vaz & Vicente (2007) for the list of default options of PSwarm and related information.

3. Tests and applications

The lower and upper bounds considered in (3) are reported in Table 1. These bounds were chosen to be representative of nearby FGK Population I stars.

3.1. The Sun

The Sun is the natural reference for modeling FGK stars. Every new stellar evolutionary code or modeling technique must be tested first using accurate solar observation data. We consequently applied the previously described approach to the Sun, using the observed data obtained from the solar calibration method described by Christensen-Dalsgaard (1991), the input physics provided in Sect. 2.1, and the commonly accepted solar age (Dziembowski et al. 1999). The results are reported in Table 2, where we display the final values estimated for the parameters M, t_★, Y, Z, α, and ov, as well as the corresponding values for T_eff, L, Z/X, and log    g compared to the observed ones. We also present the final value of the objective function f given by (2). Note that the estimate for the solar model fitted the observed luminosity, radius, and metalicity with five digits of accuracy. In this computation, 2016 function evaluations were made involving 56 poll steps, out of which 31 were successful (note that because the algorithm is stochastic, different runs can lead to different results). See also Fig. 1 for a plot representing the objective function value as a function of the number of objective function evaluations in a logarithmic scale.

Fig. 1 Optimization history of one run for the Sun.

In the next subsection we will report results of our optimization methodology on a sample of five fictitious stars. We will then incorporate the solar results in the discussion about the astrophysical quality of the solutions obtained for these stars. For the moment, we just point out that the solar mass was basically well reproduced. However, the estimated age was about 400 Myr lower, the helium was 0.013 higher, and the mixing length parameter was 0.3 lower. We also remark that we found a non-zero overshooting value for the Sun without astrophysical meaning.

3.2. Fictitious FGK stars

To further test our methodology, we generated “observed” data for five fictitious stars with the evolutionary code CESAM for predefined values of M, t _⋆ , Y, Z, α, and ov. The values of T_eff, L, Z/X, and g generated by this process are those presented in Table 3. The corresponding absolute errors (denoted by the lower scripts δ) were chosen to be representative of nearby FGK Population I stars issued from a detailed spectroscopic analysis.

We then applied our simulation-based optimization procedure (CESAM-PSwarm) to each fictitious star to estimate M, t _⋆ , Y, Z, α, and ov. As already mentioned, there is no unique solution to the corresponding optimization problems. The number of parameters to be estimated is larger than the number of available observations, and additional difficulties arise from the high non-linearity of the stellar evolution equations and the use of relatively small observational error box compared to the feasible region. To overcome these difficulties, we ran the simulation-based optimization 25 times (see below).

For the sake of illustration, we show in Table 4 the final estimated parameters M, t _⋆ , Y, Z, α, and ov for the 25 runs of the fictitious star fake1. With the exception of three cases (runs 21, 23, and 25) all final values obtained for the objective function are lower than 1, showing a good match to the central point of the observational error box. On the other hand, the solutions obtained for each stellar parameter exhibit a considerably large range of values: for instance, the mass runs between 0.72 M_⊙ and 0.94 M_⊙ and the age between 100 Myr (the lower bound) and 9806 Myr. Similar comments can be made about the other parameters.

We then computed the average values of all estimated values for each parameter (excluding the three cases mentioned above where f > 1; see next subsection). We observed that after 20 runs the average did not change much and this is why we chose to stop after 25 runs. In Table 5 we present the average of the final estimated values of M, t _⋆ , Y, Z, α, and ov for all fictitious stars, comparing them to their original predefined values (the ones used to generate the data). We also report the average of the objective function values, which stands clearly below 1.

We analyze below the differences between the estimated values and the original ones for the fake stars as well as for the Sun (see previous subsection), taking into account the estimates made by other authors.

• Mass (M). The differences obtained for the mass are lower than 0.1 M_⊙. Both Santos et al. (2004) and Laws et al. (2003) estimated the mass of FGK stars based on the interpolation of stellar evolutionary grids, with errors around 0.05 M_⊙. In this mass regime, the mass-luminosity relation has an rms of about 0.07 M_⊙. Finally, Casagrande et al. (2007), using isochrones with an helium non-solar scaled, tried to recover the mass of 12 members of spectroscopic and visual binary stars (for which the mass is well known), achieving an average accuracy of about 0.05 M_⊙.

• Age (t _⋆ ). Except for fake2, the difference between the estimated values of the age for the fictitious stars and the Sun was around or below 1 Gyr. In order to discuss this result, we point to the work of Saffe et al. (2005), who compared the age of 49 planet-hosting stars obtained by different methods: chromospheric activity, isochrones grids, lithium abundance, metalicity, and transverse velocity dispersions. The main conclusion was that the dispersion of the derived values for the age among the different methods is never lower than 2 to 3 Gyr. On the other hand, the authors also found that the internal dispersion on the isochrones method reaches from 1 to 3 Gyr (Saffe et al. 2005; Soderblom 2010). This dispersion is still marginally compatible with the greatest difference in age obtained by our approach, achieved for fake2.

• Helium abundance (Y). Our results show a difference in helium lower than 0.025. In contrast to the previous estimates, there are not too many independent helium calculations for FGK stars because of the traditional lack of the corresponding grids of values. The most recent one is given in the Casagrande et al. (2007) paper, where individual values for the helium abundance are given for the above mentioned binary members with an average accuracy of about 0.05. A few years ago, Ribas et al. (2000b) analyzed a sample of 50 detached double-lined eclipsing binaries, some of them with FGK components, by means of isochrones with non-solar scaled helium. They were able to compute individual values of helium with an accuracy of 0.04.

• Mixing length parameter (α). Our results show a difference on the mixing length parameter of about 0.3H_p. To the best of our knowledge there has never been a paper published with computed values of the mixing length parameter for more than 10 stars. Most of the calculations have been made in the context of particular objects like the Sun, binary stars (e.g., Fernandes et al. 1998; Lastennet et al. 2003; Miglio & Montalbán 2005; Torres et al. 2006), and FGK stars in Hyades (e.g., Lebreton et al. 2001; Yildiz et al. 2006). The reported accuracy in these cases was not higher than 0.2 − 0.3H_p.

• Overshooting (ov). Except for sub-solar mass models, our results show differences of about 0.1 − 0.2H_p. For sub-solar models (where no convective core is expected), our simulations exhibit a non-zero value. These non-zero values must be considered as numerical results without astrophysical relevance. Double-line eclipsing binary stars are currently used to constrain the amount of overshooting. One of the most recent works is the one by Claret (2007), where the overshooting for individual stars was estimated with an accuracy of 0.2H_p, causing a relative error of at least 50%. The overshooting seems, in deed, to be the most poorly known parameter among the six stellar ones.

The above discussion shows that our results are aligned to what has been recently achieved in stellar modeling by other techniques. Our approach has the advantage of considering as free all stellar parameters of relevance, without imposing ad-hoc assumptions on their values.

3.3. Real FGK stars

In this subsection we describe the application of our technique to a large sample of stars. Regarding the choice of the observed sample, our main concern was to select a group of stars with spectroscopic measurements for the effective temperature, metalicity, and gravity. Ideally, the observations for each star should be made by the same technique to keep the internal consistency. Proximity of the star was also a requirement to reduce luminosity errors. With these considerations in mind, we chose a sample of the FGK stars considered in Santos et al. (2004, 2005) and Sousa et al. (2006), eliminating the sub-giants and stars with errors in parallaxes higher than 10% or without HIPPARCOS data.

In order to apply our method, we needed to compute the stellar luminosity. For this purpose, we took the HIPPARCOS parallaxes and assumed that the solar bolometric magnitude was 4.75, adopting the bolometric correction in Flower (1996).

Our final sample is composed of 196 stars. As previously said, we ran our method 25 times for each star. To draw conclusions, we removed runs for which the final objective function value was greater than 1 (f > 1) or runs that lead to any variable being equal to an upper or lower bound (because the algorithm would likely converge to another solution if the bounds were removed). After this first filtering process, stars with less than 5 runs were also removed from the test set.

Because the overshooting estimation is affected by a large error (see Sect. 3.2), we decided to keep it constant in our estimations. On the other hand, for Population I stars, Z and  [Fe/H]  can be easily related to each other. Accordingly, we decided to constrain both Z and the overshooting as follows:

• Z = 0.0173  ×  10 ^[Fe/H] ;

• ov ≈ 0.0 if M < 1.1   M_⊙, and ov ≈ 0.2H_p if M ≥ 1.1   M_⊙.

In Table 6 we show the solutions (average for the remaining runs) obtained for the four stellar parameters, for each of the 115 stars identified by this selection procedure. Although there are no known results for a simultaneous estimation of these four parameters, we will attempt to compare our obtained results to previous results in the literature.

In Fig. 2 we depict a comparison between the results obtained by the mass expressed as polynomials in T_eff, log    g, and  [Fe/H] , derived by Torres et al. (2010) using binary stars. A filled line of slope one passing by the origin is plotted (and an exact match between the Torres et al. 2010, results and the obtained ones would mean that all stars would be on the filled line). Evidently, both estimates agree well.

Fig. 2 Torres et al. (2010) results vs. obtained mass results.

In Table 7 we compare the ages reported by Saffe et al. (2005) and Holmberg et al. (2009), including the errors and the percentage differences to our obtained ages for about 16 stars (consisting of those in the cited papers that also belong to our test set of 115 stars).

We can observe that the obtained results are consistent with those works. We point out that our age estimates disagree (within 1.5 Gyr accuracy) with the Saffe et al. (2005) and Holmberg et al. (2009) results for only three stars (HD 106252, HD 117618, and HD 179949). On the other hand, the mean relative differences among these three estimates are consistent pairwise: 39% (ours-Saffe), 31% (ours-Holmberg), and 36% (Saffe-Holmberg).

Casagrande et al. (2007) reports numerical results on five stars belonging to our test set. The results published there are for the estimate of the age, the helium, and the mass. We report in Table 8 our obtained results versus theirs. Casagrande et al. (2007) also report the error in the helium and the mass associated with each star. A brief analysis of Table 8 shows that our results are compatible with those published byCasagrande et al. (2007). Moreover, the star HD 3651 should not be included in this comparison, because we imposed a lower bound of 0.23 for helium, excluding from the feasible set the value reported by Casagrande et al. (2007). On the other hand, our set includes the nearest visual binary α Centauri (HD 128620/1). Our results predict 1.2   M_⊙ and 0.96   M_⊙ respectively for the mass of α Centauri A and B. These values agree with the observed values 1.105   M_⊙ and 0.934   M_⊙ (Torres et al. 2010).

To finish this section, we report the average of the errors obtained in the numerical results. Given a star j and a parameter P, the average for each parameter is computed by where is the average value reported in Table 6,  is the parameter obtained at run number i, and N_j is the number of runs (25 when all runs were considered valid). Looking at these numbers over all the stars tested, we found the following typical values for Δ_P: The stellar estimation for the FGK stars was performed in the Milipeia cluster, which is available at the University of Coimbra. The cluster is formed by two management nodes (Sun Fire X4100) and 130 processor nodes (Sun Fire X4100), running the CentOS 4.4 operating system (each node offering two double core processors). Because the objective function evaluation is obtained by simulating the evolution of a star, the corresponding CPU times can be significant for stars of different ages. While the PSwarm population was selected to be a multiple of the number of available processors, some synchronization problems can still occur (for instance, an invalid star configuration immediately detected by CESAM takes a negligible amount of time to evaluate).

4. Discussion and conclusions

We have presented a mathematical optimization methodology to estimate the stellar parameters (mass, age, helium and metal abundances, mixing length parameter, and overshooting) from the correspondent photometric, astrometric, and spectroscopic observations. The estimation was carried out by solving a simulation-based optimization problem using a global derivative-free algorithm. We made available a computational tool to the community to interface the global optimization solver (PSwarm) and the stellar evolutionary code (CESAM)¹. We tested our method on both fictitious and real FGK stars of Population I. For each star, we ran the optimizer 25 times, obtaining a sample of 25 estimated values for M, Y, Z, t _⋆ , α, and ov, from which we then took the corresponding average values. In particular, we derived the mass, age, helium, and mixing length parameters for a sample of 115 FGK stars. The comparison between our quoted mass, helium and age values and those obtained previously by different authors revealed encouragingly coherent results, mainly on the mass estimates. In a forthcoming paper, we plan to study in detail the impact of the estimated stellar parameters on the knowledge of the chemical properties of the solar neighborhood and to compare the planet-hosting stellar population and that of stars with no planet detection. We conclude from our results (in Table 7) that the average helium abundance over the sample of planet-hosting stars (around the solar value 0.29) is similar to the one over the sample of stars where no planet has yet been detected. This is a particularly interesting conclusion because it is known that the planet-hosting stars are on average richer in metal than the others. Thus, the observed metal enrichment in the solar neighborhood does not seem to be followed by one in helium. On the other hand, the averaged age is similar in both stellar populations, around 5 Gyr. Assessing the whole sample of 115 FGK stars, we did not find any clear correlation between the mass and the mixing length parameter.

Acknowledgments

Support was provided by FCT under grants POCTI-SFA-2-675, POCI/MAT/58957/2004, POCI/MAT/59442/2004, PTDC/CTE-AST/66181/2006, and PTDC/MAT/64838/2006, and by the CDS-Strasbourg Database. J.M.F. also acknowledges the support from Fundação Calouste Gulbenkian and the European Astronomical Society. We are grateful to Pierre Morel and the CESAM team for providing us the stellar evolutionary code. Finally, we would like to thank Ana Luísa Custódio (New University of Lisbon) for several interesting discussions and two anonymous referees whose comments have much improved the presentation of the paper.

References

All Tables

All Figures

	Fig. 1 Optimization history of one run for the Sun.
In the text

	Fig. 2 Torres et al. (2010) results vs. obtained mass results.
In the text