© ESO, 2015

1. Introduction

Surface layer independent model fitting seeks to find a stellar model (or models) whose interior structure is such that the oscillation properties of the model are compatible with surface layer independent combinations of an observed frequency set ν_n,ℓ with error estimates σ_nℓ.

The reason for developing such techniques is that modelling the outer layers of a star is subject to considerable uncertainties (cf. Christensen-Dalsgaard et al. 1988; Dziembowski et al. 1988) due to our poor understanding of the physical processes that determine the structure of these layers; these include modelling convection, convective overshooting, non-adiabatic effects on both convection and oscillations, turbulent pressure, the equation of state, diffusion, mild turbulence, magnetic fields, rotation, global circulation... All these factors impact on the oscillation frequencies of a model star, and therefore hinder efforts to find stellar models whose eigenfrequencies best fit an observed frequency set. Even for the solar case observations and models disagree by up to ~10 μHz (Christensen-Dalsgaard et al. 1996).

We here present details of a surface layer independent technique based on phase matching. The method draws on the result that the contribution of the outer layers of a star to its frequencies is, to a very good approximation, independent of the angular degree ℓ. In this respect it is similar to the method using the ratio of small to large separations (cf. Roxburgh & Vorontsov 2003a; 2013) but, unlike the separation ratio technique, it is not dependent on the actual value of the radial order n, the errors as a function of frequency are uncorrelated and, in principle, it can be applied to stars which have mixed modes.

In the outer layers of a star the phase shift α_ℓ(ν) of the eigenfunction of an oscillation mode is, to high accuracy, independent of angular degree ℓ, allowing us to replace α_ℓ(ν) by a degree independent function α(ν), which is determined solely by the structure of the outer layers, whereas in the inner layers the phase shifts δ_ℓ(ν) depend on ℓ, and in the outer layers of a star are determined solely by the structure of the inner layers. So for the model to fit the frequencies a function (the phase function) of the inner phase shifts, calculated at the observed frequencies, must be such that it matches on to an ℓ independent function only of frequency in the outer layers. The extent to which this is satisfied is a measure of the goodness of fit of the model to the frequencies.

A different approach to overcoming the uncertainty in our modelling of the outer layers of a star is the “frequency offset technique” (Kjeldsen et al. 2008), in which the difference between the observed solar frequencies and those of a “best solar model” is fitted by a power law aν^b, which is then scaled by a single factor (determined by observed and model average frequencies and large separations) and applied to models prior to seeking a best fit model to observed frequency sets of other stars. This technique is widely used and has been incorporated into the Asteroseismic Modeling Portal (AMP) software (Metcalfe et al. 2009). However the assumption that the many differences in the properties of the outer layers of different stars (mass, luminosity, radius, composition...) can be captured in a single scaling factor remains to be verified.

The phase matching technique was briefly outlined in Roxburgh & Vorontsov (2003b), Roxburgh (2010), but in Sect. 2 we present a more detailed description of the theory underlying the procedure and define the phase function with examples of both main sequence models and a post main sequence model with mixed modes. In Sect. 3 we give details of the phase match algorithm and the application to model stars. In Sect. 4 we apply the procedure to data on two stars: HD 49933 observed by CoRoT, and HD 177153 (aka Perky) observed by Kepler, comparing our results with previous investigations. In Sect. 5 we examine in detail the validity of the approximations made in the analysis and in Sect. 6 we compare the method with other surface independent model fitting procedures using separation ratios.

2. Theory of model fitting by phase matching

2.1. The phase shifts

For any angular degree ℓ the frequencies ν of a spherical (or nearly spherical) model star are given as the eigenvalues of the 4th order system of equations governing the oscillations, subject to the boundary conditions of regularity at the centre, and surface conditions on the gravitational potential and pressure perturbations (Pekeris 1938).

However for radial (ℓ = 0) modes the equations can be reduced to 2nd order (Eddington 1918) which we write in terms of the independent variable ψ₀(ν,t) = rp′/ (ρc)^1/2(1)is the acoustic radius, p′ the Eulerian pressure perturbation and Q₀(νt) and V₀(ν,t) are acoustic potentials that depend on both frequency and the structure variables: density ρ(r), pressure P(r), sound speed c(r) (c² = Γ₁P/ρ) and adiabatic exponent Γ₁(r). Were Q₀ and V₀ = 0, this is just the simple harmonic equation whose solution is sin(2πνt). The gravitational potential is calculated from Poisson’s equation using the solution of Eq. (1) and satisfies the surface boundary condition.

Since Eq. (1) is a 2nd order homogeneous equation it can be reduced to a 1st order inhomogeneous equation in the variable 2πνψ₀/ (dψ₀/ dt) and solved in terms of an inner phase shift δ₀(ν,t) defined by (2)where δ₀(νt) satisfies the equation (3)The central boundary condition of regularity at t = 0 requires δ₀ = 0 at t = 0; given the structure of the star δ₀(ν,t) can then be evaluated at any t for any ν (not necessarily an eigenvalue), and the value at any t_f is determined solely by the structure interior to t_f. For general ν the solution does not satisfy the surface boundary condition on ψ₀, the requirement that it do so determines the eigenfrequencies.

One can equally write Eq. (1) in terms of the acoustic depth τ and represent the solution in the form (4)where α₀ is a continuous function of ν and τ which satisfies the surface boundary condition on the pressure perturbation and is determined solely by the structure of the outer layers above τ.

For an eigenfrequency 2πνψ₀/ (dψ₀/ dt) must be continuous at all t, hence on equating the expressions in Eqs. (2) and (4) at any intermediate t_f, and recalling that if tanA + tanB = 0 then A + B = nπ for integer n, we obtain the result that an eigenfrequency ν_n0 satisfies (5)where is the total acoustic radius.

For general ℓ we solve the full 4th order equations for eigenfrequencies ν and Eulerian pressure eigenfunctions , set and define inner and outer phase shifts δ_ℓ(ν,t),α_ℓ(ν,τ) by The inclusion of the factor ℓπ/ 2 in Eq. (6a) reflects the fact that the eigenfunctions behave like spherical Bessel functions near t = 0 and leads to the approximatey equidistant separation of ℓ = 0,1 modes. Since for an eigensolution ψ_ℓ, dψ_ℓ/ dt are continuous at any t_f, τ_f, it follows that for eigenfrequencies ν = ν_nℓ(7)This is the eigenfrequency equation of Roxburgh & Vorontsov (2000; 2003a).

