© ESO, 2012

1. Introduction

Determining accurate stellar mass is crucial to several domains in astrophysics. Indeed, mass plays a dominant role in the evolution and final fate of stars. It is also a key parameter when characterising exoplanetary systems since the masses of exoplanets, essential for determining whether they are rocky or gaseous in nature, are in general determined with respect to that of the star. However, in spite of its importance, stellar mass is usually difficult to obtain accurately for single stars and is frequently model-dependent.

A method which has often been used for estimating stellar mass is to compare the position of a star in an Hertzsprung-Russell (HR) diagram with evolutionary tracks for stellar models of different masses. The main difficulties with this method are its model dependence, the large error bars, and the fact that there are regions in the HR diagram in which more than one stellar evolution track goes through each and every point. A promising alternative is to use asteroseismology. Indeed, acoustic modes are a sensitive indicator of the mean density of a star. Current space missions CoRoT (e.g. Michel et al. 2008) and Kepler (e.g. Kjeldsen et al. 2010) are measuring stellar pulsation frequencies to high accuracy. When combined with an independent determination of the radius, this yields the stellar mass. Thanks to very precise parallax measurements, radii accurate to 2% are expected from the forthcoming astrometric mission Gaia (Perryman et al. 2001). The main challenge is then finding the best way of extracting the mean density of a star from a given set of frequencies.

Various methods have been devised or proposed for estimating stellar mean density from pulsation frequencies. These include using simple scaling laws based on the large frequency separation (Ulrich 1986) or other frequency combinations (Kjeldsen et al. 2008), searching for the best-fitting model within a grid of models or a parameter space and retaining its mean density (e.g. Bazot et al. 2008; Metcalfe et al. 2010; Kallinger et al. 2010), and performing a full structural inversion to determine the density profile, then integrating it to obtain the mean density. The structural inversion may be carried out using linear or non-linear techniques (e.g. Gough & Thompson 1991; Antia 1996; Roxburgh & Vorontsov 2002). The scaling relations have the obvious advantage of simplicity over the other methods and require fewer frequencies than full structural inversions, but are less accurate. Although more expensive computationally, searching for best-fitting models in a grid is also straightforward and usually more precise, but is model-dependent, so that the conclusions are only as firm as the input physics. Full structural inversions have the advantage of being less model-dependent, but usually involve a number of free parameters which need careful adjustment and can fail to converge when non-linear. However, with the current missions CoRoT and Kepler and future asteroseismic projects, there is a growing need to provide accurate seismic stellar masses for a large number of stars in an automatic way. An ideal method would combine the advantages of the above techniques while avoiding their weaknesses.

In this paper, we focus on kernel-based linear inversion methods that directly yield the mean density without inverting for the full density profile. As such, this represents a first attempt at constructing methods which combine the simplicity of scaling laws with the accuracy of full structural inversions. Although the results are not as accurate as initially hoped for, this approach gives insight into the precise relation between the mean density of a star and pulsation modes, and provides a better way of evaluating scaling laws. Section 2 establishes a general framework in which to construct and evaluate linear inversion techniques which yield the mean density directly. The following section then describes various linear inversion techniques, including linearised scaling laws. This is followed by a section which shows how to extend these techniques beyond the linear regime. Sections 5–7 apply these different methods to the Sun, to test cases using a grid of models, and to three stars, two of which were observed by CoRoT (HD 49933 and HD 49385). Finally a discussion concludes the paper.

2. General aspects

2.1. Starting point

As for most inversion methods, the methods described below start with an observed star and a reference stellar model. At this point we will assume that the star and the model have the same radius, but will relax this constraint later on. The structure of the reference model needs to be sufficiently close to the structure of the star so that the link between the relative frequency differences, δν_nℓ/ν_nℓ, and the differences in the stellar structure can be described by a linear relation: $\begin{array}{ccc}\frac{\mathit{\delta }{\mathit{\nu }}_{\mathit{n\ell }}}{{\mathit{\nu }}_{\mathit{n\ell }}}& \mathrm{=}& \int \begin{array}{}\mathrm{1}\\ \mathrm{0}\end{array}{\mathit{K}}_{\mathit{\rho ,}{\mathrm{\Gamma }}_{\mathrm{1}}}^{\mathit{n\ell }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\frac{\mathit{\delta \rho }}{\mathit{\rho }}\mathrm{d}\mathit{x}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}{\mathit{K}}_{{\mathrm{\Gamma }}_{\mathrm{1}}\mathit{,\rho }}^{\mathit{n\ell }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\frac{\mathit{\delta }{\mathrm{\Gamma }}_{\mathrm{1}}}{{\mathrm{\Gamma }}_{\mathrm{1}}}\mathrm{d}\mathit{x}\\ & \mathrm{+}& \end{array}$(1)where ρ is the density, Γ₁ the first adiabatic exponent, x = r/R the fractional radius, n the radial order, ℓ the harmonic degree, R the radius and $\frac{{\mathit{F}}_{\mathrm{surf}}\mathrm{\left(}{\mathit{\nu }}_{\mathit{n\ell }}\mathrm{\right)}}{{\mathit{Q}}_{\mathit{n\ell }}}$ an ad hoc surface term (e.g. Christensen-Dalsgaard & Berthomieu 1991). The function F_surf is assumed to be a slowly varying function of ν only, and Q_nℓ represents the mode’s inertia (e.g. Aerts et al. 2010), normalised by the inertia of a radial mode, interpolated to the same frequency. The kernels ${\mathit{K}}_{\mathit{\rho ,}{\mathrm{\Gamma }}_{\mathrm{1}}}^{\mathit{n\ell }}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ and ${\mathit{K}}_{{\mathrm{\Gamma }}_{\mathrm{1}}\mathit{,\rho }}^{\mathit{n\ell }}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ are known functions which can be deduced from the (n,ℓ) eigenmode of the reference model, via the variational principle (Gough & Thompson 1991). Although the above relation can also be applied to other structural pairs, we will work with the pair (ρ,Γ₁) since ρ is an obvious choice for finding the mean density, and Γ₁ is expected to vary little between different models. Throughout this paper, we will use the following definition for the relative difference of a quantity f, although we note that other definitions have been used (e.g. Antia & Basu 1994): $\frac{\mathit{\delta f}}{\mathit{f}}\mathrm{=}\frac{{\mathit{f}}_{\mathrm{obs}}\mathrm{-}{\mathit{f}}_{\mathrm{ref}}}{{\mathit{f}}_{\mathrm{ref}}}\mathrm{·}$(2)With this definition, we have the relation f_obs = f_ref(1 + δf/f).

The mass difference, δM, between the reference model and the star is given by $\mathit{\delta M}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{R}}\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathit{\delta \rho }\mathrm{d}\mathit{r}\mathit{.}$(3)From this equation it is possible to deduce the relative difference in mean density, , $\frac{\mathit{\delta }\mathit{\rho ̅}}{\mathit{\rho ̅}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{4}\mathit{\pi }{\mathit{x}}^{\mathrm{2}}\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{R}}}\frac{\mathit{\delta \rho }}{\mathit{\rho }}\mathrm{d}\mathit{x,}$(4)where the quantity ρ_R = M/R³ is used for non-dimensionalisation. It is important to note that if δρ/ρ ≡ 1 then , in other words $\mathrm{1}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{4}\mathit{\pi }{\mathit{x}}^{\mathrm{2}}\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{R}}}\mathrm{d}\mathit{x}\mathit{.}$(5)If we now assume that the radius of the reference model is different from that of the star, then the model can be rescaled through homology, i.e. a transformation which preserves all dimensionless quantities, so as to have the same radius as the star. Furthermore, the mean density needs to be kept the same to avoid modifying the pulsation frequencies when rescaling the model. Equation (3) will then give the relative mass difference between the star and the rescaled reference model. However, given their non-dimensional form and because the density is preserved between the original and scaled models, Eqs. (1) and (4) remain unchanged and can be applied directly to the original (unscaled) reference model, as long as x = 1 corresponds to both the radius of the star, and the radius of the reference model. The procedures which will be described below are based on these two equations, and therefore do not require explicitly scaling the reference model to the same radius as the star. They will consequently yield the mean density rather than the mass, although the latter can be deduced if the radius is known.

2.2. Linear inversions

