© 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; Christensen-Dalsgaard & 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), SAHA-S (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 sound-speed 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 EOS-induced 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 F-stars for asteroseismic analysis.

At the same time, fine helioseismic inversion indicates that the adiabatic exponent Γ₁ in the equation of state of solar plasma can be determined from observational data with an accuracy on the order of 10⁻⁴ (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 model-theory comparisons, which need high-quality stellar models. The equation-of-state part of these models must be a high-precision realization of the theory, irrespective of the accuracy of the theory itself. Keeping the interpolation errors in equation-of-state 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 ideal-gas 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, SAHA-S 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 piece-wise 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 trade-off 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 odd-degree splines. The first type is splines with maximal number of continuous derivatives, which are called natural or B-splines. The second type is Hermite splines, which have fewer continuous derivatives but more free parameters.

While using B-splines, 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 high-precision thermodynamic quantities needed in models of the solar-structure. 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 B-splines 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 differential-geometric identities. Section 3 describes interpolation method with Hermite splines and comparatively discusses its advantages and limitations. Section 4 compares differences between the SAHA-S 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 two-dimensional 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 differential-geometric structure of the potential is a manifold with its first-order partial derivatives on the first level and the second and third-order partial derivatives at the higher levels. Specifically, the case of the free energy F is illustrated in Fig. 1.

Fig. 1. Differential-geometric structure based on the Helmholtz free energy.

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 downward-pointing arrows in the figure) and to volume ∂/∂V  (left- and downward-pointing 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:

(1)

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

(2)

For a transition from the first to the third level, the normality condition leads to the identities

(3)

(4)

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 Γ₁ because of its fundamental role in stellar modeling. It can be calculated with the help of the partial derivatives of the potential as follows:

(5)

where the dimensionless heat capacity is introduced by

(6)

with the scaled pressure defined as

(7)

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, Γ₁, χ_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 2D-manifolds (surfaces) in coordinates T and Q_s = ρ · (T₆)^−2.25, where T₆ = T/10⁶. These variables are used in the SAHA-S tables and described in details by Baturin et al. (2017). We note that the table domain of SAHA-S 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³ K, the band where the hydrogen is half-ionized passes at (10−20) × 10³ K, and helium ionization regions lies at temperatures of (20−100) × 10³ 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 SAHA-S EOS for X = 0.80 and Z = 0.02. Red lines show points (T, ρ) from the solar model.

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 finite-difference methods in the mesh points, together with pressure. This is the main justification for the construction of a spline-manifold. The behavior of the logarithmic derivatives is shown in Fig. 3. Differentiation of P emphasizes all the physical features of the scaled-pressure surface. Conditions χ_T = χ_ρ = 1 are equivalent to the ideal-gas equation P = R_gTρ/μ 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 SAHA-S EOS for X = 0.80 and Z = 0.02. Red lines indicate points (T, ρ) from the solar model.

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 SAHA-S EOS for X = 0.80 and Z = 0.02 (a) over the whole SAHA-S domain of definition, the red line indicates points (T, ρ) from the solar model, (b) for points (T, ρ) from the solar model.

The adiabatic exponent Γ₁ is defined by Eq. (5) and it depends both on the pressure derivatives and heat capacity. We used Γ₁ to estimate the accuracy of the interpolation. The Γ₁ 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 SAHA-S EOS for X = 0.80 and Z = 0.02 (a) over the whole SAHA-S domain of definition, the red line indicates points (T, ρ) from the solar model, (b) for points (T, ρ) from the solar model.

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 ρ:

(8)

Thus, the second pressure derivative with respect to temperature can be expressed via χ_T and the derivative of c_V with respect to density:

(9)

For the sake of simplicity in the spline differentiation, we rewrite (9) with dimensionless function C_Π instead of c_V:

(10)

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 third-degree polynomial splines (i.e., B-spline) 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 B-spline 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 B-spline, 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 B-splines, 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 finite-difference 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 SAHA-S 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) ∂_TB[ln P̂] is not equal to . The corresponding differences

(11)

and

(12)

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 SAHA-S 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.

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⁻⁸−10⁻⁶.

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⁻⁴ (see Fig. 6). In the rest of the region, the difference varies from 10⁻¹⁰ to 10⁻⁷. 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 Hermite-spline construction.

The accuracy of the B-spline 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:

(13)

Such mixed derivatives are not presented in the up-to-date EOS-tables, so we need to use spline-differentiation 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

(14)

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⁻⁷−10⁻⁶ in the remaining regions (Fig. 7).

Fig. 7. Difference between the second mixed derivatives of pressure (Eq. (14)) on the SAHA-S mesh for X = 0.80 and Z = 0.02.

3. Quintic Hermite spline

3.1. Algorithm

The Hermite spline (H-spline) 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 H-spline. Thus, H-splines 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 H-spline 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 two-dimensional case). To calculate the H-spline inside the interval, one does not necessarily need to know the surrounding elements and the properties of the mesh around it. This property of H-splines prevents a disturbance due to the irregularity of the mesh or the boundary conditions.