Note that the integer n is just the integer value of the left hand side of Eq. (7) and, since α_ℓ and δ_ℓ can only be determined to within an arbitrary multiple of π, its value depends on the particular algorithm used in determining the α and δ that satisfy the Eqs. (6a) and (6b). The algorithm used in the following analysis is

where j_a is the number of nodes of ψ_ℓ(τ) counting in from the surface to τ_f with j_a = 0 at the first node, and j_d is the number of nodes in ψ_ℓ(t) counting out to t_f with j_d = 0 at the centre where we set δ_ℓ(ν,0) = ℓπ/ 2.

With this definition the integer n = j_d + j_a + 1 is the number of nodes in the eigenfunction ψ_ℓ, but this is not necessarily the same as the n values assigned to an observed frequency set, since the number of nodes in the eigenfunction is not an observable.

2.2. A main sequence example

To illustrate the analysis we consider one model (Model A) out of an evolutionary sequence of models of a star of mass 1.15 M_⊙ with an initial hydrogen abundance X = 0.72, heavy element abundance Z = 0.015, GN93 relative abundances (Grevesse & Noels 1993), OPAL opacities (Iglesias & Rogers 1996) supplemented by Wichita opacities at low temperatures (Ferguson et al. 2005), NACRE reaction rates (Angulo et al. 1999), EOS5 equation of state (Rogers & Nayfonov 2002), mixing length parameter α_c = 1.8 and an Eddington atmosphere out to an optical depth of 10^-3. The models were calculated using the STAROX code (Roxburgh 2008) and do not include diffusion nor convective overshooting.

Model A is a main sequence model evolved to a central hydrogen abundance X_c = 0.246, radius 1.231 R_⊙ and effective temperature T_eff = 6339 K, with a frequency set ν_kℓ of 14 frequencies for each ℓ = 0,1,2,3 centred around ν_max = 2207 μHz given by (cf. Brown et al. 1991; Kjeldsen et al. 1995) (9)with ν_⊙ = 3050 μHz and T_⊙ = 5777 K; The frequencies are listed in Table 1 and illustrated in the echelle diagram Fig. 1; here k = 1,14 is just a label to distinguish the frequencies, and n is the number of nodes in the eigenfunction.

The number of nodes in the eigenfunction is only known for models and for an observed frequency set we have to assign n-values based on observables. For main sequence stars the frequencies show a regular pattern with distinct vertical bands for ℓ = 0,1,2,3 in the echelle diagram as in Fig. 1, and a repeated ordering of the frequencies ℓ = 2,0,3,1 in the power spectrum which allows us to assign observational n-values. We define the mean large separation Δ in terms of four ℓ = 0 frequencies around ν_max as (10)which for Model A gives Δ = 106.54 μHz, and assign an n values to one ℓ = 0 frequency close to ν_max by the condition 0.5 ≤ ν_nℓ/ Δ < 1.5 and then assign n-values to the other frequencies in the order { n − 1,n,n − 1,n } for modes of degree ℓ = 2,0,3,1. The n-values assigned in this way all satisfy the condition (11)For most main sequence models the n values assigned in this way are the same as the number of nodes in the eigenfunctions ψ_ℓ, but we cannot assume this to be true, and it is not the case for the post-main sequence model considered below, where different frequencies with the same ℓ values can have the same n values as assigned by Eq. (11) and/or the same number of nodes.

In Fig. 2 we show the phase shifts for Model A evaluated at a fractional radius x_f = r/R = 0.95 using the definitions in Eqs. (8a), (8b). The δ_nℓ depend on both degree ℓ and radial order n, whereas the α_nℓ almost all lie on the same curve α(ν) independent of degree ℓ, the maximum departure δα_ℓ/α(ν) ~ 4 × 10^-4 (see Sect. 5). It is this almost ℓ independence of α_ℓ(ν) that is the basis of the phase matching technique (and of the use of separation ratios). The solid lines in this figure are the continuous phase shifts calculated as described below.

To calculate continuous inner phase shifts δ_ℓ(ν,t) one solves the oscillation equations for any frequency ν subject to the central boundary conditions and the surface gravitational boundary condition, but not the surface pressure condition, determining the variable ψ_ℓ(ν,t) for this ν. The inner phase shift δ_ℓ(ν,t) is then given by Eq. (6a) with this value of ν and ψ_ℓ. In practice it is sufficient to impose the gravitational boundary condition anywhere in the outer layers where the density is very small so its contribution to the potential and hence to δ can be neglected.

Continuous outer phase shifts α_ℓ(ν,τ) can be determined for any ν by the (very good) approximation of solving the oscillation equations from the surface inwards subject to the pressure boundary condition and setting the oscillating gravitational potential and its derivative to zero at the surface¹. These boundary conditions completely determine the solution (which of course will not, in general, satisfy the central boundary conditions). The solution gives ψ_ℓ(ν,τ) and the outer phase shift α_ℓ(ν,τ) is then given by Eq. (6b) with this value of ν and ψ_ℓ.

In Fig. 3 we show the continuous phase shifts for Model A, calculated at a frequency ν = 2133.16 μHz, plotted against both fractional acoustic radius t/T, and fractional radius x = r/R. The inner phase shifts δ_ℓ depend on ℓ and are primarily determined by the inner structure, becoming almost constant in the intermediate layers. In the outer layers the surface phase shifts α_ℓ(ν,t) are almost independent of degree ℓ and away from the surface region are almost constant in the intermediate layers. This behaviour of the phase shifts holds in general for other models and is discussed in more detail in Sect. 5 below.

The frequency ν is the ℓ = 2,n = 18 eigenfrequency so we can compare the continuous phase shift α₂(ν,t) using the surface boundary conditions of zero oscillating gravitational potential and its derivative, with the value from the full eigensolution; the fractional difference at x = 0.95 is less than 10^-8 which exceeds the accuracy of the calculation! Likewise we can determine the departure of the α_ℓ(ν,t) from an ℓ independent function α(ν,t), the fractional departure from the average, δα_ℓ/α, at x = 0.95 and ν = 2133.16 μHz is 2.7 × 10^-4.

2.3. The phase function