The procedures described in this paper then consist in searching for a linear combination of the relative frequency differences that best reproduces the relative mean density variation, : $\frac{\mathit{\delta }\mathit{\rho ̅}\mathrm{inv}}{\mathit{\rho ̅}}\mathrm{=}\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}\frac{\mathit{\delta }{\mathit{\nu }}_{\mathit{i}}}{{\mathit{\nu }}_{\mathit{i}}}\mathit{,}$(6)where the index i represents the pair (n,ℓ). Using Eq. (1) to replace δν_i/ν_i then yields the following expression: $\begin{array}{ccc}\frac{\mathit{\delta }\mathit{\rho ̅}\mathrm{inv}}{\mathit{\rho ̅}}& \mathrm{=}& \int \begin{array}{}\mathrm{1}\\ \mathrm{0}\end{array}{\mathit{K}}_{\mathrm{avg}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\frac{\mathit{\delta \rho }}{\mathit{\rho }}\mathrm{d}\mathit{x}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}{\mathit{K}}_{\mathrm{cross}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\frac{\mathit{\delta }{\mathrm{\Gamma }}_{\mathrm{1}}}{{\mathrm{\Gamma }}_{\mathrm{1}}}\mathrm{d}\mathit{x}\\ & \mathrm{+}& \end{array}$(7)where $\begin{array}{ccc}{\mathit{K}}_{\mathrm{avg}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \\ {\mathit{K}}_{\mathrm{cross}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \end{array}$In order for to be a good estimate of the relative mean density variation, the averaging kernel, K_avg, needs to be as close as possible to 4πx²ρ/ρ_R, as can be seen by comparing Eqs. (7) and (4). Furthermore, K_cross, which represents the cross-talk coming from relative variations on the Γ₁ profile, needs to be reduced as does the last term, which represents the surface effects. Hence, these kernels provide a natural way of evaluating an inversion procedure (e.g. Christensen-Dalsgaard et al. 1990).

2.3. Various error bars

The starting point for obtaining the error on the mean density estimate is to calculate the difference : $\begin{array}{ccc}\frac{\mathit{\delta }\mathit{\rho ̅}\mathrm{inv}\mathrm{-}\mathit{\delta }\mathit{\rho ̅}\mathrm{exact}}{\mathit{\rho ̅}}& \mathrm{=}& \int \begin{array}{}\mathrm{1}\\ \mathrm{0}\end{array}\left({\mathit{K}}_{\mathrm{avg}}\mathrm{-}\mathrm{4}\mathit{\pi }\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{R}}}{\mathit{x}}^{\mathrm{2}}\right)\frac{\mathit{\delta \rho }}{\mathit{\rho }}\mathrm{d}\mathit{x}\\ & \mathrm{+}& \int \begin{array}{}\mathrm{1}\\ \mathrm{0}\end{array}{\mathit{K}}_{\mathrm{cross}}\frac{\mathit{\delta }{\mathrm{\Gamma }}_{\mathrm{1}}}{{\mathrm{\Gamma }}_{\mathrm{1}}}\mathrm{d}\mathit{x}\mathrm{+}\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}{\mathit{\epsilon }}_{\mathit{i}}\\ & \mathrm{+}& \end{array}$(10)where the ε_i are the actual errors on the relative frequency differences resulting from measurement errors. If the true density and Γ₁ profiles are known beforehand (for instance when testing the methods on stellar models rather than true stars) then it is possible to quantify separately the different contributions to the error using Eq. (10). Of course, in a realistic case, such profiles are not accessible.

Using Eq. (10), one can seek an upper bound on the error through the Cauchy-Schwartz inequality: $\begin{array}{ccc}\left|\frac{\mathit{\delta }\mathit{\rho ̅}\mathrm{inv}\mathrm{-}\mathit{\delta }\mathit{\rho ̅}\mathrm{exact}}{\mathit{\rho ̅}}\right|& \mathrm{\le }& {\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}{\mathit{K}}_{\mathrm{avg}}\mathrm{-}\mathrm{4}\mathit{\pi }\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{R}}}{\mathit{x}}^{\mathrm{2}}\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}_{\mathrm{2}}{\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}\frac{\mathit{\delta \rho }}{\mathit{\rho }}\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}_{\mathrm{2}}\\ & \mathrm{+}& {\parallel {\mathit{K}}_{\mathrm{cross}}\parallel }_{\mathrm{2}}{\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}\frac{\mathit{\delta }{\mathrm{\Gamma }}_{\mathrm{1}}}{{\mathrm{\Gamma }}_{\mathrm{1}}}\begin{array}{}\mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\\ \mathrm{\parallel }\end{array}}_{\mathrm{2}}\\ & \mathrm{+}& \sqrt{\left(\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}^{\mathrm{2}}\right)\left(\sum _{\mathit{i}}{\mathit{\epsilon }}_{\mathit{i}}^{\mathrm{2}}\right)}\\ & \mathrm{+}& \end{array}$(11)where $\mathrm{\parallel }\mathit{f}{\mathrm{\parallel }}_{\mathrm{2}}\mathrm{=}\sqrt{{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathrm{1}}{\mathit{f}}^{\mathrm{2}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}}$. The difficulty, of course, is finding reasonable estimates for  ∥ δρ/ρ ∥ ₂,  ∥ δΓ₁/Γ₁ ∥ ₂, and $\sqrt{{\sum }_{\mathit{i}}{\left[{\mathit{F}}_{\mathrm{surf}}\mathrm{\left(}{\mathit{\nu }}_{\mathit{i}}\mathrm{\right)}\mathit{/}{\mathit{Q}}_{\mathit{i}}\right]}^{\mathrm{2}}}$. In their solar structural inversions, Rabello-Soares et al. (1999) used the quantity  ∥ K_avg − T ∥ ₂ as a measure of how well the averaging kernels matched specified target functions, T (which in their case were localised around specific grid points), and  ∥ K_cross ∥ ₂ to characterise how well cross-talk was being suppressed.

It is also useful to look at the statistical error on the estimated density variation, resulting from frequency measurement error. Let us denote as σ_i the 1σ error bar on the relative frequency differences, δν_i/ν_i. The resulting 1σ error bar on the estimated mean density variation is then given by ${\mathit{\sigma }}_{\mathit{\delta }\mathit{\rho ̅}\mathit{/}\mathit{\rho ̅}}\mathrm{=}\sqrt{\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}\mathit{,}$(12)where we have assumed that the errors on the individual frequencies are independent. The quantity $\sqrt{{\sum }_{\mathit{i}}{\mathit{c}}_{\mathit{i}}^{\mathrm{2}}}$ is known as the error magnification and is the amount by which the error is amplified if the σ_i are uniform.

It is important to realise that there are other sources of error besides those given in Eq. (10). These include non-linear effects where Eq. (1) is no longer valid, non-adiabatic effects which are not included in the pulsation calculations nor in the variational principle (used to calculate the kernels), inaccuracies in the structural kernels which come from neglecting surface terms in the integration by parts (although these may be absorbed up to some extent in  ∑ _ic_iF_surf(ν_i)/Q_i) and numerical inaccuracies.

2.4. Homologous transformations

Before proceding to describe various inversion procedures, it is useful to point out simple properties which can be deduced from homologous transformations. Let us suppose that the observed star can be deduced from the model through a homologous transformation, such that the relative density variation is , where ϵ is a small quantity. As is well known from dimensional analysis, the frequencies scale with $\sqrt{\mathit{GM}\mathit{/}{\mathit{R}}^{\mathrm{3}}}$. Therefore, to first order the relative frequency variation will be ϵ/2. As a consequence, a necessary and sufficient condition to obtain the correct inversion result in the case of a homologous transformation is $\mathrm{2}\mathrm{=}\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}\mathit{,}$(13)where we have made use of Eq. (6). In what follows, inversion procedures which satisfy this condition will be called “unbiased”.

Furthermore, one can also deduce a simple property about the K_ρ,Γ1 kernels. Remembering that δΓ₁/Γ₁ = 0 in a homologous transformation and making use of Eq. (1) yields the following result: $\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\simeq }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}{\mathit{K}}_{\mathit{\rho ,}{\mathrm{\Gamma }}_{\mathrm{1}}}^{\mathit{i}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{x}\mathit{.}$(14)This relation is approximate because in deriving the kernel K_ρ,Γ1, there are several integrations by part in which surface terms are neglected (Gough & Thompson 1991). Nonetheless, Eq. (14) is satisfied to a good degree of accuracy.

Using Eqs. (14), (8) and (5), it is possible to re-express Eq. (13) as follows: ${\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}{\mathit{K}}_{\mathrm{avg}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{d}\mathit{x}\mathrm{\simeq }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{4}\mathit{\pi }{\mathit{x}}^{\mathrm{2}}\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{R}}}\mathrm{d}\mathit{x,}$(15)which is another way of expressing that an inversion procedure is unbiased.

3. Inversion procedures

We will now describe three different methods for estimating the mean density of a star. The first method is an inversion procedure that directly yields inversion coefficients, c_i. The next two methods are non-linear scaling laws rather than inversion procedures. In order to put them in the same form as an inversion procedure, we linearise these laws and obtain associated coefficients c_i. This then allows a systematic comparison of the three methods using the averaging and cross-term kernels. Hence, although the second and third method are not initially inversion procedures, they will be treated and described as such below.

3.1. SOLA inversion

