Issue 
A&A
Volume 626, June 2019



Article Number  A108  
Number of page(s)  11  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201935669  
Published online  21 June 2019 
Interpolation of equationofstate data
^{1}
Sternberg Astronomical Institute, M.V. Lomonosov Moscow State University, 13, Universitetskij pr., 119234 Moscow, Russia
email: avo@sai.msu.ru
^{2}
Department of Physics and Astronomy, University of Southern California, Los Angeles, CA 90089, USA
Received:
12
April
2019
Accepted:
26
April
2019
Aims. We use Hermite splines to interpolate pressure and its derivatives simultaneously, thereby preserving mathematical relations between the derivatives. The method therefore guarantees that thermodynamic identities are obeyed even between mesh points. In addition, our method enables an estimation of the precision of the interpolation by comparing the Hermitespline results with those of frequent cubic (B) spline interpolation.
Methods. We have interpolated pressure as a function of temperature and density with quintic Hermite 2Dsplines. The Hermite interpolation requires knowledge of pressure and its first and second derivatives at every mesh point. To obtain the partial derivatives at the mesh points, we used tabulated values if given or else thermodynamic equalities, or, if not available, values obtained by differentiating Bsplines.
Results. The results were obtained with the grid of the SAHAS equationofstate (EOS) tables. The maximum lgP difference lies in the range from 10^{−9} to 10^{−4}, and Γ_{1} difference varies from 10^{−9} to 10^{−3}. Specifically, for the points of a solar model, the maximum differences are one order of magnitude smaller than the aforementioned values. The poorest precision is found in the dissociation and ionization regions, occurring at T ∼ 1.5 × 10^{3}−10^{5} K. The best precision is achieved at higher temperatures, T > 10^{5} K. To discuss the significance of the interpolation errors we compare them with the corresponding difference between two different equationofstate formalisms, SAHAS and OPAL 2005. We find that the interpolation errors of the pressure are a few orders of magnitude less than the differences from between the physical formalisms, which is particularly true for the solarmodel points.
Key words: equation of state / methods: numerical / Sun: evolution / Sun: interior / stars: evolution / stars: interiors
© ESO 2019
1. Introduction
The equilibrium equation of state (EOS) of the plasma is a chief physical component in modeling the evolution and structure of stars. At present, several equations of state are widely used: CEFF (Eggleton et al. 1973; ChristensenDalsgaard & Däppen 1992), MHD (Hummer & Mihalas 1988; Mihalas et al. 1988; Däppen et al. 1988), OPAL (Rogers et al. 1996; Rogers & Nayfonov 2002), FreeEOS (Irwin 2012), SAHAS (Gryaznov et al. 2006, 2013) and others. Different EOSs lead to slightly different results of stellar evolution modeling. For example, Morel et al. (1997) and Buldgen et al. (2019) have built solar models with the aforementioned equations of state and found differences, for example, in the soundspeed profile, helioseismic inversions, and the position of the bottom of the convection zone. Turning to stellar modeling, Somers & Pinsonneault (2014) investigate the influence on lithium abundance in stars at an early stage of evolution; the EOSinduced difference reaches 0.81 dex. Joyce & Chaboyer (2018) state that the choice of the equation of state affects the evolutionary tracks of stars with a mass of less than 0.65 M_{⊙}. Brito & Lopes (2018) demonstrate the importance of accurate EOS calculations in ionization regions to model convective envelopes of Fstars for asteroseismic analysis.
At the same time, fine helioseismic inversion indicates that the adiabatic exponent Γ_{1} in the equation of state of solar plasma can be determined from observational data with an accuracy on the order of 10^{−4} (Vorontsov et al. 2013). Therefore, a theoretical EOS for the solar interior that comes as close to reality as possible is a crucial question for astrophysics. Currently, there are many theoretical formalisms and practical equations of state available. To evaluate their absolute accuracy, currently astrophysics alone can deliver the observational data that constrain the theories. However, the use of astrophysical data relies on modeltheory comparisons, which need highquality stellar models. The equationofstate part of these models must be a highprecision realization of the theory, irrespective of the accuracy of the theory itself. Keeping the interpolation errors in equationofstate formalisms at a minimum is therefore an important part of the modeling task.
The plasma inside stars is close to a reacting mixture of ideal gases. Therefore, accurate calculations of thermodynamic functions require laborious computations and extensive computational resources.
An equation of state is a set of thermodynamic functions which consists of thermodynamic potential and its first and second partial derivatives. Namely, these derivatives can be physically measurable, in contrast to the potential itself. In the framework of the chemical picture (Krasnikov 1977), a suitable potential is Helmholtz free energy F, which is convenient because it has relatively simple expressions for mixture of idealgas components. Expression for free energy describes also reactions of conversion for ideal components and can include small corrections for particle interaction. Explicit expressions for pressure and entropy can be obtained as first derivatives of free energy from explicit form of F.
The second derivatives of free energy describe response functions, that is, specific capacity, isothermal susceptibility and thermal expansivity (Reichl 1998). The most thorough explicit analytical derivatives of the free energy were obtained for MHD (Däppen et al. 1988). However, using such explicit expressions does not solve the problem of time consuming computations, because they are based on the degrees of reactions (Reichl 1998), which have to be obtained from large nonlinear system of Saha equations. Thus, the analytical approach is not unique. For example, SAHAS EOS uses analytical first derivatives and numerical differentiation to estimate the second derivatives.
As a result, analytical function of EOS is represented on practice by tabulated functions and its derivatives. Interpolation of tabulated data is a replacement of the complex function by more simple piecewise polynomial function.
Interpolation of tabulated data solves the problem of computational speed because interpolation polynomials are elementary functions and they have all the required derivatives. The polynomial derivatives are continuous and obey to differential relations, similar to Maxwell conditions, just like analytical thermodynamic potential. The tradeoff for simplification is discontinuity at knots and mesh boundaries. Polynomial degree and the maximal order of continuous derivative are determined by type of interpolation spline. In our work, we have considered two types of odddegree splines. The first type is splines with maximal number of continuous derivatives, which are called natural or Bsplines. The second type is Hermite splines, which have fewer continuous derivatives but more free parameters.
While using Bsplines, the common approach consists in independent spline interpolation of thermodynamic functions. The analytical structure of splines differs from the structure of free energy. Hence, even high accuracy interpolation of the potential does not guarantee good accuracy of its derivatives obtained by differentiation of the natural spline. In practice, the accuracy of derivatives is often more important than accuracy of the potential. So, interpolation by independent splines seems relevant. The main disadvantage of the independent interpolation is violation of consistency between the derivatives (the Maxwell relations).
Consistency can be saved through construction of Hermite polynomials. The Hermite spline interpolates not only function but also simultaneously its derivatives at knots. This method preserves information about tabulated thermodynamic functions, and algebraic structure of the spline becomes closer to the original. Equation of state is described by one spline, which interpolates free energy, in contrast to independent interpolation approach. Pressure and response functions are obtained by differentiation of this spline, and thermodynamic consistency is automatically fulfilled.
This approach has been considered in the work by Swesty (1996), in which quintic Hermite splines were used to interpolate a relatively simple equation of state for the purpose of hydrodynamical simulations. This method does not provide the highprecision thermodynamic quantities needed in models of the solarstructure. Therefore we used a completely different approach in our paper. Instead of free energy, we interpolate pressure and its partial derivatives with respect to temperature and density. Thus, input parameters are namely those functions which are of particular interest for practical computations. 2D Hermite splines need more derivatives in mesh knots than in 1D case. Not all the derivatives have clear physical meaning and are presented in tabulated or analytical form. To estimate the missing derivatives in mesh knots, we differentiated Bsplines of the appropriate surface.
Our aim is to describe a type of interpolator that preserves the differential and thermodynamic identities between the functions of EOS. Section 2 briefly describes the main thermodynamic definitions of the equation of state. Accuracy of a traditional interpolation by independent splines is estimated using differentialgeometric identities. Section 3 describes interpolation method with Hermite splines and comparatively discusses its advantages and limitations. Section 4 compares differences between the SAHAS and OPAL 2005 EOS to provide a requested level of accuracy for a spline interpolation. Section 5 provides the main conclusions of the work.
2. Equation of state: Basic concepts and principles of classical interpolation
2.1. Thermodynamic structures and relationships
To describe the equilibrium thermodynamics of plasma within a star, a couple of functions, such as pressure P(T,ρ) and entropy S(T, ρ), are required, each of which can be represented by a twodimensional manifold of temperature T and density ρ. Since under the relevant conditions there are no phase transitions, these functions can be considered as differentiable and their derivatives are exist at least up to second order. However, the functions and their partial derivatives are not completely independent since they obey to thermodynamic identities. A thermodynamic potential allows a compact description of all thermodynamic properties by a single function, for instance, the Helmholtz free energy F(T,V). For definitions and properties of thermodynamics potentials details, see, for example, the book by Reichl (1998). The differentialgeometric structure of the potential is a manifold with its firstorder partial derivatives on the first level and the second and thirdorder partial derivatives at the higher levels. Specifically, the case of the free energy F is illustrated in Fig. 1.
Fig. 1. Differentialgeometric structure based on the Helmholtz free energy. 