The eigenfrequency equation (Eq. (7)) contains the acoustic radius T, which is not an observable, so we write Eq. (7) in terms of an observable mean large separation Δ (e.g. the value 106.5 μHz in Fig. 1) and a phase function defined by

and the phase shifts α_ℓ,δ_ℓ are evaluated at any radius x_f.

The last term on the right hand side Eq. (12b) is just a linear function of frequency ν_nℓ. If we calculate the phase shifts in the outer layers where the departure of the surface phase shifts α_ℓ from an ℓ independent function of frequency can be neglected (e.g. at x_f = 0.95), then the whole of the right hand side is an ℓ independent function of frequency, α^∗(ν), and hence so too is , which only depends on the structure interior to x_f. In Fig. 4 we show for Model A; it is an ℓ independent function of frequency. This result provides the basis of model fitting by phase matching. If the model star has the same interior structure as an observed star then the function of the internal phase shifts of the model calculated using the observed frequencies should collapse to an ℓ independent function of frequency. The extent to which this is satisfied is a measure of the goodness of fit of a model to the interior structure of the observed star.

The collapse of the phase function to a function only of frequency is independent of the addition of a constant to the n values, and of a constant multiple of π to the inner phase shifts, since both just add a constant to the value of α^∗(ν). But it does depend on knowing the correct ordering of the n values and is therefore only applicable to main sequence stars where the n values can be assigned as in Eq. (11) above.

However the collapse of a phase function to a function only of frequency is not dependent on knowing the n values nor on the specific algorithm used to determine the δ_ℓ. For example were we to defining a phase match function periodic in π e.g. (13)where k is just a label not the n values, this is independent of n and of any multiples of π in the definition of the δ_ℓ, and is determined solely by the structure interior to x_f. However we shall not use this definition due to the difficulty in assigning errors on due to errors on the frequencies.

2.4. Mixed modes

The above example (Model A) is a main sequence star where the frequencies all follow a regular pattern and can therefore be ordered in n values, but for later stages of evolution, where there are mixed modes, the situation is more complicated and the regularity in n values only applies to the ℓ = 0 modes. We give an example of such a star in our evolutionary sequence, Model C, which is well beyond the main sequence and into the hydrogen shell burning phase as the star moves over to the red-giant branch.

The frequency set has 14 ℓ = 0 modes centred on ν_max = 942 μHz, but for ℓ ≠ 0 has many mixed modes all in the range 565 − 1297 μHz. These are displayed in the echelle diagram in Fig. 2 and listed in Table 2, the different frequencies being identified by the label k. The average large separation Δ as defined in Eq. (10) is 55.89 μHz.

In Table 2 we also list the number of nodes in the ψ eigenfunction, which of course is not an observable, and the n values determined on using the observational condition given in Eq. (11). For the ℓ = 0 modes these n values are the same as the number of nodes but for other ℓ values the observationally assigned n values bear little relation to the number of nodes and different frequencies can have the same n value. Moreover there is no reason to assume that the algorithm for assigning n values to main sequence stars is valid for mixed modes, and the n values depend on the specific limits imposed on ϵ in Eq. (11), with different limits, e.g. if we impose 0 <ϵ ≤ 1, some of the n values would be changed and some remain unchanged.

We follow the analysis described above and calculate the inner and outer phase shifts, δ_kℓ,α_kℓ, of the eigenfrequencies ν_kℓ using the definitions in Eqs. (8a) and (8b). These phase shifts are shown in Fig. 6; as the inner phase shifts δ_kℓ range between 0.7 and 36 we have subtracted an integral number of π to bring all the inner phase shifts within the interval { 0,π} on this diagram. The outer phase shifts α_kℓ again all lie on a single curve α(ν).

For frequency sets with mixed modes the algorithm for determining the phase function given in Eq. (12) needs to be modified to be independent of any assignment of n values and any multiples of π in the determination of the inner phase shifts δ_ℓ. An algorithm that does not rely on assigning n values, nor for that matter on the k values assigned to frequencies ν_kℓ with ℓ ≥ 1, is given below.

Since there are no mixed modes for ℓ = 0, the frequencies ν_k0 follow a regular pattern and we can label them in order of increasing frequency (the absolute value does matter), and evaluate the phase function as a function of frequency using only the ℓ = 0 frequencies, (14)For frequencies ν_kℓ, ℓ ≠ 0, we then interpolate in to determine the value g_kℓ at ν_kℓ and form the function (15)If the observed frequencies fit the model then is an integer; the departure from an integer being a measure of the goodness of fit of the model to the observed frequencies. We then determine the integer closest to the value of and define the phase function as (16)which is independent of an addition of an arbitrary multiple of π to the δ_ℓ and reduces to an ℓ independent function ν if . The extent to which this is satisfied is a measure of the goodness of fit of a model to the interior structure of the observed star.

Applying this algorithm ito Model C we derive the phase function which is displayed in Fig. 7, it collapses to an ℓ independent function of ν.

3. Model fitting by phase matching

3.1. The algorithm

Given a set of observed frequencies with error estimates , and a model to be tested for goodness of fit, solve the oscillation equations for the model for each observed frequency subject to the central boundary conditions of regularity and the outer boundary condition on the gravitational potential (but not the pressure). Form the scaled Eulerian pressure perturbation ψ_ℓ(t) = rp′/ (ρc)^1/2 and hence the inner phase shift δ_kℓ at some fitting point in the outer layers (e.g. x_f = 0.95)² given by the solution of (6a) e.g. by using Eq. (8b) above. The phase function is then determined using the algorithm give in Eqs. (12) or (14)−(16), where Δ is an estimate of the average large separation of the observed frequencies.

We determine error estimates e_kℓ on by calculating and subtracting from , and then determine , a best fit function of frequency to the with errors e_kℓ (M being the number of parameters in ), and determine the goodness of fit in terms of (17)where N is the number of frequencies. In practice is modelled as a series in Chebyshev polynomials and M ≤ N₀ is taken to minimise χ², where N₀ is the number of ℓ = 0 frequencies. All models with χ² ≤ 1 have an interior structure that is consistent with the observed frequencies.

We emphasise that the outer layers of the model above x_f do not enter into this algorithm, all that is assumed is that the unknown surface phase shifts α_ℓ(ν) can be modelled by a single function of frequency. We also note that the error estimates on each are uncorrelated (assuming that the error estimates on the frequencies are uncorrelated!), in contrast to the situation in model fitting by comparing separation ratios.