The SOLA procedure (subtractive optimally localised averages, Pijpers & Thompson 1992, 1994) consists in optimising K_avg by minimising its difference with a specified target function, T, and reducing the effects of the remaining terms. This is achieved by minimising the following cost function: $\begin{array}{ccc}\mathit{J}\mathrm{\left(}{\mathit{c}}_{\mathit{i}}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{}\mathrm{1}\\ \mathrm{0}\end{array}{\mathrm{\{}{\mathit{K}}_{\mathrm{avg}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{-}\mathit{T}\mathrm{(}\mathit{x}{\mathrm{)}}^{\mathrm{\}}}}^{\mathrm{2}}\mathrm{d}\mathit{x}\mathrm{+}\mathit{\beta }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}{\left\{{\mathit{K}}_{\mathrm{cross}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\right\}}^{\mathrm{2}}\mathrm{d}\mathit{x}\\ & \mathrm{+}& \sum _{\mathit{m}\mathrm{=}\mathrm{1}}^{{\mathit{M}}_{\mathrm{surf}}}{\mathit{a}}_{\mathit{m}}\sum _{\mathit{i}}\frac{{\mathit{c}}_{\mathit{i}}{\mathit{\psi }}_{\mathit{m}}\mathrm{\left(}{\mathit{\nu }}_{\mathit{i}}\mathrm{\right)}}{{\mathit{Q}}_{\mathit{i}}}\mathrm{+}\mathit{\lambda }\left\{\mathrm{2}\mathrm{-}\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}\right\}\\ & \mathrm{+}& \end{array}$(16)where $\mathrm{⟨}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\frac{\mathrm{1}}{\mathit{N}}{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}$, N being the number of observed frequencies. In order to obtain the relative mean density variation, the appropriate target function is $\mathit{T}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{x}}^{\mathrm{2}}\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{R}}}\mathrm{·}$(17)The first term on the right-hand side of Eq. (16) is used to reduce the differences between K_avg and T, and the next two terms minimise the effects of cross-talk from Γ₁ and pollution from unknown near-surface effects. The functions ψ_m are a basis of slowly varying functions, such as the first M_surf Legendre polynomials, and the a_m are Lagrange multipliers, which are used to suppress these terms. The fourth term ensures that the inversion is not biased, λ being a Lagrange multiplier. Finally, the last term is a regularisation term which reduces (see Eq. (12)). The free parameters in this approach are β, which regulates between optimising K_avg and minimising K_cross, θ, which is used to adjust the amount of regularisation, and M_surf, which controls the number of terms used in suppressing unknown surface effects.

Using Eq. (17) to define the target function avoids inverting for the entire density variation profile, δρ, and integrating it to obtain the mean density variation, which simplifies the calculations and reduces the number of free parameters (such as the choice of a grid on which to carry out the inversion, or the width of the target functions at each grid point). Furthermore, the averaging and cross-term kernels are calculated for the inverted mean density variation, rather than for the density at specific grid points, which is more useful in evaluating the inversion errors. Finally, it is worth noting that Pijpers (1998) had already used a similar direct approach to determine the Sun’s moment of inertia rather than integrating the rotation profile, and had found similar benefits.

3.2. Large frequency separation

An alternative approach is to use the well-known scaling law between the mean density and the large frequency separation (e.g. Ulrich 1986) $\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}\mathrm{\propto }\sqrt{\mathit{\rho ̅}}\mathit{,}$(18)where  ⟨ Δν ⟩  is an average of the large frequency separation, ν_{n + 1,   ℓ} − ν_n,   ℓ. This law has been applied by many authors to obtain the mean density of a variety of stars (e.g. Kjeldsen & Bedding 1995, Mosser et al. 2010). Recently, White et al. (2011) have investigated its accuracy for a grid of models and showed that the error is around 2 to 3.5% for stars with T_eff below 6700 K, relevant to solar-like oscillations, but can reach 10% in the worst cases, when calibrated to solar values.

Equation (18) stems from the fact that the frequencies, and consequently any linear combination thereof, also scale with $\sqrt{\mathit{\rho ̅}}$ for homologous transformations, as was pointed out in Sect. 2.4. This still approximately holds for non-homologous transformations, although some dependence on the precise stellar structure will appear. Asymptotic analysis shows that for high-order p-modes, $\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}\mathrm{=}{\left[\mathrm{2}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{R}}\frac{\mathrm{d}\mathit{r}}{\mathit{c}}\right]}^{-1}$, i.e. the inverse of the time it takes for a sound wave to travel from the surface of a star to its centre and back again (Vandakurov 1967). Given that sound waves travel more slowly near the surface, the large frequency separation is sensitive to surface conditions.

As such, this law is not directly comparable to the above procedure. However, one can apply this law in differential form, i.e. chose a reference model which is close to the star, calibrate the law using this model, and evaluate the mean density of the star. Given that the reference model is close to the star, one can simply perturb Eq. (18) to relate the relative variation of the large frequency separation to the relative mean density variation: $\mathrm{2}\frac{\mathit{\delta }\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}}{\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}}\mathrm{\simeq }\frac{\mathit{\delta }\mathit{\rho ̅}}{\mathit{\rho ̅}}\mathrm{·}$(19)If we assume that the average large frequency separation can be expressed as a linear combination of the frequencies,  ⟨ Δν ⟩  =  ∑ _id_iν_i, the left-hand side of Eq. (19) becomes $\mathrm{2}\frac{\mathit{\delta }\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}}{\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}}\mathrm{=}\frac{\mathrm{2}{\sum }_{\mathit{i}}{\mathit{d}}_{\mathit{i}}\mathit{\delta }{\mathit{\nu }}_{\mathit{i}}}{{\sum }_{\mathit{i}}{\mathit{d}}_{\mathit{i}}{\mathit{\nu }}_{\mathit{i}}}\mathrm{=}\frac{\mathrm{2}{\sum }_{\mathit{i}}{\mathit{d}}_{\mathit{i}}{\mathit{\nu }}_{\mathit{i}}\frac{\mathit{\delta }{\mathit{\nu }}_{\mathit{i}}}{{\mathit{\nu }}_{\mathit{i}}}}{{\sum }_{\mathit{i}}{\mathit{d}}_{\mathit{i}}{\mathit{\nu }}_{\mathit{i}}}\mathrm{=}\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}\frac{\mathit{\delta }{\mathit{\nu }}_{\mathit{i}}}{{\mathit{\nu }}_{\mathit{i}}}\mathit{,}$(20)where ${\mathit{c}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{2}{\mathit{d}}_{\mathit{i}}{\mathit{\nu }}_{\mathit{i}}}{{\sum }_{\mathit{j}}{\mathit{d}}_{\mathit{j}}{\mathit{\nu }}_{\mathit{j}}}\mathrm{·}$(21)Equation (19) then takes on the form given in Eq. (6). Consequently, it is also possible to construct averaging and cross-term kernels for this approach, which can then be compared with those coming from the SOLA approach. By construction, this approach is unbiased and will therefore provide the correct answer for homologous transformations. One can also easily verify that  ∑ _ic_i = 2.

It is important to realise that although the large frequency separation is typically applied to p-modes in the asymptotic regime, Eqs. (18) and (19) have a more general application. Indeed, what is important is the $\sqrt{\mathit{\rho ̅}}$ scaling which applies to all modes, not the fact that the large separation asymptotically gives the inverse sound travel time. Hence, unless specified otherwise, we will typically use all of the available modes when applying Eqs. (18) and (19).

In what follows, we will use the following estimate for the large frequency separation to obtain the coefficients d_i: $\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}\mathrm{=}\frac{{\sum }_{\mathit{\ell }}\mathrm{(}{\mathit{N}}_{\mathit{\ell }}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{⟩}}{{\sum }_{\mathit{\ell }}{\mathit{N}}_{\mathit{\ell }}\mathrm{-}\mathrm{1}}\mathit{,}$(22)where $\begin{array}{ccc}\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{⟩}& \mathrm{=}& \frac{{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{\ell }}}\left[{\mathit{\nu }}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{-}\mathrm{⟨}\mathit{\nu }\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{⟩}\right]\left[{\mathit{n}}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{-}\mathrm{⟨}\mathit{n}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{⟩}\right]}{{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{\ell }}}{\left[{\mathit{n}}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{-}\mathrm{⟨}\mathit{n}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{⟩}\right]}^{\mathrm{2}}}\mathit{,}\\ \mathrm{⟨}\mathit{\nu }\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{⟩}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{N}}_{\mathit{\ell }}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{\ell }}}{\mathit{\nu }}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{,}\\ \mathrm{⟨}\mathit{n}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{⟩}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{N}}_{\mathit{\ell }}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{\ell }}}{\mathit{n}}_{\mathit{i}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathit{.}\end{array}$Here, N_ℓ is the number of modes with a given harmonic degree, ℓ, and the n_i(ℓ) and ν_i(ℓ) are their associated radial orders and frequencies, respectively. The above definition is a weighted average of the least-squares estimates of the mean large frequency separation from Kjeldsen et al. (2008) for each ℓ value. Using a least-squares approach has the advantage of not requiring consecutive radial orders to estimate the large frequency separation.

3.3. The KBCD approach

Recently, Kjeldsen et al. (2008) proposed a method for correcting unknown surface effects which they treated as a power law. This method uses a combination of the large frequency separation and an average of the frequencies to estimate either the exponent of the power law or the mean density of the star. Having established the value of the exponent in the solar case, Kjeldsen et al. (2008) then go on to use this value to estimate the mean density of other stars. In what follows, we will refer to this method as the “KBCD approach”. Like the scaling law based on the large frequency separation, this approach is non-linear and therefore not directly comparable to linear inversion methods. However, if one assumes that the reference model is close to the observed star, it is possible to linearise this approach as was done above. To do so, we start with Eq. (6) of Kjeldsen et al. (2008), subtract 1 from both sides, and retain only first-order terms $\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\delta }\mathit{\rho ̅}}{\mathit{\rho ̅}}\mathrm{\simeq }\frac{\mathit{b}\frac{\mathit{\delta }\mathrm{⟨}\mathit{\nu }\mathrm{⟩}}{\mathrm{⟨}\mathit{\nu }\mathrm{⟩}}\mathrm{-}\frac{\mathit{\delta }\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}}{\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}}}{\mathit{b}\mathrm{-}\mathrm{1}}\mathit{,}$(26)where b is the exponent involved in the power law describing the unknown surface effects, and  ⟨ ν ⟩  an average of the frequencies. Interesting limiting cases of the above equation are