Open with DEXTER 
Each higher level of the structure of the potential in Fig. 1 is obtained as a result of differentiation with respect to temperature ∂/∂T (designated by right and downwardpointing arrows in the figure) and to volume ∂/∂V (left and downwardpointing arrows). Obtained after differentiation, thermodynamic functions are shown in the appropriate cells: the pressure P and the entropy S on level 1, the logarithmic pressure derivatives χ_{T}, χ_{ρ}, and the specific heat capacity c_{V} on the second level:
The functions F, P, and S provided on the first level, must each obey the condition of “normality” (Cauchy condition or Maxwell relations), which is the condition of the equality of mixed derivatives. These relations appear when one follows each arrow in Fig. 1. For example, in the case of a transition from the F to level 2, the normality condition leads to the identity
For a transition from the first to the third level, the normality condition leads to the identities
Thus, a single thermodynamic potential F(T, V) gives all the necessary quantities, including the fundamental relations between them.
For the sake of completeness, we have added a definition of the adiabatic exponent Γ_{1} because of its fundamental role in stellar modeling. It can be calculated with the help of the partial derivatives of the potential as follows:
where the dimensionless heat capacity is introduced by
with the scaled pressure defined as
The central idea of our paper is the construction of a geometrical structure, which represents simultaneously the function and its derivatives. The construction of a differentiable manifold leads to an automatic obedience of the thermodynamic relations. We realize the manifold with spline functions of pressure P(T,ρ) alone. Of course pressure and its derivatives only cover a part of all thermodynamic quantities, but they are everything needed in stellar model calculations. In particular, the free energy F itself is not used in such applications.
2.2. Behavior of interpolated thermodynamic functions
To estimate the complexity of the interpolation effort, let us consider the structural behavior of the functions P, Γ_{1}, χ_{T}, χ_{ρ}, c_{V}, computed for a fixed representative chemical composition (hydrogen mass fraction X = 0.8; mass fraction of elements heavier than helium Z = 0.02). Functions are shown as 2Dmanifolds (surfaces) in coordinates T and Q_{s} = ρ · (T_{6})^{−2.25}, where T_{6} = T/10^{6}. These variables are used in the SAHAS tables and described in details by Baturin et al. (2017). We note that the table domain of SAHAS is rectangular.
The surface of the total pressure P (Fig. 2a) shows a wide range of variation: about 20 orders of magnitude, which prevents the discussion of any characteristic details. Therefore, a plot of the scaled pressure Π is much more revealing (Fig. 2b). In the case of the ideal gas law, the scaled pressure Π is inversely proportional to a molecular weight, and it shows up as a horizontal surface. In the regions of ionization and dissociation, the molecular weight rapidly increases with temperature. These areas look like steps on the surface. The dissociation of molecular hydrogen takes place at temperatures of about 1.5 × 10^{3} K, the band where the hydrogen is halfionized passes at (10−20) × 10^{3} K, and helium ionization regions lies at temperatures of (20−100) × 10^{3} K. Ionization of hydrogen is going in asymmetrical way: while hydrogen is rapidly ionized from neutral state, reaching almost complete ionization spreads over wide band in temperature. As result, the hydrogen ionization region and two helium ionization regions overlap, and are hardly distinguishable. The scaled pressure increases at high temperatures and low densities (low Q_{s}) due to the contribution of radiative pressure. Also, the scaled pressure increases over its classical ideal value in the corner of high temperatures and high densities (or high Q_{s}) due to partially degenerate electrons.
Fig. 2. Pressure P (a) and scaled pressure Π (b) in SAHAS EOS for X = 0.80 and Z = 0.02. Red lines show points (T, ρ) from the solar model. 