The major problem is the assignment of angular degree ℓ. For main sequence stars with a regular pattern of frequencies – as is the case for Model A (Fig. 1) – this poses little difficulty, but for stars with mixed modes it is not so straight forward. Indeed applying this phase matching algorithm to stars with mixed modes is fraught with problems; a small change in the model produces large changes in the inner phase shifts and renders finding a best fit model exceedingly difficult.

Surface layer independent model fitting algorithms necessarily cannot give information on the outer layers of the star. One may impose other fitting restrictions on the models: mass, luminosity, radius, surface gravity, surface composition, large separation... but these are not tested by surface layer independent model fitting,

3.2. Examples

To demonstrate that the collapse of the phase function to a function only of ν depends only on the structure of the inner layers we take two models with the same interior structure but different outer layers and subject them to the model fitting test. One model is model A as used above, the second is Model G which has the same structure as Model A but we have modified the adiabatic exponent Γ₁, setting values above x_f = 0.95 to the value at x_f = 0.95. The frequencies of Model G are listed in Table 3; as can be seen by comparing these frequencies with those of Model A listed in Table 1, the difference between the two frequency sets is up 57 μHz and the large separation of Model G is 108.15 μHz compared to 106.54 μHz for model A.

We then calculated the phase functions for model A using Model G frequencies, and conversely for model G using Model A frequencies, assuming an error estimate of 0.2 μHz on all frequencies. The results are shown in Fig. 8: both phase functions collapses to a function only of frequency with very small χ², demonstrating that the phase matching procedure is independent of the outer layers of a star. is identical to that shown in Fig. 4 (which is Model A frequencies through Model A) as it should be since the interior structure of Model G is is identical to that of Model A; has a different functional form as the frequencies and large separation of Model G are different from those of Model A. But the test is the ℓ independence of the phase function, not its particular functional form

To demonstrate model fitting by phase matching we take the frequencies of Model A to represent a set of observed frequencies , with error estimates Hz on all frequencies, and test for a fit to other models in the same evolutionary sequence. The first example is Model B, the 1.15 M_⊙ model evolved to a central hydrogen abundance of X_c = 0.009. The results are shown in Fig. 9; does not collapse to a function only of ν. The continuous curve is the least squares fit to a single function of ν, the error estimates on are too small to be seen in this figure. The χ² of the fit is 136; the observed frequencies of model A do not fit model B.

We then repeated this analysis for a set of models in the evolutionary sequence of the 1.15 M_⊙ star from the initial main sequence X_c = 0.718, to terminal main sequence with X_c = 10^-5 at an age of 5.02 × 10⁹ ys. Figure 10 shows the resulting χ² for this series of models. Obviously models that are far removed from Model A have a very poor fit, so to examine more closely the models in the neighbourhood of Model A we used a finer output mesh in the evolutionary calculation and the models with χ²< 1 are shown in Fig. 11. With an estimated error on the frequencies of 0.2 μHz all these models fit the observed input frequencies.

To further restrict the models we need to impose additional constraints. Since the phase match algorithm only tests the inner structure one first looks for additional constraints that are independent of the structure of the outer layers. The two obvious such constraints are mass and luminosity. For an increasing number of stars bolometric flux measurements and parallaxes give the luminosity to high precision, often better that 1% (cf. Boyajian et al. 2013). Imposing the luminosity constraint much reduces the set of models consistent with Model A; this is shown in Fig. 11 where the red squares match the luminosity of Model A to within 0.5% whilst the blue triangles match it to within 1%.

One can of course impose additional constraints such as radius and large separation but these are not surface layer independent. We discuss this below where we seek models that fit the observations of HD 49933 and HD 177153.

4. Application to HD 49933 and HD 177153

4.1. Constraints on model fitting

Surface layer independent model fitting by phase matching only seeks to find models whose internal structure out to some fitting radius r_f = x_fR (e.g. where x_f = 0.95) is compatible with an observed frequency set. The oscillation properties of the model out to r_f are determined by the run of density ρ(r), pressure P(r), and adiabatic exponent Γ₁(r) for r ≤ r_f, and for main sequence stars, dominantly just by the sound speed c(r) (c² = Γ₁P/ρ). We cannot constrain the internal profile of chemical composition or temperature,

As mentioned above in addition to the phase matching condition there are two model constraints that only depend on the interior structure: the mass M and luminosity L. The luminosity can be estimated using the parallax and either the magnitude or measurements of the bolometric flux. The mass may be dynamically determined for stars in a binary system (e.g. α Cen A&B) or, in principle, from from the radius R combined with a spectroscopically determined surface gravity log  g. However the determination of log  g is clearly dependent on modelling stellar atmospheres and large differences can exist between values based on LTE and non-LTE models (Serenelli et al. 2013), and even between different determinations using similar techniques (cf. HD 49933 in Table 5 of Bruntt 2009), so we shall not impose values of spectroscopically determined mass or log  g on our model fitting; we will however will compare our best fit models with the spectroscopically determined log  g. Nor will we use the spectroscopically determined atmospheric abundance as a constraint since the relation between atmospheric and interior composition is fraught with problems of diffusion, gravitational settling and mixing due to circulation and mild turbulence as well as with atmospheric modelling problems.

The radius R and average large separation Δ are dependent on the structure of the outer layers, the outer 5% of radius contributing of the order of 30% to the large separation (cf. Roxburgh 2014). An error of 10% in the structure of the outer layers can give an error of ~ 3% in Δ and ~ 0.5% in radius. We will use estimated values of R as a constraint both with and without this enhancement.

However the average large separation Δ poses many problems – the value being determined by the specific algorithm taken as the definition. For both observed and model frequency sets the large separation Δ_nℓ = ν_nℓ − ν_{n − 1,ℓ} varies with frequency – for HD 49933 the large separation Δ_n0 varies between 82.7 and 89.3 μHz and for HD 177153 it varies between 102.2 and 105.5 μHz. We here define an average Δ for both observations and models as the value at ν_max, the peak of p-mode power in the observed star, obtained by fitting a straight line to N_Δ values of Δ_n0, where N_Δ is the number of values for the observed star. We do not impose a constraint on Δ in our initial model fitting but as a filter on models that fit all other constraints – both with and without the additional 3 μHz uncertainty from possible errors in the contribution of the outer layers.