• b = 0: this corresponds to fitting the large frequency separation and brings us back to the previous section;

• b → ∞: this corresponds to fitting the average frequency,  ⟨ ν ⟩ , directly.

One will also note that b = 1 is a singular case which corresponds to a degeneracy in the equations. Kjeldsen et al. (2008) found b = 4.9 for the Sun and consequently used this value for other stars. We will do likewise in what follows.

In practice, we will use Eq. (22) to obtain the large frequency separation, and a similar equation to obtain the average frequency: $\mathrm{⟨}\mathit{\nu }\mathrm{⟩}\mathrm{=}\frac{{\sum }_{\mathit{\ell }}\mathrm{(}{\mathit{N}}_{\mathit{\ell }}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}\mathrm{⟨}\mathit{\nu }\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{⟩}}{{\sum }_{\mathit{\ell }}{\mathit{N}}_{\mathit{\ell }}\mathrm{-}\mathrm{1}}\mathit{,}$(27)thereby extending this method to non-radial modes. Once more, it is straightforward to see that this approach is unbiased.

4. Non-linear extension

The above procedures are linear and are therefore only adapted to reference models which are sufficiently close to the observed star. However, one may try to extend their applicability by pre-scaling the reference model so as to approach the mean density of the star and hopefully bring the relative differences into the linear regime. We therefore introduce a scale factor s, such that the mean density of the scaled reference model is and the corresponding frequencies $\mathit{s}{\mathit{\nu }}_{\mathit{i}}^{\mathrm{ref}}$. With this rescaled model, the frequency shifts then take on the following expression: $\frac{{\mathit{\nu }}_{\mathit{i}}^{\mathrm{obs}}\mathrm{-}\mathit{s}{\mathit{\nu }}_{\mathit{i}}^{\mathrm{ref}}}{\mathit{s}{\mathit{\nu }}_{\mathit{i}}^{\mathrm{ref}}}\mathrm{=}\frac{\mathrm{1}}{\mathit{s}}\left(\frac{\mathit{\delta }{\mathit{\nu }}_{\mathit{i}}}{{\mathit{\nu }}_{\mathit{i}}}\mathrm{+}\mathrm{1}\right)\mathrm{-}\mathrm{1.}$(28)The inverted mean density becomes $\mathit{\rho ̅}\mathrm{inv}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\mathit{\rho ̅}\mathrm{ref}{\mathit{s}}^{\mathrm{2}}\left\{\mathrm{1}\mathrm{+}\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}\left[\frac{\mathrm{1}}{\mathit{s}}\left(\frac{\mathit{\delta }{\mathit{\nu }}_{\mathit{i}}}{{\mathit{\nu }}_{\mathit{i}}}\mathrm{+}\mathrm{1}\right)\mathrm{-}\mathrm{1}\right]\right\}\mathit{.}$(29)If one assumes that the linear inversion procedure is unbiased, then Eq. (29) can be simplified as follows, using Eq. (13): $\mathit{\rho ̅}\mathrm{inv}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\mathit{\rho ̅}\mathrm{ref}\left\{\mathrm{-}{\mathit{s}}^{\mathrm{2}}\mathrm{+}\mathit{s}\left(\mathrm{2}\mathrm{+}\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}\frac{\mathit{\delta }{\mathit{\nu }}_{\mathit{i}}}{{\mathit{\nu }}_{\mathit{i}}}\right)\right\}\mathit{.}$(30)Equation (30) is a second-order equation with the following maximum value: $\mathit{\rho ̅}\max \mathrm{\equiv }\mathit{\rho ̅}\mathrm{inv}\mathrm{\left(}{\mathit{s}}_{\max }\mathrm{\right)}\mathrm{=}\mathit{\rho ̅}\mathrm{ref}{\mathit{s}}_{\max }^{\mathrm{2}}\mathit{,}$(31)where ${\mathit{s}}_{\max }\mathrm{=}\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}\frac{\mathit{\delta }{\mathit{\nu }}_{\mathit{i}}}{{\mathit{\nu }}_{\mathit{i}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}\frac{{\mathit{\nu }}_{\mathit{i}}^{\mathrm{obs}}}{{\mathit{\nu }}_{\mathit{i}}^{\mathrm{ref}}}\mathrm{·}$(32)It is this maximum value which we will use as the best mean density estimate. Indeed, if one were to start with a given reference model, obtain a first mean density estimate through a given linear inversion, scale the reference model using this mean density, re-apply the inversion to obtain a second mean density, re-scale the reference model, and iterate till convergence, then the final mean density would be . Another way of looking at this is that s_max is the scale factor for which the inversion yields no further correction to the mean density: $\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}\left[\frac{\mathrm{1}}{{\mathit{s}}_{\max }}\left(\frac{\mathit{\delta }{\mathit{\nu }}_{\mathit{i}}}{{\mathit{\nu }}_{\mathit{i}}}\mathrm{+}\mathrm{1}\right)\mathrm{-}\mathrm{1}\right]\mathrm{=}\mathrm{0.}$(33)It is then interesting to compare the result from a non-linear procedure to the value obtained from the linearised version of the procedure. If we start with the scaling law between the large frequency separation and the mean density, then the value from the linearised procedure is $\mathit{\rho ̅}\max \mathrm{=}\mathit{\rho ̅}\mathrm{ref}{\left(\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathrm{2}\mathit{\delta }\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}}{\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }\mathrm{⟩}}\right)}^{\mathrm{2}}\mathrm{=}\mathit{\rho ̅}\mathrm{ref}{\left(\frac{\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }{\mathrm{⟩}}_{\mathrm{obs}}}{\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }{\mathrm{⟩}}_{\mathrm{ref}}}\right)}^{\mathrm{2}}\mathit{,}$(34)where we have used Eq. (2). This is identical to what is given directly by the original, non-linearised, scaling law, hence providing additional justification for choosing as the mean density estimate. If we calculate for the linearised version of the KBCD approach, the result is $\mathit{\rho ̅}\max \mathrm{=}\mathit{\rho ̅}\mathrm{ref}{\left(\frac{\mathit{b}\frac{\mathrm{⟨}\mathit{\nu }{\mathrm{⟩}}_{\mathrm{obs}}}{\mathrm{⟨}\mathit{\nu }{\mathrm{⟩}}_{\mathrm{ref}}}\mathrm{-}\frac{\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }{\mathrm{⟩}}_{\mathrm{obs}}}{\mathrm{⟨}\mathrm{\Delta }\mathit{\nu }{\mathrm{⟩}}_{\mathrm{ref}}}}{\mathit{b}\mathrm{-}\mathrm{1}}\right)}^{\mathrm{2}}\mathrm{·}$(35)This time the result is different from what is given by the original non-linear procedure. However, this difference turns out to be minute. Indeed, Eq. (35) corresponds to the result which would be obtained had the unknown surface effects been represented as a power law in terms of the reference frequencies rather than the observed ones. Using this formulation, the exponent obtained in the solar case becomes 4.86 instead of 4.89, and mean density estimates are hardly affected (for b fixed at 4.9). Hence, in what follows, we will use this formulation as representative of the KBCD approach.

Finally, the associated 1σ error bar on the estimated mean density is now given by ${\mathit{\sigma }}_{\mathit{\rho ̅}\max }\mathrm{=}\mathit{\rho ̅}\mathrm{ref}{\mathit{\sigma }}_{{\mathit{s}}_{\max }^{\mathrm{2}}}\mathrm{=}\mathit{\rho ̅}\mathrm{ref}{\mathit{s}}_{\max }\sqrt{\sum _{\mathit{i}}{\mathit{c}}_{\mathit{i}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}\mathit{,}$(36)where we have assumed that σ_i ≪ 1. This is different from what is given in Eq. (12) due to the extra term s_max.

By construction, does not depend on the mean density of the initial reference model as opposed to a linear inversion result (this can also be shown explicitly through Eqs. (31) and (32)). The same also applies to ${\mathit{\sigma }}_{\mathit{\rho ̅}{\max }_{}}$. Furthermore, the various individual kernels, the averaging and cross-term kernels, and the target function remain unchanged since they are non-dimensional and consequently do not depend on the scale factor. As such, they remain a useful diagnostic of the inversion procedure. However, when calculating the error as based on Eq. (10), the relative difference on the density profile needs to take into account the scaling by s_max: ${\left(\frac{\mathit{\delta \rho }}{\mathit{\rho }}\right)}_{\mathrm{scaled}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{s}}_{\max }^{\mathrm{2}}}\left[{\left(\frac{\mathit{\delta \rho }}{\mathit{\rho }}\right)}_{\mathrm{unscaled}}\mathrm{+}\mathrm{1}\right]\mathrm{-}\mathrm{1.}$(37)