Open with DEXTER 
The red line shows points of temperature and density taken from a solar model. The model used was calculated with the CESAM2k stellar evolution code (Morel & Lebreton 2008) and represents one of the various standard solar models.
The pressure derivatives exhibit highly significant thermodynamical features. Since pressure itself is only given at discrete points, its values cannot provide the derivatives. Therefore, the EOS tables usually also include the pressure derivatives (typically as logarithmic derivatives). They are computed with finitedifference methods in the mesh points, together with pressure. This is the main justification for the construction of a splinemanifold. The behavior of the logarithmic derivatives is shown in Fig. 3. Differentiation of P emphasizes all the physical features of the scaledpressure surface. Conditions χ_{T} = χ_{ρ} = 1 are equivalent to the idealgas equation P = R_{g}Tρ/μ with constant molecular weight μ. R_{g} is gas constant. If these conditions are fulfilled in some domain of temperatures and densities, then the scaled pressure is represented by horisontal surface in this domain (Fig. 2b). Ionization regions now look like narrow ridges, with maximum temperatures that weakly depend on density (or on Q_{s}). The differentiation of pressure with respect to temperature leads to a much larger amplitude of the ridges. In addition, the contribution of radiative pressure is clearly seen on the derivatives χ_{T} and χ_{ρ} as well as corresponding effect of electron degeneracy.
Fig. 3. Logarithmic derivatives of pressure χ_{T} (a) and χ_{ρ} (b) in SAHAS EOS for X = 0.80 and Z = 0.02. Red lines indicate points (T, ρ) from the solar model. 

Open with DEXTER 
Besides pressure P as a function of T and ρ, the specific heat capacity c_{V} (or equivalent caloric value) is needed to get a complete thermodynamic description. Figure 4 shows the surface of heat capacity c_{V} in units of the scaled pressure Π (see Eq. (6)). The behavior of C_{Π} resembles that of χ_{T}, but there are some differences. First, the amplitude of the ridges varies with density (or with Q_{s}) more rapidly than in the case of χ_{T}. Second, electron degeneracy does not significantly affect the C_{Π} surface.
Fig. 4. Specific heat capacity in units of scaled pressure in SAHAS EOS for X = 0.80 and Z = 0.02 (a) over the whole SAHAS domain of definition, the red line indicates points (T, ρ) from the solar model, (b) for points (T, ρ) from the solar model. 