A further much used constraint on model fitting is the mass estimated using the so called scaling relation Δ ∝ M^1/2/R^3/2, but this is not an exact relation and is subject to considerable errors which, as estimated from model comparison, vary between 5−30% for stars in the mass range 1−1.5 M_⊙ depending on the heavy element abundance Z and other factors that enter the model calculations (cf. Roxburgh 2014). It would anyway be inconsistent to use this method when the large separation is so strongly influenced by the structure of outer layers of the star, so we consider the mass as an unknown to be determined by the model fitting.

4.2. HD 49933

The star HD 49933 was observed by CoRoT in two separate runs; the initial short run was described in Appourchaux et al. (2008) and the combined date from the two runs in Benomar et al. (2009) who tabulated a total of 51 frequencies in the range 1200−2600 μHz for modes of degree ℓ = 0,1,2. These are the frequencies used in our analysis together with an estimated ν_max ~ 1800 μHz. The average large separation at ν_max as given by fitting a straight line through the 16 frequency dependent values Δ_n0 is formally 86.0 ± 0.5 μHz, but the error estimate is unreasonably small – if we fit to even and odd subsets of Δ_n0 (when the errors are uncorrelated), we obtain Δ = 85.0 μHz (odd) and 86.9 μHz (even), so a more realistic estimate is Δ = 86.0 ± 1.0 μHz. However as mentioned above there is added uncertainty of ~ 3 μHz due to possible errors in the contribution of the outer layers. We do not impose a constraint on Δ in our initial model fitting but nevertheless find many models that fit all other constraints and have Δ = 86 ± 1 μHz.

Constraints on the luminosity and radius were derived from angular diameter measurements of θ_LD = 0.445 ± 0.012 mas by Bigot et al. (2011), bolometric flux F_bol = 12.780 ± 0.080 × 10^-8 erg/s/cm² by Boyajian et al. (2013) and the Hipparcos parallax π = 33.69 ± 0.42 mas (van Leeuwen 2007), where we have taken solar values as R_⊙ = 6.960 × 10¹⁰ cm, L_⊙ = 3.846 × 10³³ erg/s/cm². These are given in Table 4. Note that the value for L/L_⊙ in Table 5 is slightly different from that in Boyajian et al. (2013), the difference presumably being due to the values taken for L_⊙ which is not specified in their paper. As described above we estimate an additional uncertainty in R of ~0.005R due to possible errors in the contribution of the outer layers.

Spectroscopic values of log  g have been estimated by several authors and several values are given in Bruntt (2009); they vary between 4.0 ± 0.15 (Kallinger et al. 2010) and 4.28 ± 0.06 (Brunt 2009). As stated above we do not use log  g as a constraint, but calculate its value for our best fit models.

We first searched for best fit models using only the phase match criterion plus the luminosity constraint within a model set, GN, constructed using the STAROX code (Roxburgh 2008) with GN93 relative abundances (Grevesse & Noels 1993), OPAL opacities (Iglesias & Rogers 1996) supplemented by Whichita opacities at low temperatures (Ferguson et al. 2005), EOS5 equation of state (Rogers & Nayfonov 2002) and NACRE reaction rates (Angulo et al. 1999), the mixing length model of convection with mixing length =αH_p, two models of overshooting and mixing from a convective core; one with just chemical mixing a distance of α_oH_p with α_o = 0,0.1,0.2, the second with full entropy and chemical mixing a distance α_od_I where d_I is maximum extension of a convective core as given by the Integral Constraint (Roxburgh 1978; 1989); diffusion was not incorporated. The model set has 31 values of M/M_⊙ = 0.9 − 1.50, 4 values of the initial hydrogen abundance X_H = 0.68 − 0.74, 6 values of heavy element abundance Z = 0.009 − 0.24, 3 values of the mixing length α = 1.6 − 2.4, 4 values of the overshoot parameter α_o = 0,0.1,0.2,0.5 and in 14 different stages of main sequence evolution with the central hydrogen abundance X_c = 0.6 − 0.001.

Some 633 models satisfy both the luminosity constraint and phase match criterion with reduced χ² ≤ 1, with masses in the range M = 1.08 − 1.44 M_⊙. Imposing the radius constraint teduces this to 261, 236 of which satisfy the radius constraint without the enhancement of 0.5%R; the models have masses in the range 1.08 − 1.44 M_⊙ and all have convective core overshooting. Some 38 of these have an average large separation Δ in the range 86 ± 1 μHz, with masses in the range 1.10 − 1.30 M_⊙ a selection of which are given in Table 65 and in Fig. 12 we show the phase match for the 1.20 M_⊙ model. We emphasise that all these models with χ² ≤ 1 fit the data, there is no justification for taking the models with the lowest χ² as the best fit model.

The wide range of models that fit the data is compatible with the wide range of model fits reported by Piau et al. (2009), Benomar et al. (2010), Kallinger et al. (2010), Bigot et al. (2011), Creevey & Bazot (2011), and Liu et al. (2014) whose best fit models vary in mass between 1.12 − 1.35 M_⊙, ages between 2.15 − 3.25 × 10⁹ yr and Z values between 0.005 − 0.20.

In their analyses Kallinger et al. (2010), and Liu et al. (2014) used only modes of degree ℓ = 0,1, since the error estimates on ℓ = 2 modes are very large (> 2 μHz). We therefore repeated our model search using only ℓ = 0,1 modes. We found 18 models which satisfied all the constraints (without the enhancement of 0.5%R) whose details are listed in Table 6. These models are a subset of the 38 models that similarly satisfied all constraints but with ℓ = 0,1,2 modes.There is still a wide spread of models that fit the data.

All the models listed in Tables 5 and 6 satisfy the constraints to within 1σ of the observed values, and the phase match criterion to χ² ≤ 1. But this is too rigid a requirement; there will be many more models with slightly higher χ² whose probability of fitting the data is large enough to consider them as possible models. The obvious conclusion from these results is that the precision on the frequencies is not good enough to constrain the models.

4.3. HD 177153 (aka Perky)