5. The solar case

In order to illustrate the above methods, we applied them to the Sun. A set of 104 GOLF frequencies and associated error bars (Lazrek et al. 1997) were used as observations. We only used the frequencies with a signal-to-noise ratio greater than 1, i.e. those which were not in brackets in Table 1 of Lazrek et al. (1997). The n and ℓ values of these modes are given in Table 1, and the average of the error bars on these frequencies is 0.694   μHz. Model S (Christensen-Dalsgaard et al. 1996) was used as the reference model. It is well known that there are surface effects in the Sun which are not accounted for in Model S. This causes the high-order modes to differ in frequency from the observed values, as is illustrated in Fig. 1. For the sake of comparison, the large frequency separation is also indicated in the figure. The pulsation modes were calculated with ADIPLS (Christensen-Dalsgaard 2008), using an isothermal outer boundary condition, and the inversions were carried out using an improved (and corrected) version of InversionKit¹.

Before giving the inversion results, it is useful to discuss what the latest and hopefully most accurate value of the solar mean density is. The mass and radius used in Model S correspond to the values given in Allen (1973). Since then, the solar radius underwent a downward revision (Brown & Christensen-Dalsgaard 1998) following a discrepancy with the radius inferred from solar f-modes (Schou et al. 1997; Antia 1998). Further revisions can be found in Haberreiter et al. (2008), who also point out the need for a generally accepted definition. The solar mass is determined by calculating the ratio GM_⊙/G, where GM_⊙ is known from planetary motions and G, the gravitational constant, from various experiments. The product GM_⊙ is known to a high degree of accuracy so the main source of uncertainty comes from the value of G. A summary of the latest values as well as those used in Model S can be found in Table 2. As can be seen in Table 2 (fifth row), the latest mean density represents a 0.14 ± 0.09% increase over the value used in Model S. However, given that the pulsation frequencies scale with $\sqrt{\mathit{GM}\mathit{/}{\mathit{R}}^{\mathrm{3}}}$, it is not possible to distinguish between a variation on the solar mass and a modification of the gravitational constant from frequency data alone. Hence, it is more appropriate to compare between Model S and the Sun. As can be seen in the last row, the latest value represents a 0.16 ± 0.08% increase over the value from Model S. This corresponds to what one would expect from inversions if these were exact. Given the excellent agreement between GM_⊙ from Model S and the latest value (Table 2, first row), the main contribution to this difference comes from the discrepancy on the radius.

Table 3 lists the results for SOLA inversions with various parameter settings, and for a  ⟨ Δν ⟩  scaling. The corresponding averaging and cross-term kernels are shown in Fig. 2. The first two columns give the values of β and θ, when relevant. No surface corrections are included in the SOLA approach (i.e., M_surf = 0 in Eq. (16)). The third column gives the relative mean density variation (taking into account the non-linear generalisation described in Sect. 4). Evidently, all inversions predict a slight decrease of the mean density, as opposed to the 0.16% increase expected from the latest solar parameters. This is likely to be caused by poorly modelled surface effects. However, we do note that the  ⟨ Δν ⟩  scaling is much more affected, which is not surprising given that the sound travel-time integral is dominated by the near-surface contribution. The next column gives the 1σ error bar around resulting from the measurement uncertainties on δν_i/ν_i. Finally the last two columns give  ∥ K_avg − T ∥ ₂ and  ∥ K_cross ∥ ₂, which intervene in the upper bound on the error, as given in Eq. (11).

As for most inversion methods, the results are a trade-off between minimising the effects of measurement errors and cross-talk from δΓ₁/Γ₁ and optimising the averaging kernel. The parameters β = 10^-8, θ = 10^-4 seem to be a good compromise. If β is increased to 10^-1, the inversion wastes much effort trying to reduce the cross-talk from δΓ₁/Γ₁, which is expected to be small. Although  ∥ K_cross ∥ ₂ is somewhat reduced, the inversion result is worse and has noticeably increased. If instead θ is reduced to 10^-6, the inversion result is slightly better, as is  ∥ K_avg − T ∥ ₂, but both and  ∥ K_cross∥₂ take on high values.

5.1. Near-surface effects

It is then interesting to look at different ways of correcting surface effects. Traditionally, these effects have been modelled as an additive term made up of an unknown function of frequency only, F_surf, normalised by Q_i, the mode inertia divided by the inertia of a radial mode interpolated to the same frequency (Christensen-Dalsgaard & Berthomieu 1991). One way of dealing with this term in a SOLA inversion is to treat the function F_surf as a slowly varying function of ν, such as a linear combination of the first M_surf Legendre polynomials, and then to constrain the inversion to cancel out any such function (Däppen et al. 1991, also see Eq. (16)). The first row of Table 4 gives the inversion result and characteristics for M_surf = 1, which amounts to treating F_surf as a constant. The corresponding averaging and cross-term kernels are given in the top panels of Fig. 3. Obviously, the results are quite poor. A likely reason for this behaviour is that the function F_surf is not satisfactorily described by a constant. Instead, judging from Fig. 1, and remembering that Q_i is normalised in such a way as to be close to 1 for low ℓ values, |F_surf| starts at low values and increases with frequency.