Open with DEXTER 
The adiabatic exponent Γ_{1} is defined by Eq. (5) and it depends both on the pressure derivatives and heat capacity. We used Γ_{1} to estimate the accuracy of the interpolation. The Γ_{1} surface is shown in Fig. 5a and, separately, its profile in the solar interior in Fig. 5b. All the effects listed above except the electron degeneracy are clearly manifested.
Fig. 5. Adiabatic exponent Γ_{1} in SAHAS EOS for X = 0.80 and Z = 0.02 (a) over the whole SAHAS domain of definition, the red line indicates points (T, ρ) from the solar model, (b) for points (T, ρ) from the solar model. 

Open with DEXTER 
2.3. The second derivative of pressure with respect to temperature
Constructing a Hermite interpolator that takes into account thermodynamic relations requires partial derivatives up to the second order. Some of these derivatives with respect to temperature may be expressed through the derivatives with respect to density (see Eqs. (2) and (4)). The expression (4) is particularly important, because it provides the second derivative with respect to temperature, which shows sharp features and is difficult for numerical interpolation. The classical thermodynamic relation Eq. (4) (see for example Landau & Lifshitz (1980), Sect. 16, Eq. (16.1)) may be rewritten by a transformation from volume V to density ρ:
Thus, the second pressure derivative with respect to temperature can be expressed via χ_{T} and the derivative of c_{V} with respect to density:
For the sake of simplicity in the spline differentiation, we rewrite (9) with dimensionless function C_{Π} instead of c_{V}:
2.4. Accuracy of conventional interpolation
Thermodynamic computations in astrophysics require a globally smooth interpolation of the functions of the equation of state. Conventionally, standard thirddegree polynomial splines (i.e., Bspline) are commonly used for this purpose. The basic definitions of spline theory are given in Appendix A. We have used the notation B[f̂] for this Bspline interpolation of a function f. The cap over the function f̂ denotes that f is given on a set of discrete points.
To construct a Bspline, one needs only the values of the function itself at the nodes of mesh (and some sort of boundary conditions). But once the spline is constructed, the partial derivatives are also available. In the case of cubic Bsplines, we were able to estimate the first and second partial derivatives.
Commonly, the spline derivative ∂B[f̂] differs from the true derivative of the function f′. Moreover, the difference Δ_{∂B} = f′ − ∂B[f̂] is not zero even at the nodes. In practice, we did not have an exact analytical estimation of f′; but instead, the derivatives f̂′ are obtained by finitedifference method at the stage of calculating the EOS tables. Their spline interpolation is denoted by B[f̂′]. The difference Δ = B[f̂′] − ∂B[f̂] can be calculated and used to estimate accuracy of spline differentiation ∂B[f̂].
We illustrate this with the example of SAHAS tables. The tabulated pressure ln P̂ and its derivatives , are interpolated with B[ln P̂], , and (Baturin et al. 2017). As a result, the spline derivative (with respect to temperature, for example) ∂_{T}B[ln P̂] is not equal to . The corresponding differences
and
are shown in Fig. 6.
Fig. 6. Differences between the pressure derivatives calculated as the interpolator of the derivative and the derivative of the interpolator (Eqs. (11) and (12)) on the SAHAS mesh for X = 0.80 and Z = 0.02. Please note the difference of an order of two between the scales in the top and bottom panels. 

Open with DEXTER 
The greatest differences of Δ_{T} are in the ionization and dissociation regions, since the functions ln P̂ and vary there significantly. The maximum error is equal to 0.01, although the typical error in most part of the domain of applicability varies within the range of 10^{−8}−10^{−6}.
The differences of the derivatives with respect to density are remarkably smaller. The maximum deviations occur at the edges of the tables and do not exceed 10^{−4} (see Fig. 6). In the rest of the region, the difference varies from 10^{−10} to 10^{−7}. Thus, differentiation with respect to density is more accurate (that is, it leads to smaller deviations) than with respect to temperature. This fact is used in the production of the necessary derivatives in the Hermitespline construction.
The accuracy of the Bspline interpolator can be estimated also in other way, based on the reciprocity relation used in the thermodynamic identities. For a function of two variables f(x,y) (smooth with second derivatives), the matching condition for mixed derivatives must be fulfilled:
Such mixed derivatives are not presented in the uptodate EOStables, so we need to use splinedifferentiation for this value. There are two ways to obtain the derivative of pressure
The first is by differentiating with respect to temperature, and the second one by differentiating with respect to density. According to Eq. (13), the difference
has to be zero. In practice, it reaches values on the order of 0.04 in the ionization regions and lies in the range of 10^{−7}−10^{−6} in the remaining regions (Fig. 7).
Fig. 7. Difference between the second mixed derivatives of pressure (Eq. (14)) on the SAHAS mesh for X = 0.80 and Z = 0.02. 