The star HD 177153 was observed by Kepler (KIC 6106415, aka Perky) and analysed by Silva Aguirre et al. (2013) who list a total of 33 frequencies in their Tables 1, 11 for each ℓ = 0,1,2, and give ν_max ~ 2210 μHz and an average large separation, determined from the autocorrelation of the time series, as Δ = 104.0 ± 0.5 μHz. The errors on these frequencies are an order of magnitude smaller than those for HD 49933 so we expect to find much better results than for HD 49933. Silva Aguirre et al. searched for best fit models by two procedures: comparing frequencies with a “surface offset” (cf. Kjeldsen et al. 2008), and comparing separation ratios (cf. Roxburgh & Vorontsov 2003, 2013). Their set of best fit models (Tables 4 in Silva Aguirre et al. 2013). has masses in the range 1.05 − 1.17 M_⊙, ages in the range 3.8 − 5.5 × 10⁹ yr, and Z in the range 0.014 − 0.021.

The average large separation at ν_max, determined by a straight line fits to the 10 frequency dependent values Δ_n0, is formally Δ ~ 104.0 ± 0.1 μHz but again this is unrealistic since neighbouring values of Δ_n0 differ by between 0.14 and 1.60 μHz; we therefore take the error estimate of 0.5 μHz given by Silva Aguirre et al. (2013), to which we may add ± 3 μHz due to possible errors in the contribution of the outer layers. We do not use the observed value of Δ as a constraint on model fitting, but do compare values for best fit model values.

Table 7 lists constraints on luminosity and radius derived from angular diameter measurements of θ_LD = 0.289 ± 0.006 mas and bolometric flux F_bol = 3.39 ± 0.07 × 10^-8 erg/s/cm² given by Huber et al. (2012), and the Hipparcos parallax π = 24.11 ± 0.44 mas (van Leeuwen 2007). We have, as before, added 0.5%R to the error estimates on R.

A spectroscopic value of log  g = 4.40 ± 0.08 is given in Bruntt et al. (2012), which differs from the value of log  g = 4.27 ± 0.1 given in Silva Aguirre et al. As explained above we do not use log  g as a constraint, hence we have no constraint on the mass, but we list the values of log  g for our best fit models.

These constraints are not the same as those in the study by Silva Aguirre et al. (2013) who did not use the angular diameter and flux measurements of Huber et al. (2012). Their luminosity constraint derived from parallax and magnitude 1.66 ≤ L/L_⊙ ≤ 2.00 is less severe than that from the bolometric flux; and they used both the large separation estimated as Δ = 104 ± 0.5 μHz and log  g = 4.27 ± 0.1 as constraints on their model search.

We then searched two model sets to find models that fitted the above constraints, using the full set of 33 frequencies, one set, GN, had the same properties as the set used above but with a finer mesh restricted to 0.9 ≤ M/M_⊙ ≤ 1.20, the second set, GS, was the same as GN except that we used GS98 relative abundances (Grevesse & Sauval 1998) and corresponding opacities.

66 models satisfied the luminosity constraint and the phase match criterion χ² ≤ 1 with masses in the range 1.11−1.19 M_⊙, of which 10 satisfied the radius constraint without the 0.5% addition, and a further 5 with the addition; the properties of these models are listed in Table 8, none have a convective core nor core overshooting. We emphasise that we did not impose a mass, log  g or large separation constraint when fitting these models.

Nevertheless the first 8 models in Table 8, which satisfy the luminosity constraint and (unenhanced) radius constraint to within 1σ and have χ² ≤ 1 all have a large separation within 104.0 ± 0.5 μHz. The following 1.18 M_⊙ model is essentially just as good, but the 1.12 M_⊙ model has a Δ3 μHz off the observed value and perhaps should be excluded. The last 5 models, which all satisfy the enhanced radius constraint have Δ within 1.1 μHz of the observed value. All models except the anomalous 1.12 M_⊙ model have log  g = 4.30 − 4.31. One example of the phase matching is given in Fig. 13.

The 8 best fit models have masses in the range 1.15 − 1.17 M_⊙, ages in the range 4.45 − 4.70 × 10⁹ yr, Z in the range 0.021 − 0.024, XH = 0.71 − 0.72, helium abundance Y = 0.256 − 0.266 and (surprisingly) mixing length parameter α = 1.8. However if one were to relax the constraints a little one would find additional models whose probability of fitting the data is large enough to consider them as possible models.

Our best fit models are not inconsistent with the results of Silva Aguirre et al. (2013) but favour masses at the higher level of their range of values and in general have higher values of Z. But, as pointed out above, the constraints they imposed on their searches for best fit models are not the same as ours, in particular they did not use the bolometric flux and radius from Huber et al. (2012), and used constraints on log  g = 4.27 ± 0.1, and their 8 different models were computed with different assumptions on the physical processes governing stellar evolution. Also they only used a subset of the 33 frequencies in their model fitting.

We therefore ran our model fitting procedure using the constraints from Silva Aguirre et al. (2012; hereafter referred to as SA) which we list in Table 9. As a pointed out above the model value of Δ depends on the definition used so we again used the fit to the 10 values of Δ_n0 surrounding ν_max.

There were 178 models that satisfied the SA luminosity constraint and the phase match criterion χ² ≤ 1, all of which also satisfied the SA log  g constraint; the models had masses in the range 1.07 ≤ M/M_⊙ ≤ 1.19, with average large separations Δ in the range 98 ≤ Δ ≤ 120 μHz. This is in sharp contrast with the 15 models listed in Table 8 which satisfied the Huber constraints on L and R. 24 of the 178 models satisfied the constraint on Δ = 104 ± 0.5 μHz and a further 21 had Δ within the range 104 ± 1.0 μHz.

Since the SA luminosity constraint of 1.66 ≤ L/L_⊙ ≤ 2.00 encompasses the more restrictive constraint from Huber et al. (2012) of 1.738 ≤ L/L_⊙ ≤ 1.890, and as the log  g constraint imposes no restrictions, all the models that satisfy the Huber constraints necessarily satisfy the SA constraints; the major difference being that smaller mass models satisfy the SA constraints. Table 10 lists all the 24 models that satisfy the SA criteria including Δ in the range 104 ± 0.5 μHz, and below the line a few of the additional models that have Δ in the range 104 ± 1.0 μHz. The masses now range from 1.07 − 1.19 M_⊙, ages from 4.42 − 5.11 × 10⁹ yr, and Z from 0.018 − 0.024.