One has to note, however, that H-splines 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 H-splines. In the two-dimensional 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 H-spline is represented in the form

(15)

Thirty-six 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:

(16)

The independent variables for 2D SAHA-S tables are

(17)

The interpolated function is the logarithm of pressure

(18)

Below we provide all expression for the derivatives of the matrix (16) via the tabulated functions and their spline-differentiation. Two first derivatives (16) are calculated using and :

(19)

(20)

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:

(21)

(22)

(23)

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): ∂²H/∂u², ∂³H/∂u²∂v and ∂⁴H/∂u²∂v² contain a double differentiation of lgP with respect to lgT. One way to calculate them is to differentiate B-splines 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

(24)

as well as the higher derivatives, which are obtained from it by differentiation with respect to lgQ_s:

(25)

(26)

The formulas for these derivatives in terms of (T,Q_s) are given in Appendix B.

3.2. Comparison of B- and H-splines

Figure 8 shows the difference in lgP, calculated with B- and H-splines. Panel a displays the data for the entire domain of applicability of SAHA-S. The greatest deviations take place in ionization and dissociation zones, where the discrepancy is on the order of 10⁻⁴. In the higher-temperature 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⁻⁶ in the low-temperature region while in the high-temperature region they are 3 orders of magnitude smaller.

Fig. 8. Difference between lgP obtained with B- and H-splines (a) over the whole SAHA-S domain of definition, (b) for points (T,  ρ) from the solar model. X = 0.80, Z = 0.02.

Thus we notice that for pressure under solar conditions, the precision of B- and H-spline interpolation is better than (2−3) × 10⁻⁶, even in ionization regions. Therefore, the level of accuracy of SAHA-S is sufficient for high-precision solar modeling.

To estimate the interpolation precision of pressure derivatives, we compare the calculation of Γ₁ using two different methods. The first is with a B-spline that directly interpolates the discrete set of Γ₁. 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:

(27)

The function C_Π in Eq. (5) was obtained with B-splines- 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⁻³ and it is attained in the temperature range of (1.5−300)  ×  10³ K. The precision is much better for the points of the solar model. The maximum difference of the adiabatic exponent is 2  ×  10⁻⁴ in the low-temperature ionization region, but it is smaller by two orders of magnitude in the high-temperatures region of the radiative core.

Fig. 9. Difference between the adiabatic exponents Γ1, computed using B- and H-splines (a) over the whole SAHA-S domain of definition, (b) for points (T, ρ) from the solar model. X = 0.80, Z = 0.02.

Rapid oscillations of ΔΓ₁ may be caused by different degrees of polynomials. In the case of B-splines, Γ₁ is represented by third-degree polynomials, whereas the derivatives χ_T and χ_ρ from Eq. (5) are described the fourth-degree polynomials of the H-splines.

4. Discussion

To demonstrate the significance of our numerical difference between B- and H-splines, we compare our result with the difference between the SAHA-S and OPAL 2005 equations of state. For this, we have interpolated the SAHA-S 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. SAHA-S tables are interpolated using H-spline method. The difference between the logarithms of pressure is shown in Fig. 10. It attains 1.3 × 10⁻³ 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⁵ K and it is generally coincides with profiles of Coulomb correction parameter inside the Sun.

Fig. 10. Comparison of pressure in SAHA-S and OPAL 2005 tables interpolated at points (T, ρ) from the solar model. X = 0.80, Z = 0.02.