Open with DEXTER 
3. Quintic Hermite spline
3.1. Algorithm
The Hermite spline (Hspline) is sometimes defined as a spline that simultaneously interpolates both the function and its derivatives. The interpolation conditions for the function and its derivatives (see Eqs. (A.2) and (A.3) in Appendix A) uniquely define an Hspline. Thus, Hsplines solve the problem of representing the thermodynamic structure of pressure together with its derivatives. Hermite spline can be defined for any odd degree of polynomial. We considered an Hspline of the fifth degree. The choice of the fifth degree instead of third in classical approach allows not only first, but also second derivatives to be continuous at all boundaries of the interpolation cells.
The number of interpolation conditions (A.2) is equal to that of (A.3), which means that Hermite spline can be calculated locally without additional boundary conditions. To interpolate a function at a given point, information is needed only at the boundary nodes of the corresponding interval (or cell in the twodimensional case). To calculate the Hspline inside the interval, one does not necessarily need to know the surrounding elements and the properties of the mesh around it. This property of Hsplines prevents a disturbance due to the irregularity of the mesh or the boundary conditions.
One has to note, however, that Hsplines cannot provide estimations of the derivatives at the nodes and boundaries. Instead, a system of derivatives must be given at the corner knots in order to define Hsplines. In the twodimensional case, it becomes necessary to specify mixed derivatives listed in Appendix B.
The algorithms for constructing Hermite splines are presented, for example, in the work of Hsu (2010), as well as in the lectures by Finn (2004). For a function of two variables, a quintic Hspline is represented in the form
Thirtysix coefficients a_{ij} must be available for its construction. They are determined from the conditions at the nodes of the cell into which the given point falls. Nine values are needed in each of the four nodes:
The independent variables for 2D SAHAS tables are
The interpolated function is the logarithm of pressure
Below we provide all expression for the derivatives of the matrix (16) via the tabulated functions and their splinedifferentiation. Two first derivatives (16) are calculated using and :
The computation of mixed and higher derivatives in (16) is based on the spline differentiation of and . The functions in the lgQ_{s} direction are smoother than along lgT (Fig. 3). So, differentiation with respect to lgQ_{s} is preferable whenever it is possible:
Detailed formulae that explicitly relate to the higher derivatives of pressure with χ_{T} and χ_{ρ} are given in Appendix B.
The three derivatives in the third column of (16): ∂^{2}H/∂u^{2}, ∂^{3}H/∂u^{2}∂v and ∂^{4}H/∂u^{2}∂v^{2} contain a double differentiation of lgP with respect to lgT. One way to calculate them is to differentiate Bsplines of and with respect to lgT, but this leads to an inaccurate result, because and vary greatly with temperature. A better way is to avoid differentiation with respect to lgT altogether by using c_{V} and χ_{T}, which are derived from the thermodynamic relations in Sect. 2.3. Equation (10) allows one to estimate the second derivative
as well as the higher derivatives, which are obtained from it by differentiation with respect to lgQ_{s}:
The formulas for these derivatives in terms of (T,Q_{s}) are given in Appendix B.
3.2. Comparison of B and Hsplines
Figure 8 shows the difference in lgP, calculated with B and Hsplines. Panel a displays the data for the entire domain of applicability of SAHAS. The greatest deviations take place in ionization and dissociation zones, where the discrepancy is on the order of 10^{−4}. In the highertemperature regions, they are five orders of magnitude lower. The red line shows the conditions (T, ρ) of the solar model. The deviations with respect to these points are also shown in panel b. They do not exceed 2.5 × 10^{−6} in the lowtemperature region while in the hightemperature region they are 3 orders of magnitude smaller.
Fig. 8. Difference between lgP obtained with B and Hsplines (a) over the whole SAHAS domain of definition, (b) for points (T, ρ) from the solar model. X = 0.80, Z = 0.02. 

Open with DEXTER 
Thus we notice that for pressure under solar conditions, the precision of B and Hspline interpolation is better than (2−3) × 10^{−6}, even in ionization regions. Therefore, the level of accuracy of SAHAS is sufficient for highprecision solar modeling.
To estimate the interpolation precision of pressure derivatives, we compare the calculation of Γ_{1} using two different methods. The first is with a Bspline that directly interpolates the discrete set of Γ_{1}. The other is based on Eq. (5), which requires the interpolation of the three values χ_{T}, χ_{ρ}, C_{Π}. The values χ_{T} and χ_{ρ} are calculated in accordance with Eq. (1) as derivatives of the Hermite spline of the pressure function:
The function C_{Π} in Eq. (5) was obtained with Bsplines interpolation.
Figure 9 shows the difference between these two methods. The adiabatic exponent is sensitive to the method of calculation. The maximum difference is on the order of 10^{−3} and it is attained in the temperature range of (1.5−300) × 10^{3} K. The precision is much better for the points of the solar model. The maximum difference of the adiabatic exponent is 2 × 10^{−4} in the lowtemperature ionization region, but it is smaller by two orders of magnitude in the hightemperatures region of the radiative core.
Fig. 9. Difference between the adiabatic exponents Γ_{1}, computed using B and Hsplines (a) over the whole SAHAS domain of definition, (b) for points (T, ρ) from the solar model. X = 0.80, Z = 0.02. 