One way of trying to deal with this problem is to normalise the function F_surf in a different way, before treating it as a linear combination of slowly varying functions. Two alternatives are $\frac{{\mathit{G}}_{\mathrm{surf}}\mathrm{\left(}{\mathit{\nu }}_{\mathit{i}}\mathrm{\right)}}{{\mathit{E}}_{\mathit{i}}}\mathit{,} \mathrm{and} \frac{{\mathit{H}}_{\mathrm{surf}}\mathrm{\left(}{\mathit{\nu }}_{\mathit{i}}\mathrm{\right)}}{{\mathit{\nu }}_{\mathit{i}}^{\mathrm{-}\mathit{b}\mathrm{+}\mathrm{1}}{\mathit{Q}}_{\mathit{i}}}\mathit{,}$(38)where we have used the notation G_surf and H_surf to distinguish these functions from the original function F_surf, and where E_i is the mode inertia and b the exponent from the KBCD approach. The first possibility is based on the following reasoning. In the original description, the quantity Q_i was introduced to replace E_i as it does not vary by several orders of magnitude over the p-mode frequency range (Christensen-Dalsgaard & Berthomieu 1991). However, factoring out this strong variation from E_i also means introducing it into F_surf. Therefore, normalising by E_i instead of Q_i avoids this problem. The second possibility is based on the idea that ${\mathit{F}}_{\mathrm{surf}}\mathrm{\left(}{\mathit{\nu }}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{\mathit{Q}}_{\mathit{i}}\mathrm{(}\mathit{\delta }{\mathit{\nu }}_{\mathit{i}}\mathit{/}{\mathit{\nu }}_{\mathit{i}}{\mathrm{)}}_{\mathrm{surf}}\mathrm{\simeq }\mathit{a}{\mathit{\nu }}_{\mathit{i}}^{\mathit{b}\mathrm{-}\mathrm{1}}$, based on Kjeldsen et al. (2008). Hence by dividing both sides by ν^b − 1, one is left with a function H_surf which is slowly varying.

Rows 2–4 of Table 4 give the results from applying these two alternative normalisations, and the corresponding kernels are shown in Fig. 3. In the first two cases, only one (constant) term is used to correct surface effects. The results are slightly better than when no surface corrections are included, with a slight preference for the ν ^− b + 1Q_i normalisation. At the same time, both the averaging kernel and are slightly worse. In the last case, M_surf = 7 in accordance with the optimum number of terms found in Rabello-Soares et al. (1999). However, this leads to poor results. The main difference between the inversions carried out here and those carried out in Rabello-Soares et al. (1999) is the number of modes involved in the inversion. Given the limited number used here (which still remains generous in comparison with stars other than the Sun), one cannot afford to over-correct for surface effects. However, even when no surface corrections are explicitly included, Table 3 shows that a SOLA inversion is less affected by surface effects than a  ⟨ Δν ⟩  scaling, probably because as K_avg approaches the target function T, its surface amplitude is reduced.

Another way to try to limit surface effects is to perform a SOLA inversion using only kernels based on frequency separation ratios. Indeed, Roxburgh & Vorontsov (2003) showed that these ratios are insensitive to surface conditions through analytical developments based on internal phase shifts and through numerical simulations. Likewise, Otí Floranes et al. (2005) showed that the associated kernels have a small amplitude near the surface. However, they also pointed out that both the numerator and denominator of a frequency separation ratio scale with $\sqrt{\mathit{GM}\mathit{/}{\mathit{R}}^{\mathrm{3}}}$, meaning that the ratio is insensitive to the mean density of the underlying model. As a result, the c_i coefficients associated with an inversion based on these kernels add up to 0, making it impossible to perform an unbiased inversion. Likewise, it can be shown using Eq. (14) that the integral of the kernels ${\mathit{K}}_{\mathit{\rho ,}{\mathrm{\Gamma }}_{\mathrm{1}}}^{\mathrm{ratio}}$ is approximately 0.

Finally, the last row of Table 4 gives the results from applying the KBCD approach. The associated averaging and cross-term kernels are displayed in the lower panels of Fig. 3. The mean density estimate and compare favourably with the other inversion results and the associated kernels are well-behaved, thus explaining why this method works. This approach is found to be superior to the simple scaling law based on the large frequency separation, as can be seen both in terms of the results it produces and in terms of the kernels. The likely reason for this is that the averaged frequency, which intervenes in this approach, uses a more uniform weighting on the individual frequencies, whereas the large frequency separation favours frequencies at either end of the frequency range, thus amplifying surface effects, which are strongest at high radial orders. Overall, the best results are obtained with the ν ^− b + 1Q_i normalisation, using one term for surface effects, and the KBCD approach (third and fifth line of Table 4).

6. Test cases with a grid of models

In order to conduct a more systematic study, we applied the different inversion procedures to a grid of 93 main-sequence and pre-main-sequence models. These models were downloaded from grid B of the CoRoT-ESTA-HELAS website², and were chosen within the (T_eff,log (g)) box (5291 ± 100   K,4.563 ± 0.1   dex), the error bars being representative of what can be expected from spectroscopic observations. Their masses range from 0.80   M_⊙ to 0.92   M_⊙ and their ages from 28 Myr to 17.6 Gyr after their appearance on the birthline. A description of the physical ingredients, chemical composition and initial conditions used in these models is given in Marques et al. (2008). In particular, the initial chemical composition is (X,Y,Z) = (0.70,0.28,0.02), where the metal abundances correspond to the solar mixture given in Grevesse & Noels (1993). Figure 4 shows the positions of these models in a (T_eff,log (g)) diagram, along with some evolutionary tracks, the dashed parts corresponding to the pre-main-sequence. Three additional models, nicknamed “Models A, A′ and B”, were selected near the centre of this box to play the role of observed stars. These models are also shown in Fig. 4, and Table 5 gives some of the characteristics of Models A and B. The frequencies of these models are shown in an échelle diagram in Fig. 5. In all cases tested here, the 1σ error bar on the frequencies was set at 0.3   μHz, although no errors were added to the frequencies. The goal of this study is then to see how well the density of these three models can be reproduced by applying the various inversion procedures to the models from the grid.

6.1. Model A

Model A represents a 0.9   M_⊙ main-sequence star and was downloaded from grid A of the CoRoT-ESTA-HELAS website³ Models from grid A are started off as fully convective spheres, instead of using a birthline as in grid B. Apart from this difference, this model has the same chemical composition and physical ingredients as the models in the grid.

Each one of the models from the grid was used as a reference model from which to apply the different inversion techniques to the “observed” frequencies from Model A. Figure 6 shows the results from these inversions for a set of 33 modes with radial orders n = 15−25 and harmonic degrees ℓ = 0−2. The parameter β was set to 10^-8 for all SOLA inversions. Although the mean densities of the reference models differ from that of Model A by up to 30% (see the x-axis in figures), the inversion results are within 0.3 to 0.6% of the true mean density except for the worst cases where the departure can exceed 1.0%, which shows the sensitivity of these modes to density. A more detailed look at the figures shows that overall, the KBCD approach and SOLA inversions produce the best results, provided the parameter θ is not too small. The two are very similar if M_surf = 1, θ = 10^-2 and a ν ^− b + 1Q_nℓ normalisation is used for surface effects. The large frequency separation, for the most part, produces results which are worse. One may also notice that in most cases, the difference between the inverted and true mean densities is much larger than the 1σ error bar from the inversions. This behaviour can be expected, since the 1σ error bar from the inversion is only based on the frequency measurement errors and does not take into account the other sources of error mentioned in Sect. 2.3, such as differences between the averaging kernel and the target function or cross-talk from the δΓ₁/Γ₁ profile.

Evidently, reducing θ mainly seems to increase the 1σ error bar in SOLA inversions without improving the results. This can be understood by looking at the trade-off curves illustrated in Fig. 7. Indeed, reducing θ only has a small effect on  ∥K_avg − T∥ ₂ while increases substantially. Furthermore, the same mechanism which amplifies observational errors can also amplify errors which are not accounted for in Eq. (10), such as non-linear effects not included in Eq. (1) or numerical errors on kernel calculations. Based on Figs. 6 and 7, the value θ = 10^-2 seems to be a good compromise between the two effects.

Furthermore, given that the true structure of Model A is known, it is possible to analyse the different sources of error based on Eq. (10). These are displayed in Fig. 8 for SOLA inversions. The dominant contribution comes from a mismatch between the averaging kernel and the target function and is represented by a solid line. The dotted line represents the cross-talk from the δΓ₁/Γ₁ profile and is the smallest contribution. Lastly, the dashed curve represents remaining effects such as surface contributions, measurement errors on the frequencies and effects not included in Eq. (10), such as non-linear effects not accounted for in Eq. (1). Of these three, the non-linear effects are likely to be the strongest because the data is error-free (except for numerical round-off errors) and no attempt was made to simulate surface effects in Model A. Obviously, the non-linear effects are not negligible compared to the other sources of error. This is not entirely surprising given that the only criteria used for selecting the models from the grid was the (T_eff,log (g)) error box, and consequently shows the limits of Eq. (1). It must also be noted that the errors in Fig. 8 tend to compensate for each other. However, one cannot always expect such a fortuitous cancellation to arise in other situations.

Figure 9 shows what happens if the non-linear extension described in Sect. 4 is not applied to the inversion procedures. As expected, inversion results are substantially worse when the mean density of the reference models is far from the true mean density.

Finally, the wiggles observed in Fig. 6, especially for reference models with a lower mean densities, can be partially explained by the fact that the grid contains a mixture of pre-main-sequence and main-sequence models. Indeed, when the two types of models are separated as in Fig. 10, two separate trends can be observed which likely reflect their differences in structure.

6.2. Near-surface effects (Model A′)

Model A′ is nearly identical to Model A, except for an ad hoc 50% increase of the density in the outer 1% of the model, the transition being represented by a hyperbolic tangent function. The pressure was then recalculated using the hydrostatic equation, thus resulting in a very slight offset from the original profile. The difference in mean density between the two models is negligible, but the frequencies are noticeably different (see Fig. 5) because sound waves spend more time in the outer layers of a star.

Figure 11 shows the inversion results for this model, using the same set of modes as in the previous case. Overall, the results are worse than in the previous case. Furthermore, the scaling law based on the large frequency separation yields results which are substantially offset from the true mean density, when compared to the other procedures. This could be expected because the large frequency separation, which asymptotically gives the travel time of a sound wave from centre to surface, is very sensitive to surface conditions. This high sensitivity to surface conditions is also reflected in the associated averaging kernels which tend to have large amplitudes near the surface. The SOLA inversions without surface corrections also seem to be affected, but to a lesser degree. In particular, it produces worse results than the KBCD approach for the higher density models, which is the opposite from the previous case. If surface corrections are included in the SOLA approach, they produce slightly better results than the KBCD approach. Once more, reducing θ does not improve the SOLA inversions, but only serves to increase the 1σ error bar.

6.3. Uncertainties in the stellar physics (Model B)

Model B is different in a number of ways from Model A and also from the models in the grid. It has a different chemical composition and includes the effects of diffusion and rotational mixing, thereby providing a way to test the effects of unknown physics on mean density inversions. Specifically, the chemical abundances are (X,Y,Z) = (0.735,0.250,0.015) and follow the more recent solar mixture given in Asplund et al. (2005).

The results are shown in Fig. 12, using once more the same set of modes. These results are indeed similar to those obtained for Model A and show that all inversion procedures presented here seem to be fairly robust to unknown physics.

6.4. The effects of mode misidentification

We then investigated the effects of misidentifying the radial orders of the modes. Figure 13 (left panels) shows mean density estimates for Model A, using the same set of observed modes as previously, but identified as n = 16−26 (upper panel) and n = 14−24 (lower panel). As can be seen in these figures, the scaling law based on the large frequency separation provides the best results, whereas the other methods are off by 10%. A simple explanation is that the relative change over one radial order of Δν is smaller than that of ν, given that Δν goes to a constant when n goes to infinity, whereas ν is roughly proportional to n.

Another way of misidentifying modes is by labelling the ℓ = 0 modes as ℓ = 1 and vice-versa (e.g. Benomar et al. 2009a). This, however, also leads to detecting a different set of ℓ = 2 frequencies because these will always be close to the ℓ = 0 frequencies. In order to investigate this scenario, we swapped the ℓ = 0 and 1 identifications and created a fictitious set of ℓ = 2 frequencies in such a way as to preserve the original small frequency separation, averaged over adjacent radial orders. Two possible ways of assigning radial orders exist and are given in Table 6. The corresponding mean density estimates are given in the right-hand panels of Fig. 13. Once more, the large frequency separation gives the best results, the other methods being off by about 5%, i.e. by about half the amount compared to the previous case when only n was misidentified. One way of interpreting this is to view this as a  ± 1/2 offset on the radial order, and then to apply the same argument as above.

These observations provide a way of identifying modes: one can calculate the mean density using the large frequency separation and another method such as the KBCD approach or a SOLA inversion and search for the values of ℓ and n, for which the two methods agree to within 1% or so. This procedure, although different from the methods described in Bedding & Kjeldsen (2010) and White et al. (2011), is based on the same principles. Indeed, the method presented here and the approach from Bedding & Kjeldsen (2010) both amount to scaling a model (or a comparison star) so that the large frequency separation agrees with that of the observed star, then adjusting the identification so that the offset between the two sets of frequencies is minimised. White et al. (2011) achieve the same result by directly comparing their model-based values for the large frequency separation and ϵ, the additive offset from the asymptotic formula for p-mode frequencies, to those from observed stars. Of course, surface effects acting upon the large frequency separation may lead to an offset even for the correct identification, although we expect such an offset to be smaller than what is shown in Fig. 11, which was calculated for a fairly drastic modification of the surface. White et al. (2011) have nonetheless pointed out situations where ϵ can vary rapidly, which makes mode identification more difficult.

6.5. Kernel mismatch

It it also interesting to test the various inversion techniques on different sets of modes. We therefore considered the following 2 sets:

• Set 1: n = 2−10 for ℓ = 0,2,3, n = 5−10 for ℓ = 1, n = 1−10 for ℓ = 4.

• Set 2: the same as Set 1 without the n = 1, ℓ = 4 mode.

The mean density estimates are shown in Fig. 14. As can be seen by comparing the two figures, including the n = 1, ℓ = 4 mode has a drastic effect on the SOLA inversion result for one of the models, which we shall call the “worst model”. The explanation for this is quite simple. Figure 15 shows the kernels associated with this mode for the worst model and for Model A. Clearly, the kernels in the worst model are dominated by a g-mode behaviour (even if the eigenfunctions still seem to be dominated by a p-mode character), whereas those of Model A are typical of p-modes. This difference on the kernels inevitably leads to a large discrepancy between the actual frequency shift for that particular mode and the frequency shift which would be deduced by applying the linear approximation given in Eq. (1). This, of course, affects any mean density estimates, particularly a SOLA inversion which explicitly optimises the averaging and cross-term kernels using the individual kernels.

7. Observed stars

7.1. α Cen B

Having tested the inversion techniques on stellar models in a number of ways, we now try them out on observed stars. A particularly interesting example is the nearby K1 main-sequence star α Cen B. Its close proximity allows precise interferometric measurements of its angular radius (Kervella et al. 2003). When combined with a reprocessed Hipparchos parallax (Söderhjelm 1999), this yields a radius R = 0.863 ± 0.005   R_⊙. Given that α Cen B is part of a multiple system, it is possible to determine its mass through orbital parameters. Hence, Pourbaix et al. (2002) obtained a mass of M = 0.934 ± 0.0061   M_⊙ using the parallax from Söderhjelm (1999). This then yields a mean density of 2.046 ± 0.058 g   cm^-3. A time series of radial velocities was obtained by Kjeldsen et al. (2005) using the Very Large Telescope (VLT) and the Anglo-Australian Telescope (AAT). A Fourier analysis of this time series yielded a set of 37 frequencies with harmonic degrees ℓ = 0 to 3, and radial orders, n, between 17 and 32.

Given that the mass of α Cen B is close to the mass range covered by the grid from the previous section, we simply reused the models from this grid as reference models from which to estimate the mean density of α Cen B. The results are shown in Fig. 16. These turn out to be similar to the results obtained in the previous sections. The tendencies are the same, as is the spread in inverted mean densities, except for the results from the  ⟨ Δν ⟩  scaling where the spread is somewhat larger. Nonetheless, a noticeable difference appears between the estimates from the  ⟨ Δν ⟩  scaling and the two other methods. This suggests that surface effects play an important role in this star.

It is then interesting to see how these densities compare with previous studies. Table 7 gives a summary of previous and current estimates of the mean density of α Cen B. As was mentioned above, the mass from Pourbaix et al. (2002) and the radius from Kervella et al. (2003) lead to a mean density of 2.046 ± 0.058 g   cm^-3. The results from the KBCD and SOLA approaches are 2 to 3% lower than this value, but remain within the error bars. They agree with the range of values, delimited by the triple-dot-dashed horizontal lines in Fig. 16, obtained in Kjeldsen et al. (2008) where the KBCD approach is developed. The seismic study in Eggenberger et al. (2004) leads to a mean density of 1.998 g cm^-3, which also falls within this range, although they did increase the error bar on the radius to reach this value (the error bar on the mass seems to have been ignored in their study). On the other hand, the  ⟨ Δν ⟩  scaling leads to values below all of these previous studies and which mostly lie outside the error bars of the non-seismic mean density. Hence, using the KBCD or SOLA approaches leads to better agreement with the non-seismic constraints and tends to confirm the idea that the  ⟨ Δν ⟩  scaling is affected by surface effects.

7.2. HD 49933

We now turn our attention to stars observed by CoRoT. A particularly interesting example is the F5 main-sequence star HD 49933, which has been observed in two separate CoRoT runs, thereby yielding a rich frequency spectrum. In what follows, we will use the set of 51 frequencies and associated error bars and identifications given in Benomar et al. (2009b), which takes into account the two observational runs. An alternative identification had been favoured before data from the second observational run were available (e.g., Appourchaux et al. 2008). The grid of models used in the previous sections is not suitable for this star. We therefore downloaded a set of 42 pre-main-sequence and main-sequence models from the CoRoT-ESTA-HELAS website, which lie in the error box T_eff = 6665 ± 215 K and log (L/L_⊙) = 0.54 ± 0.02 dex, based on the values quoted in Benomar et al. (2010) and references therein. No attention was paid to possible differences in chemical composition between these models and HD 49933.

The mean density estimates are given in Fig. 17. They fall in the interval 0.580 ± 0.015 g   cm^-3, which corresponds to a  ± 2.6% range around the mid-value. This interval is on the lower end of the range of values covered by the reference (or input) mean densities. Furthermore, the estimated mean density slightly increases with the reference mean density, in much the same way as in the previous examples but on a larger scale. If we limit ourselves to reference models with a closer mean density to the results, this would favour the lower values, but this reasoning must be taken with some caution since the inversions use the non-linear extension and are therefore insensitive to the mean density of the input models. Various other groups have performed detailed asteroseismic interpretations of his star and have come up with best-fitting models using a variety of techniques. Here, we will only focus on studies for which both the mass and radius of the best-fitting model are provided. Kallinger et al. (2010) searched for a best-fitting model in a dense grid using a χ² criteria, whereas Creevey & Bazot (2011) found two solutions using both a Levenberg-Marquardt and a Markov Chain Monte Carlo approach. Their results are summarised in Table 8 and represented as triple-dot-dashed horizontal lines in Fig. 17. The first solution uses a mode identification in which the ℓ values are swapped. Its mean density lies in the upper half of the range of values found here. The two other solutions use the same identification as the one used here. They have very similar mean densities in spite of having different radii and masses, and are at the lower end of the interval of solutions found here. Although deviating more from our mid-value, these solutions are located where the mean densities are least modified by the inversions.

Given the difficulties in labelling the modes prior to the second CoRoT run, we tested the swapped ℓ identification, using only the ℓ = 0 and 1 modes. The results are shown in Fig. 18 and are very similar to the results presented in Fig. 13. Clearly, the models used in this study favour the currently accepted identification. In particular, this agrees with the identification favoured by Kallinger et al. (2010) and with the results from Bedding & Kjeldsen (2010), who compared HD 49933 to other observed stars whose identification is known.

Finally, it is worth noting that surface effects have not manifested themselves in the results presented here: all three methods yield similar results for the correct identification. This suggests that surface effects are less important in this star compared to other effects deeper within the star. Interestingly, Bedding et al. (2010) recently reported on an asteroseismic investigation of Procyon, also of spectral type F5. They showed that the KBCD approach with a solar exponent, b = 4.9, did not provide a good description of the differences between the frequencies of various best-fitting models and those of the star. Instead, these were better described as a constant offset, which corresponds to b = 0. Because these differences do not go to zero when the frequency is reduced, the authors also interpreted this as signifying that the structural differences between their best-fitting models and Procyon extend beyond the surface layers.

7.3. HD 49385

The star HD 49385 is an evolved G0 type star which was observed by CoRoT for 137 days during its second long run. Deheuvels et al. (2010) were able to detect 27 ℓ = 0−2 frequencies as well as what seems to be mixed modes and some unconfirmed ℓ = 3 frequencies. In what follows, we will not use the ℓ = 3 modes because of their low signal-to-noise ratio, and will postpone dealing with the mixed modes to the end of this section, because they vary rapidly from one model to another, thereby limiting the validity of Eq. (1). A new set of 18 models within the (T_eff,log (g)) = (6095 ± 65   K,   4.00 ± 0.06) error box (Deheuvels et al. 2010) was downloaded from the CoRoT-ESTA-HELAS website. Given the position of this box in the HR diagram, these models correspond to either pre- or post-main-sequence stars⁴.

The mean density estimates are given in Fig. 19. The horizontal triple-dot-dashed lines are the mean densities corresponding to the four families of solutions in Deheuvels & Michel (2011). In the upper panel, the results are particularly poor. This failure is caused by differences in mode labelling between the pre- and post-main-sequence models owing to the presence of mixed modes in the latter. The inversion results were based on an identification which corresponds to some of the best-fitting models from Deheuvels & Michel (2011), who applied Newton’s method to a grid of models specifically designed to handle mixed modes. Such an identification is ill-adapted for pre-main-sequence models and to a lesser extent for some of the post-main-sequence models because of the jumps it introduces in the n values. Consequently, the results from the five pre-main-sequence models, indicated by the circles, stand out as being the worst. In order to circumvent these difficulties, a second inversion was carried out, only using the ℓ = 0 modes. This time, the results are far better and are within 2.5% of the mean densities from Deheuvels & Michel (2011) for the SOLA inversions and for the KBCD approach, and 3.5% for the scaling law based on the large frequency separation. Paradoxically, the pre-main-sequence models now give the results which are closest to Deheuvels & Michel (2011). Once more, the scaling law is somewhat different from the two other approaches, which is probably indicative of surface effects.

In order to investigate the influence of mixed modes on SOLA mean density inversions, we used a low-overshoot model from Deheuvels & Michel (2011) with the Asplund et al. (2005) mixture. Inversions were carried out for three different sets of modes:

• Set 1: 26 p-modes with ℓ = 0 − 2.

• Set 2: Set 1 + the (n = 14,ℓ = 1) p-mode.

• Set 3: Set 1 + the (n = 10,ℓ = 1) mixed mode⁵.

Two different sets of SOLA parameters were used:

• Case 1:β = 10^-8, θ = 10^-2 and M_surf = 1. These parameters were chosen in accordance with typical inversion parameters (the number of modes was sufficient to allow M_surf = 1).

• Case 2: β = 0, θ = 0, M_surf = 0. These parameters correspond to optimising K_avg only.

The results are summarised in Table 9. As can be seen in case 1, adding a p-mode or a mixed mode has a similar impact on both the averaging kernel and the cross-term kernel. The reduction of the statistical error on the inverted mean density, , is higher for the mixed mode than for the p-mode, indicating that the mixed mode may be yielding more information on the mean density. This seems to be confirmed in case 2, where only the averaging kernel is being optimised. Indeed, the averaging kernel including the mixed mode is closer to the target function than the averaging kernel which uses the extra p-mode. However, the extra benefit from the mixed mode remains small and is probably of little use because the domain for which Eq. (1) is valid is likely to be much smaller for these modes than for ordinary p-modes. In some cases, non-linear effects on these modes can even cause a substantial error on the mean density, as was shown in Sect. 6.5.

Does this mean that mixed modes are therefore not a useful diagnostic of core conditions, and that previous studies somehow exaggerated their importance? Certainly not. Rather, it means that the methods presented here are not as efficient at exploiting the full diagnostic potential of these modes. Other methods, such as that presented in Deheuvels & Michel (2011), or inversions localised around the stellar core, are needed to fully benefit from these modes. A possible explanation for this failure can be described as follows. Structural kernels associated with mixed modes and g-modes, such as those illustrated in the upper panels of Fig. 15, show a strong oscillatory behaviour near the centre of the star, which does not match the target function, T, defined Eq. (17). As a result, their contribution to SOLA inversions tends to be reduced, as reflected, to some extent, in their inversion coefficients. On the other hand, their contribution is not reduced in the KBCD approach or the  ⟨ Δν ⟩  scaling law, and this consequently leads to averaging kernels with highly oscillatory behaviour near the centre, which can also be a source of error.

8. Discussion and conclusion

We have investigated different techniques for estimating directly the mean density of a star. These include the usual scaling law based on the large frequency separation, an approach described in Kjeldsen et al. (2008) (the “KBCD approach”), and a SOLA inversion technique with a non-linear extension. We have shown how all three techniques can be put in a similar form, which allows the construction of associated averaging and cross-term kernels, useful for characterising what is actually being measured and for quantifying the error.

A comparison of the three techniques, applied to the Sun and to various test cases, shows that the KBCD approach and SOLA inversions yield similar results, whereas the scaling law based on the large frequency separation produces sub-optimal results. The averaging kernels produced by the KBCD approach can come close in quality to those produced by a SOLA inversion and are substantially better than the averaging kernels associated with the simple scaling law based on  ⟨Δν⟩, thereby explaining its relative success. Hence, if one only aims to estimate a mean density and does not need a full diagnostic on the accuracy of the result, the KBCD approach should be sufficient. However, for a more detailed analysis, the SOLA method seems ideal because it is specifically designed to optimise the averaging kernel while reducing various sources of error. Furthermore, possible improvements to the SOLA method, such as introducing weight functions in the cost function Eq. (16) (e.g. Rabello-Soares et al. 1999), may help to fine-tune this method and should therefore be explored.

As expected from asymptotic theory, near-surface effects are shown to have a stronger impact on the  ⟨Δν⟩  scaling law than on the two other methods, even when surface corrections are not included in the SOLA approach. However, these two methods require a correct identification of the modes, whereas the  ⟨ Δν ⟩  scaling law is much more tolerant of mode misidentification. Hence, looking for an ℓ and n assignment which provides the best agreement between the simple scaling law and one of the two other methods can be a useful way of identifying modes. This approach to mode identification, although different from the methods presented in Bedding & Kjeldsen (2010) and White et al. (2011), is also based on the same principles.

Quantitatively, a precision around 0.5% or better was reached for the Sun, for α Cen B (as based on the spread of the results), and for various test cases using SOLA inversions (with β = 10^-8 and θ = 10^-2) and the KBCD approach, whereas the scaling law was less accurate with a precision of 1% or worse. However, when applied to two stars observed by CoRoT, the spread in values was larger, up to  ±2.5% for the first two methods, and  ± 3.5% for the  ⟨Δν⟩  scaling law. The reason for this remains an open question, but does not seem to stem from the fact that these are observed stars, as indicated by additional tests involving models instead (results not included) and the results obtained for α Cen B. Rather, it may be that in the part of the HR diagram which contains these stars, there are larger variations in the density profile with stellar evolution. A  ± 3.5% spread is comparable to the spread obtained in White et al. (2011) for stars with T_eff below 6700 K. Above this temperature, White et al. (2011) obtained larger differences, which in some cases reached 10%, between the true density and that obtained from the  ⟨Δν⟩  scaling law. This raises the question as to how accurate the KBCD and SOLA approaches would be in a similar situation, although one also has to bear in mind that this temperature range is less relevant for solar-like oscillations. One notable difference between the results described here and those presented in White et al. (2011) is that the latter are calibrated to solar values rather than nearby reference models, which is likely to be an additional source of error, especially for stars with a large temperature difference with respect to the Sun.

If combined with a 2% error bar on measured radii, as is expected from Gaia, a 2.5% to 3.5% uncertainty on the mean density yields an 8.5−9.5% error bar on the mass of the star, where the dominant contribution to the error comes from the uncertainty on the radius. Such an error bar will then directly translate itself into a comparable uncertainty on the inferred mass of any potential exoplanet (e.g. Guillot & Havel 2011, and references therein). Accordingly, one may try to obtain more accurate results through more complicated approaches such as searching for best-fitting models in a grid or applying non-linear inversion techniques. These approaches will overcome the non-linear effects shown in Fig. 8 and which are beyond the control of the methods presented here. However, given the difficulties in optimising the averaging kernel beyond a certain point for a given set of modes, one is still left with the question as to how much information is intrinsically available in the pulsation frequencies and how much is being assumed in these more sophisticated methods.

Finally, it was shown that the methods presented here do not fully exploit the diagnostic potential of mixed modes. Indeed, mixed modes have kernels with a highly oscillatory behaviour near the centre, which strongly differs from the target function, T, defined in Eq. (17). This reduces the contribution of mixed modes to SOLA inversions and leads to averaging kernels with an oscillatory behaviour near the centre for the  ⟨ Δν ⟩  scaling law and the KBCD approach. Also, the linear domain for these modes is smaller and can lead to substantial errors when outside this domain, as shown in Sect. 6.5. However, one may expect these modes to be very useful when applying localised inversion methods to models which are close to the true stellar structure.

---

Acknowledgments

We thank the referee for helpful comments and suggestions which have improved the manuscript. We also thank I. Roxburgh, S. Vorontsov and S. Basu for interesting discussions and answers to questions. D.R.R. is supported by the CNES (“Centre National d’Études Spatiales”) through a postdoctoral fellowship. This article made use of InversionKit, an inversion software developed in the context of the European Helio- and Asteroseismology Network (HELAS), which is funded by the European Commission’s Sixth Framework Programme. NCAR is sponsored by the National Science Foundation. D.R.R. also acknowledges financial support through a postdoctoral fellowship from the “Subside fédéral pour la recherche 2011”, University of Liège.

References

All Tables

All Figures