Figure 11 shows the difference in the adiabatic exponent Γ₁. The largest values, up to 0.032, are in the outer low-temperature 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⁻³, 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 SAHA-S and OPAL 2005 tables interpolated at points (T, ρ) from the solar model. X = 0.80, Z = 0.02.

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 present-age Sun, the part with a low-temperature (T <  10⁵ 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 pre-main sequence stage (pMS), the low-temperature 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 Γ₁ of SAHA-S and OPAL 2005. The red line shows the case of the present-day 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 low-temperature plasma.

Fig. 12. Difference in adiabatic exponent Γ1 in SAHA-S 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.

5. Conclusion

The paper describes an algorithm for interpolating pressure in EOS tables using quintic Hermite splines. The advantage of the H-spline method over the standard B-spline interpolation of the different thermodynamic quantities is a physical and geometrical self-consistency. This consistency requirement is built in by construction in the H-splines. In other words, the H-spline interpolation operates both on the function and its derivatives and simultaneously preserves thermodynamic and differential-geometric 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 2D-irregular. In other words, if one would have all necessary derivatives in knots of a cell, then computation of H-spline is irrelevant to global properties of the mesh (regular or irregular, even or not). In contrast, considered earlier by Baturin et al. (2017) multidimensional B-splines can be constructed on N-dimensional regular mesh only. Principally, the specific pseudo-local algorithm of N-dimensional 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 H-splines leads to the following conclusions. (1) The precision of the interpolation methods is estimated: the difference in lgP is up to 10⁻⁴, the difference in adiabatic exponent Γ₁ is up to 10⁻³; 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³ K, and highest at (T >  100 × 10³ K). (3) The equations of state SAHA-S 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 Γ₁ 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

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]:

(A.1)

In addition to the spline parameter m, another number is needed – the so-called spline defect d. The spline defect is equal to the number of discontinuous higher (non-zero) 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}ⁱⁿ. We denote natural splines by B₃ or simply B, because only cubic splines of this type are considered in the paper. An example of B-splines 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}ⁱⁿ = {x₁,…,x_N−1} of the mesh. Following Varga (1971), we define Hermite splines (H-splines) 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₅, in which d = m = 3. Thus, H₅-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}ⁱⁿ of the mesh:

(A.2)

and m values at the boundary points a and b:

(A.3)

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₃ are completely determined for given values of the function and their first derivatives at the boundaries of an elementary interval. Similarly, quintic Hermite-splines H₅ are defined on each interval if the function, first and second derivatives are given at the boundaries.

The transition from the one-dimensional piecewise-polynomial interpolation to the two-dimensional 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 one-dimensional 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

(A.4)

In other words, partial derivatives must be continuous (and given at the boundary knots in the case of H-splines), 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 χ_ρ.

(B.1)

(B.2)

(B.3)

The higher derivatives are obtained from Eqs. (B.2) and (B.3) by interpolation of tabulated functions and using B-splines and subsequent calculation of their derivatives

(B.4)

(B.5)

(B.6)

The second derivative of pressure with respect to temperature is calculated from Eq. (10), which is itself obtained from a thermodynamic identity.

(B.7)

where

(B.8)

The two remaining derivatives in Eq. (16) are calculated by differentiation of Eq. (B.7) with respect to lgQ_s

(B.9)

where

(B.10)

Similarly, the last derivative is

(B.11)

where

(B.12)

All Figures

	Fig. 1. Differential-geometric structure based on the Helmholtz free energy.
In the text

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

	Fig. 3. Logarithmic derivatives of pressure χT (a) and χρ (b) in SAHA-S EOS for X = 0.80 and Z = 0.02. Red lines indicate points (T, ρ) from the solar model.
In the text

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

	Fig. 5. Adiabatic exponent Γ1 in SAHA-S EOS for X = 0.80 and Z = 0.02 (a) over the whole SAHA-S domain of definition, the red line indicates points (T, ρ) from the solar model, (b) for points (T, ρ) from the solar model.
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 SAHA-S 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.
In the text

	Fig. 7. Difference between the second mixed derivatives of pressure (Eq. (14)) on the SAHA-S mesh for X = 0.80 and Z = 0.02.
In the text

	Fig. 8. Difference between lgP obtained with B- and H-splines (a) over the whole SAHA-S domain of definition, (b) for points (T,  ρ) from the solar model. X = 0.80, Z = 0.02.
In the text

	Fig. 9. Difference between the adiabatic exponents Γ1, computed using B- and H-splines (a) over the whole SAHA-S domain of definition, (b) for points (T, ρ) from the solar model. X = 0.80, Z = 0.02.
In the text

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

	Fig. 11. Comparison of adiabatic exponent Γ1 in SAHA-S and OPAL 2005 tables interpolated at points (T, ρ) from the solar model. X = 0.80, Z = 0.02.
In the text

	Fig. 12. Difference in adiabatic exponent Γ1 in SAHA-S 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.
In the text