Open with DEXTER 
Rapid oscillations of ΔΓ_{1} may be caused by different degrees of polynomials. In the case of Bsplines, Γ_{1} is represented by thirddegree polynomials, whereas the derivatives χ_{T} and χ_{ρ} from Eq. (5) are described the fourthdegree polynomials of the Hsplines.
4. Discussion
To demonstrate the significance of our numerical difference between B and Hsplines, we compare our result with the difference between the SAHAS and OPAL 2005 equations of state. For this, we have interpolated the SAHAS and OPAL 2005 tables to obtain values at the points (T,ρ) of the solar model. The interpolation of the OPAL 2005 tables is done by the original subroutine “esac” as provided with the OPAL tables. SAHAS tables are interpolated using Hspline method. The difference between the logarithms of pressure is shown in Fig. 10. It attains 1.3 × 10^{−3} and exceeds the difference due to the interpolation method by three orders of magnitude (Fig. 8b). The largest discrepancy between the two sets of tables is found at the temperature range of (0.3−1.0) × 10^{5} K and it is generally coincides with profiles of Coulomb correction parameter inside the Sun.
Fig. 10. Comparison of pressure in SAHAS and OPAL 2005 tables interpolated at points (T, ρ) from the solar model. X = 0.80, Z = 0.02. 

Open with DEXTER 
Figure 11 shows the difference in the adiabatic exponent Γ_{1}. The largest values, up to 0.032, are in the outer lowtemperature layer, at lgT < 3.7. This feature is an object of our next work. In the bulk of the solar model, the deviations do not exceed 2 × 10^{−3}, which is an order of magnitude higher than those coming from the different interpolation methods (Fig. 9b). Thus, the main factor of uncertainty in the equations of state is the difference in physics, and the influence of interpolation algorithms is an order of magnitude smaller.
Fig. 11. Comparison of adiabatic exponent Γ_{1} in SAHAS and OPAL 2005 tables interpolated at points (T, ρ) from the solar model. X = 0.80, Z = 0.02. 

Open with DEXTER 
Besides the interpolation method and the physics of an EOS, the choice of the table grid and its spacing will affect the resulting accuracy of the computations. This question has been examined, for instance, by Dorman et al. (1991) from the point of view of stellar modeling. Here, however, we have focussed on the interpolation method, while keeping the same table grid in the computations.
In the presentage Sun, the part with a lowtemperature (T < 10^{5} K) plasma is small, only 0.002% by mass. However, in the early stages of solar evolution, the opposite is the case. For example, at the beginning of the premain sequence stage (pMS), the lowtemperature region extends up to 90% of the solar mass. The two situations are compared in Fig. 12, which illustrates in both cases the difference between the adiabatic exponents Γ_{1} of SAHAS and OPAL 2005. The red line shows the case of the presentday Sun, and the blue line the pMS. The differences are shown as functions of the mass fraction m_{r}, where m_{r} = 0 corresponds to the center and m_{r} = 1 to the surface. Clearly, the deviations are several times greater in the early stage. The stellar models are more sensitive to the choice of the equation of state at an early stage of evolution than on the main sequence because the discrepancies in the equations of state are maximal in moderate and lowtemperature plasma.
Fig. 12. Difference in adiabatic exponent Γ_{1} in SAHAS and OPAL 2005 tables, interpolated at points (T, ρ) from the solar model on the pMS stage (blue line) and on the MS stage (red line). X = 0.80, Z = 0.02. 