The mass range and ages are similar to those in Silva Aguirre which are 1.05 − 1.17 M_⊙ and ages from 3.38 − 5.53 × 10⁹ yr, whilst their Z values range from 0.014 − 0.021. But we emphasise that our best fit models are drawn from our set of models all with the same input data (except for the 2 different relative abundances, GN93 and GS98) and one particular set of modelling assumptions which are not the same as the 8 different sets of modelling assumptions for the 8 different best fit models in Silva Aguirre et al. (2013).

The tentative conclusion we can draw from this is that the more restrictive constraints on L,R from angular diameter and bolometric flux measurements by Huber et al. (2012) favour larger masses and give best fit models with no convective core nor overshooting, with masses in the range 1.15 − 1.17 M_⊙, ages in the range 4.45 − 4.70 × 10⁹ yr, Z in the range 0.021 − 0.024, XH = 0.71 − 0.72, helium abundance Y = 0.256 − 0.266 and mixing length parameter α = 1.8.

But we have only explored model sets without diffusion and with one specific set of model physics; there may well be other models constructed with different model physics that also fit the data, but we need higher precision data to differentiate between them.

5. Accuracy of the approximations

Surface – layer – independent model – fitting techniques, phase matching or separation ratios, are based on the approximation that the frequency dependent outer phase shifts α_ℓ(ν) are, very nearly, independent of ℓ in the outer layers, and can all be replaced by an ℓ independent function of frequency α(ν).

To test the validity of this approximation we calculate δα_ℓ, the difference between the α_ℓ and their average over values ℓ = 0,1,2,3, as a function of frequency and fractional radius x = r/R for four models; the three 1.15 M_⊙ Models A, B, C as defined above, and also for a solar model (Model S of Christensen-Dalsgaard et al. 1996). All frequency sets, except Model C, have 14 frequencies of each degree (ℓ = 0,1,2,3) centred on a frequency normalised to a solar value of 3050 μHz. Model C, with many mixed modes, has 14 ℓ = 0, 18 ℓ = 1, 18 ℓ = 2, and 20 ℓ = 3 frequencies (see Fig. 5). For all frequency sets the maximum value of δα_ℓ(ν) almost always occurs for the lowest frequency ℓ = 3 mode, as expected on theoretical grounds since the dominant ℓ dependent term in the oscillation equations ∝ ℓ(ℓ + 1) c²/ (ν²r²).

Figure 14 shows the results – giving the maximum departure δα/α as a function of fractional radius on using the continuous phase shifts, the error decreases with increasing radius but is quite large at smaller radii. Clearly we should choose a fitting point x_f for the phase matching well into the outer layers.

A further requirement is that the error in fitting the phase function for the discrete values ν_nℓ with a series of Chebyshev polynomials (or any other set of basis functions), is sufficiently small not to be important in the phase matching algorithm. One way of testing this is to compare the δα_ℓ from the continuous α_ℓ(ν) with those from using the discrete values. This is shown in Fig. 15 for Model A at a fitting radius x_f = 0.95, where for the discrete values the average has been computed by fitting the full set of 56 α_ℓ(ν_nℓ) with a series of Chebyshev polynomials up to degree 13. The errors are essentially the same as those from the continuous α_ℓ(ν) with | δα/α | < 4 × 10^-4.

We then repeated the calculation for Models B, C and S and for mode sets of degree ℓ = { 0,1 } , { 0,1,2 } , { 0,1,2,3 }. The fractional errors are given in Table 11, all the errors are small but increase with increasing ℓ.

We checked the accuracy with which the continuous partial wave outer phase shifts α_ℓ(ν), computed using the boundary condition that the oscillating potential and its derivative are zero at the surface, reproduce the exact values for the eigenfrequencies ν_nℓ; the fractional error is absolutely negligible: zero for ℓ = 0,1 and < 10^-6 for ℓ = 2,3 for all 4 models. We also checked the error in using the Cowling approximation (neglect of the oscillating gravitational potential); here the error is of course zero for ℓ = 0, but is still less than 10^-6 for ℓ = 1,2,3 for all 4 models.

We then checked that the continuous partial wave inner phase shifts δ_ℓ(ν) (ℓ = 0,1,2,3) computed applying just the boundary condition on the gravitational potential at the surface, or at a fitting point x_f = 0.95, reproduce the exact values of δ_ℓ for the eigenfrequencies. The errors Δδ are small for all models, the largest errors are when applying the boundary condition at x_f and for the main sequence Models A, B, S are less 10^-6 and for Model C, which has mixed modes, are less than 10^-5, smaller with the errors in α(ν).

The second criterion for the choice of the fitting point x_f is that the δs should not be influenced by the structure of the outer layers and in particular the HeII ionisation layer, since they then have a short period quasi periodic modulation of large amplitude, the choice of x_f = 0.95 satisfies this criterion, as do other values in this region.

The error in taking α_ℓ(ν) = α(ν) limits the accuracy with which the phase function can collapse to a function only of frequency. For model A with 0.5 <δ_ℓ< 1.5 and 2 <α< 4 (see Fig. 2), the ℓ dependent contribution to from the δs is of order 1 /π ~ 0.3 (with errors < 10^-5), whereas the ℓ dependence of α/π is of order 3 × 10^-4, much smaller than the contribution from the δs.

For real data, with errors on the frequencies σ ~ 0.2 μHz, and a large separation Δ ~ 100 μHz, the dominant error in comes from σ/ Δ ~ 2 × 10^-3, larger than the errors from δα_ℓ. We conclude that the approximation α_ℓ(ν) = α(ν) in the outer layers is valid.

The one exception is the Sun; with the very high precision achieved by space and ground based experiments it is necessary to take the outer phase shifts to be represented by two functions α_ℓ(ν) = α(ν) + ℓ(ℓ + 1) α₂(ν) (cf. Roxburgh & Vorontsov 1994).

6. Relation to separation ratios

Surface layer independent model fitting, phase matching and separation ratios, are both based on the result that in the outer layers of a star the surface phase shift α_ℓ(ν) is, to a very good approximation, independent of ℓ so that α_ℓ(ν) = α(ν) which can be considered a continuous function of ν.

For p-modes in main sequence stars where the radial order n can be assigned (such as Model A) and the average large separation Δ estimated, the oscillation frequencies ν_nℓ can always be expressed in terms of functions ϵ_ℓ(ν) in the form (18)and their values are known at the discrete set of frequencies ν_nℓ.