Open with DEXTER 
5. Conclusion
The paper describes an algorithm for interpolating pressure in EOS tables using quintic Hermite splines. The advantage of the Hspline method over the standard Bspline interpolation of the different thermodynamic quantities is a physical and geometrical selfconsistency. This consistency requirement is built in by construction in the Hsplines. In other words, the Hspline interpolation operates both on the function and its derivatives and simultaneously preserves thermodynamic and differentialgeometric identities.
Another specific feature of the Hermite interpolation is locality, which allows effective calculation on irregular mesh, for example, in case of OPAL EOS, when the mesh is uneven and 2Dirregular. In other words, if one would have all necessary derivatives in knots of a cell, then computation of Hspline is irrelevant to global properties of the mesh (regular or irregular, even or not). In contrast, considered earlier by Baturin et al. (2017) multidimensional Bsplines can be constructed on Ndimensional regular mesh only. Principally, the specific pseudolocal algorithm of Ndimensional cubic spline exists and has been used in the esac procedure of OPAL tables. But uneven mesh could cause degradation of quality of computed derivatives.
A comparison of B and Hsplines leads to the following conclusions. (1) The precision of the interpolation methods is estimated: the difference in lgP is up to 10^{−4}, the difference in adiabatic exponent Γ_{1} is up to 10^{−3}; for the points taken from a solar model, these discrepancies are an order of magnitude smaller. (2) Two different regions are identified, with the accuracy being lowest in dissociation and ionization regions, T ∼ (1.5−100) × 10^{3} K, and highest at (T > 100 × 10^{3} K). (3) The equations of state SAHAS and OPAL 2005 are compared; the differences between them are shown to be some orders of magnitude greater than differences caused by the different interpolation methods. (4) Since the interpolation errors in Γ_{1} turn out to be significantly smaller than the differences between competing physical formalisms, the main issue is the choice of the equation of state.
References
 Baturin, V. A., Däppen, W., Morel, P., et al. 2017, A&A, 606, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brito, A., & Lopes, I. 2018, ApJ, 853, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Buldgen, G., Salmon, S. J. A. J., Noels, A., et al. 2019, A&A, 621, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ChristensenDalsgaard, J., & Däppen, W. 1992, A&ARv, 4, 267 [NASA ADS] [CrossRef] [Google Scholar]
 Däppen, W., Mihalas, D., Hummer, D. G., & Mihalas, B. W. 1988, ApJ, 332, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Dorman, B., Irwin, A. W., & Pedersen, B. B. 1991, ApJ, 381, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Eggleton, P. P., Faulkner, J., & Flannery, B. P. 1973, A&A, 23, 325 [NASA ADS] [Google Scholar]
 Finn, D. 2004, MA 323 Geometric Modelling. Course Notes: Day 09. Quintic Hermite Interpolation [Google Scholar]
 Gryaznov, V. K., Ayukov, S. V., Baturin, V. A., et al. 2006, J. Phys. A Math. Gen., 39, 4459 [NASA ADS] [CrossRef] [Google Scholar]
 Gryaznov, V. K., Iosilevskiy, I. L., Fortov, V. E., et al. 2013, Contrib. Plasma Phys., 53, 392 [NASA ADS] [CrossRef] [Google Scholar]
 Hsu, Y. L. 2010, Optimal Design Laboratory. Geometric Properties of Surfaces [Google Scholar]
 Hummer, D. G., & Mihalas, D. 1988, ApJ, 331, 794 [NASA ADS] [CrossRef] [Google Scholar]
 Irwin, A. W. 2012, Astrophysics Source Code Library [record ascl: 1211.002] [Google Scholar]
 Joyce, M., & Chaboyer, B. 2018, ApJ, 856, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Krasnikov, Y. G. 1977, Sov. J. Exp. Theor. Phys., 46, 270 [NASA ADS] [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1980, Statistical Physics. Pt.1, Pt.2 (Oxford: Pergamon Press) [Google Scholar]
 Mihalas, D., Däppen, W., & Hummer, D. G. 1988, ApJ, 331, 815 [NASA ADS] [CrossRef] [Google Scholar]
 Morel, P., & Lebreton, Y. 2008, Ap&SS, 316, 61 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Morel, P., Provost, J., & Berthomieu, G. 1997, A&A, 327, 349 [NASA ADS] [Google Scholar]
 Reichl, L. E. 1998, A Modern Course in Statistical Physics, 2nd edn. (WileyVCH) [Google Scholar]
 Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, F. J., Swenson, F. J., & Iglesias, C. A. 1996, ApJ, 456, 902 [NASA ADS] [CrossRef] [Google Scholar]
 Schultz, M. H., & Varga, R. S. 1967, Numer. Math., 10, 345 [CrossRef] [Google Scholar]
 Somers, G., & Pinsonneault, M. H. 2014, ApJ, 790, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Swesty, F. D. 1996, J. Comput. Phys., 127, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Varga, S. R. 1971, Functional Analysis and Approximational Theory in Numerical Analysis (Kent, Ohio: Kent State University) [CrossRef] [Google Scholar]
 Vorontsov, S. V., Baturin, V. A., Ayukov, S. V., & Gryaznov, V. K. 2013, MNRAS, 430, 1636 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Some definitions of the spline theory
Here we give some definitions of the standard (B)spline theory. We refer to the introduction of the book by Varga (1971) and restrict ourselves to the cases used in this paper. Piecewise polynomial splines of odd degree p = 2m − 1 are considered on the mesh {x_{i}} of the interval [a,b]:
In addition to the spline parameter m, another number is needed – the socalled spline defect d. The spline defect is equal to the number of discontinuous higher (nonzero) derivatives. Splines with defect d = 1 are referred to as natural. For example, a natural cubic spline (p = 3, m = 2) has continuous first and second derivatives while the third are discontinuous at the internal knots {x}^{in}. We denote natural splines by B_{3} or simply B, because only cubic splines of this type are considered in the paper. An example of Bsplines for the case of independent interpolations of the thermodynamic quantities is given by (Baturin et al. 2017).
By definition, splines have p − d = 2m − d − 1 continuous derivatives (including the zeroth order, that is, the function itself) at the internal knots {x}^{in} = {x_{1},…,x_{N−1}} of the mesh. Following Varga (1971), we define Hermite splines (Hsplines) by the relation d = m. In the case of the cubic splines, m = 2. Then the defect is 2, that is, the third and second spline derivatives have discontinuities at the internal knots. The first derivative and the spline function itself are continuous. We consider quintic splines H_{5}, in which d = m = 3. Thus, H_{5}splines have the first and second continuous derivatives, and all higher derivatives are discontinuous, beginning with the third.
A set of splines with parameters m, d on the mesh {x_{i}} (Eq. (A.1)) forms a linear space with dimension 2m + d(N−1). Schultz & Varga (1967) have shown that there is a single spline s, which interpolates a function f in such a linear space. That is, d values of the function and its derivatives are interpolated at the internal knots {x}^{in} of the mesh:
and m values at the boundary points a and b:
The following are direct consequences of Eqs. (A.2) and (A.3). Firstly, using Hermite splines allows us to solve the problem of interpolating a function together with its derivatives. This is necessary for the construction of consistent thermodynamic structures. Secondly, in the case d < m (in particular, in the case of natural splines), we must specify additional conditions on the derivatives at the boundaries of the interval. As a result, it becomes necessary to solve a system of equations consisting of continuity conditions and boundary conditions in order to find spline coefficients. Thirdly, in the case of Hermite splines, the interpolation conditions are the same at internal and boundary knots. We can take any interval as a boundary, for example [x_{i},x_{i + 1}]. Thus, construction of Hermite splines is local. Cubic Hermite splines H_{3} are completely determined for given values of the function and their first derivatives at the boundaries of an elementary interval. Similarly, quintic Hermitesplines H_{5} are defined on each interval if the function, first and second derivatives are given at the boundaries.
The transition from the onedimensional piecewisepolynomial interpolation to the twodimensional case is based on the principle of a regular 2D mesh {x,y}_{2D} of the form {x,y}_{2D} = {x_{i}} ⊗ {y_{i}}, that is, the rectangle of the domain of definition is separated by the coordinate lines (Varga 1971). The construction of the interpolating functions becomes the direct product of the interpolation in onedimensional spaces. This remains valid both for natural splines (see an example of constructing 3D splines in the paper by Baturin et al. 2017), and for Hermite splines.
The only thing that needs to be clarified is a system of partial derivatives that determine the continuity conditions (or specify the interpolated derivatives) for each elementary square. This system has the form
In other words, partial derivatives must be continuous (and given at the boundary knots in the case of Hsplines), for which each individual index of derivatives varies from zero to m − 1. Four functions , and f″_{xy} must be given for the cubic 2D Hermite spline, nine functions are needed for the quintic spline.
Appendix B: Derivatives in T, Q_{s}variables for calculating conditions in mesh knots
Nine values must be specified in each mesh knots to interpolate a 2D function by quintic Hermite splines. They are listed in Sect. 3.1, Eq. (16). For interpolation of lgP(lgT, lgQ_{s}), we used logarithmic derivatives χ_{T} and χ_{ρ}.
The higher derivatives are obtained from Eqs. (B.2) and (B.3) by interpolation of tabulated functions and using Bsplines and subsequent calculation of their derivatives
The second derivative of pressure with respect to temperature is calculated from Eq. (10), which is itself obtained from a thermodynamic identity.
where
The two remaining derivatives in Eq. (16) are calculated by differentiation of Eq. (B.7) with respect to lgQ_{s}
where
Similarly, the last derivative is
where
All Figures
Fig. 1. Differentialgeometric structure based on the Helmholtz free energy. 

Open with DEXTER  
In the text 
Fig. 2. Pressure P (a) and scaled pressure Π (b) in SAHAS EOS for X = 0.80 and Z = 0.02. Red lines show points (T, ρ) from the solar model. 

Open with DEXTER  
In the text 
Fig. 3. Logarithmic derivatives of pressure χ_{T} (a) and χ_{ρ} (b) in SAHAS EOS for X = 0.80 and Z = 0.02. Red lines indicate points (T, ρ) from the solar model. 

Open with DEXTER  
In the text 
Fig. 4. Specific heat capacity in units of scaled pressure in SAHAS EOS for X = 0.80 and Z = 0.02 (a) over the whole SAHAS domain of definition, the red line indicates points (T, ρ) from the solar model, (b) for points (T, ρ) from the solar model. 

Open with DEXTER  
In the text 
Fig. 5. Adiabatic exponent Γ_{1} in SAHAS EOS for X = 0.80 and Z = 0.02 (a) over the whole SAHAS domain of definition, the red line indicates points (T, ρ) from the solar model, (b) for points (T, ρ) from the solar model. 

Open with DEXTER  
In the text 
Fig. 6. Differences between the pressure derivatives calculated as the interpolator of the derivative and the derivative of the interpolator (Eqs. (11) and (12)) on the SAHAS mesh for X = 0.80 and Z = 0.02. Please note the difference of an order of two between the scales in the top and bottom panels. 

Open with DEXTER  
In the text 
Fig. 7. Difference between the second mixed derivatives of pressure (Eq. (14)) on the SAHAS mesh for X = 0.80 and Z = 0.02. 

Open with DEXTER  
In the text 
Fig. 8. Difference between lgP obtained with B and Hsplines (a) over the whole SAHAS domain of definition, (b) for points (T, ρ) from the solar model. X = 0.80, Z = 0.02. 

Open with DEXTER  
In the text 
Fig. 9. Difference between the adiabatic exponents Γ_{1}, computed using B and Hsplines (a) over the whole SAHAS domain of definition, (b) for points (T, ρ) from the solar model. X = 0.80, Z = 0.02. 

Open with DEXTER  
In the text 
Fig. 10. Comparison of pressure in SAHAS and OPAL 2005 tables interpolated at points (T, ρ) from the solar model. X = 0.80, Z = 0.02. 

Open with DEXTER  
In the text 
Fig. 11. Comparison of adiabatic exponent Γ_{1} in SAHAS and OPAL 2005 tables interpolated at points (T, ρ) from the solar model. X = 0.80, Z = 0.02. 

Open with DEXTER  
In the text 
Fig. 12. Difference in adiabatic exponent Γ_{1} in SAHAS and OPAL 2005 tables, interpolated at points (T, ρ) from the solar model on the pMS stage (blue line) and on the MS stage (red line). X = 0.80, Z = 0.02. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.