From the eigenfrequency equation (Eqs. (7), (12)), ϵ_ℓ(ν) can be expressed in terms of the inner and outer phase shifts as (19)With α_ℓ(ν) = α(ν) independent of ℓ, the term in square brackets in this equation (which we call α^∗) is a function only of ν.

The phase match algorithm uses this equation in the form (20)is a function only of ν; the separation ratios algorithm uses it to subtract off the α^∗(ν) by interpolating in the discrete set of the ϵ_ℓ(ν_nℓ) for values at the same frequencies, and eliminating α^∗ by subtraction to determine the inner phase shift differences, e.g. (21)which only depend on the inner structure of the star. If one uses linear interpolation one readily obtain the result that where r₀₂,r₀₁ (and its relative r₁₀) are the separation ratios.

A source of error in the separation ratio technique is the error in interpolation; this can be reduced by using higher order interpolation algorithms, such as the 5 point expressions for r₀₁,r₁₀ (Roxburgh & Vorontsov 2003a), but they still have error estimates which can be significant if the quasi periodic modulation of the δ_ℓ is on a scale smaller than twice the large separation.

In Fig. 16 we illustrate this by interpolating for the values of δ₀(ν) at frequencies ν_n1 from the values at ν_n0 for Model A (all values are multiplied by ν^∗ = ν/ (2200 μHz) for clarity of presentation and evaluated at x_f = 0.95). The curve is the continuous phase shift δ₀(ν), the black points are the discrete phase shift δ_n0 at the frequencies ν_n0, the red triangles the values δ₀(ν) at frequencies ν_n1, and the blue stars the values at ν_n1 determined from the δ_n0 by cubic interpolation.

The quasi periodicity in the δs is due to the rapid variation in acoustic variable at the base of a convective envelope, the period being determined by the acoustic radius of the interface and the amplitude by the strength of the quasi discontinuity (Roxburgh 2009). It is only an important source of error if the error from interpolation is larger than that due to errors on the frequencies. This is not the case for Model A, Fig. 17 shows the values of the 5 point small separations r₁₀ compared with to actual phase shift difference, the error bars corresponding to an error of 0.1 μHz on the frequencies.

A difficulty with the separation ratios technique is that since it uses several frequencies to determine each value of the ratios, the errors estimates on these ratios, derived from the errors on the frequencies, are strongly correlated. This is not the case with the phase matching technique since depends solely on the frequency ν.

We note here that one can interpolate for the value of ϵ₀ and ϵ_ℓ at any ν and, on subtraction, obtain an approximation to the phase shift difference δ_ℓ − δ₀ at any ν, but this does not increase the number of independent variables. If one has N values of each of ν_n0 and ν_nℓ, then essentially N of them are used to model α^∗(ν) leaving N independent data on δ₀ − δ_ℓ.

7. Conclusions

We have defined a phase function in terms of inner phase shifts δ_ℓ(ν), determined solely by the interior structure of a model star, and shown that it collapses to an ℓ independent function of frequency when the frequencies are the eigenfrequencies of the model. This then gives the surface layer independent phase match algorithm (defined in Sect. 3) for seeking models whose interior structure is compatible with the observed frequencies of a star by the goodness of fit of the phase function , calculated in the model star at the observed frequencies , to an ℓ independent function of frequency. The algorithm does not require one to know the radial order n of the observed frequencies, the subscript k here just being a label to identify individual frequencies and, in principle, it can be applied to stars with mixed modes; it does however require identification of the angular degree ℓ, in common with other model fitting procedures.

The phase shifts δ_ℓ(ν) that enter into the calculation of the phase function are defined in Sect. 2 and an algorithm given for calculating their values for a model. In Sect. 3 we define the phase match algorithm and apply it to models to demonstrate its diagnostic properties.

In Sect. 4 the algorithm has been applied to two stars: HD 49933 observed by CoRoT, and HD 177153 (aka Perky) observed by Kepler, where we search model sets constructed with particular input physical properties (GN93 and GS98 relative abundances, and no diffusion). We find a very wide range of models that satisfy the phase match criterion alone (i.e. χ²< 1). The set of models is reduced if we require the models to be consistent with observational constraint on the luminosity, which is determined solely by the inner structure of the star. The mass is also determined solely by the inner structure of the star but we have no constraint on the mass. The set is further reduced if we require compatibility with the observed radius. For HD 49933 there are many models that satisfy all these constraints, including matching the large separation, with masses between 1.10 − 1.30 M_⊙; the errors on the frequencies are too large to discriminate between these models.

For HD 177513 (aka Perky) the uncertainty on the frequencies is considerably smaller than for HD 49933 and we find much tighter constraints on models that satisfy the luminosity and radius constraints as well as the phase match criterion; again we have no constraint on the mass. The 8 best fit models also satisfy the large separation constraint and have masses in the range 1.15 − 1.17 M_⊙, ages of 4.45 − 4.70 × 10⁹ yr, Z in the range 0.021 − 0.024, XH = 0.71 − 0.72, Y = 0.256 − 0.266 and mixing length parameter α = 1.8; none have convective cores. If we relax the radius constraint a little (to allow for the possible uncertainties in the contribution of the outer layers) a slightly larger set of models with masses in the range 1.14 − 1.18 M_⊙ having similar properties are compatible with the data.

The phase match algorithm rests on the fact that the contribution of the outer layers of a star to its frequencies is, to a good approximation, independent of angular degree ℓ, a result which enables one to subtract out this contribution. The accuracy of this approximation is discussed in Sect. 5

using stellar models and shown to be more than good enough to validate the approximation.

In Sect. 6 we compare the phase match technique to model fitting using separation ratios, which also rests on the ℓ independence of the contribution from the outer layers. The separation ratio algorithm requires the radial order of modes, uses local interpolation to cancel out the effects of the outer layers, which introduces errors, and the error estimates on the ratios are strongly correlated, none of which are present in the phase match algorithm. This will be important when using high precision data.

---

Acknowledgments

The author thanks Sergei Vorontsov for many valuable contributions to joint work on phase shift analysis. He also gratefully acknowledges support from the Leverhulme Foundation under grant EM-2012-035/4.

References

All Tables

All Figures