Interpolation of equation-of-state data

TL;DR

This paper introduces a Hermite spline interpolation method for equation-of-state data that preserves thermodynamic relations and estimates interpolation accuracy, demonstrating high precision in pressure and derivative calculations across a temperature-density grid.

Contribution

The authors develop a Hermite spline interpolation technique for EOS data that maintains thermodynamic identities and provides error estimates, improving upon traditional spline methods.

Findings

Maximum pressure difference is between 10^{-9} and 10^{-4}

Interpolation errors are much smaller than differences between EOS models

Best precision at temperatures above 10^5 K

Abstract

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 Hermite-spline results with those of frequent cubic (B-) spline interpolation. Methods. We have interpolated pressure as a function of temperature and density with quintic Hermite 2D-splines. 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 B-splines. Results. The results were obtained with the grid of the SAHA-S equation-of-state (EOS) tables. The maximum $lg P$ difference lies in the range from $10^{-9}$ to $10^{-4}$, and $\Gamma_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 \sim 1.5\cdot 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 equation-of-state formalisms, SAHA-S 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 solar-model points.

Equations

$${\chi }_T={\frac{\partial \ln P}{\partial \ln T}|}_{\rho },{\chi }_{\rho }={\frac{\partial \ln P}{\partial \ln \rho }|}_T,c_V=T{\frac{\partial S}{\partial T}|}_{\rho }.$$

$${\frac{\partial P}{\partial T}|}_V={\frac{\partial S}{\partial V}|}_T.$$

$$\frac{\partial }{\partial T}\left(\frac{P}{\rho }{\chi }_{\rho }\right)=\frac{\partial }{\partial \rho }\left(\frac{P}{T}{\chi }_T\right),$$

$$T{\left(\frac{{\partial }^{2}P}{\partial T^2}\right)|}_V={\left(\frac{\partial c_V}{\partial V}\right)|}_T.$$

$${\Gamma }_1={\frac{\partial \ln P}{\partial \ln \rho }|}_S={\chi }_{\rho }+\frac{{\chi }_T^2}{C_{\Pi }},$$

$$C_{\Pi }=\frac{c_V}{\Pi },$$

$$\Pi =\frac{P}{T\rho }.$$

$$V={\rho }^{-1},dV=-{\rho }^{-2}d\rho .$$

$${\frac{{\partial }^2\ln P}{{\partial \left(\ln T\right)}^2}|}_{\rho \rho }={\chi }_T\left(1-{\chi }_T\right)-\frac{1}{\Pi }{\frac{\partial c_V}{\partial \ln \rho }|}_T.$$

$${\frac{{\partial }^2\ln P}{\partial {(\ln T)}^2}|}_{\rho \rho }={\chi }_T\left(1-{\chi }_T\right)-{\frac{\partial C_{\Pi }}{\partial \ln \rho }|}_T+\left({\chi }_{\rho }-1\right)C_{\Pi }.$$

$${\Delta }_T=B\left[{\hat{\chi }}_T\right]-{\frac{\partial }{\partial \ln T}|}_{\rho }B\left[\ln \hat{P}\right]$$

$${\Delta }_{\rho }=B\left[{\hat{\chi }}_{\rho }\right]-{\frac{\partial }{\partial \ln \rho }|}_{T}B\left[\ln \hat{P}\right]$$

$$\frac{{\partial }^{2}f}{\partial x\partial y}=\frac{{\partial }^{2}f}{\partial y\partial x}.$$

$$\frac{{\partial }^2\ln P}{\partial \ln T\partial \ln \rho }.$$

$${\Delta }_{\left\{T,\rho \right\}}={\frac{\partial }{\partial \ln T}|}_{\rho }B\left[{\hat{\chi }}_{\rho }\right]-{\frac{\partial }{\partial \ln \rho }|}_{T}B\left[{\hat{\chi }}_T\right]$$

$$H\left(u,v\right)=\sum _{i=0}^5\sum _{j=0}^5{a}_{ij}{u}^i{v}^j.$$

$$\left[\begin{array}{ccc}H& \frac{\partial H}{\partial u}& \frac{{\partial }^{2}H}{\partial u^2}\\ \frac{\partial H}{\partial v}& \frac{{\partial }^{2}H}{\partial u\partial v}& \frac{{\partial }^{3}H}{\partial u^2\partial v}\\ \frac{{\partial }^{2}H}{\partial v^2}& \frac{{\partial }^{3}H}{\partial u\partial v^2}& \frac{{\partial }^{4}H}{\partial u^2\partial v^2}\end{array}\right].$$

$$u=\lg T,v=\lg Q_s.$$

$$H(u,v)=\lg P(\lg T,\lg Q_s).$$

$$\frac{\partial H}{\partial u}{=\frac{\partial \lg P}{\partial \lg T}|}_{Q_s}={\hat{\chi }}_T+2.25{\hat{\chi }}_{\rho },$$

$$\frac{\partial H}{\partial v}={\frac{\partial \lg P}{\partial \lg Q_s}|}_T={\hat{\chi }}_{\rho }.$$

$$\frac{{\partial }^{2}H}{\partial u\partial v}=\frac{\partial }{\partial v}\left(\frac{\partial H}{\partial u}\right)={\frac{\partial }{\partial \lg Q_s}|}_T\left({\frac{\partial \lg P}{\partial \lg T}|}_{Q_s}\right),$$

$$\frac{{\partial }^{3}H}{\partial u\partial v^2}=\frac{{\partial }^2}{\partial v^2}\left(\frac{\partial H}{\partial u}\right)={\frac{{\partial }^2}{\partial {(\lg Q_s)}^2}|}_T\left({\frac{\partial \lg P}{\partial \lg T}|}_{Q_s}\right),$$

$$\frac{{\partial }^{2}H}{\partial v^2}=\frac{\partial }{\partial v}\left(\frac{\partial H}{\partial v}\right)={\frac{\partial }{\partial \lg Q_s}|}_T\left({\frac{\partial \lg P}{\partial \lg Q_s}|}_T\right).$$

$$\frac{{\partial }^{2}H}{\partial u^2}={\frac{{\partial }^2\lg P}{\partial {(\lg T)}^2}|}_{Q_s{Q}_s}$$

$$\frac{{\partial }^{3}H}{\partial u^2\partial v}=\frac{\partial }{\partial v}\left(\frac{{\partial }^{2}H}{\partial u^2}\right)={\frac{\partial }{\partial \lg Q_s}|}_T\left({\frac{{\partial }^2\lg P}{\partial {(\lg T)}^2}|}_{Q_s{Q}_s}\right),$$

$$\frac{{\partial }^{4}H}{\partial u^2\partial v^2}=\frac{{\partial }^2}{\partial v^2}\left(\frac{{\partial }^{2}H}{\partial u^2}\right)={\frac{{\partial }^2}{\partial {(\lg Q_s)}^2}|}_T\left({\frac{{\partial }^2\lg P}{\partial {(\lg T)}^2}|}_{Q_s{Q}_s}\right).$$

$${\chi }_{\rho }={\frac{\partial }{\partial \lg Q_s}|}_{T}H\left[\lg P\right],{\chi }_T={\frac{\partial }{\partial \lg T}|}_{Q_s}H\left[\lg P\right]-2.25{\chi }_{\rho }.$$

$$\left\{x_i\right\}:a=x_0<x_1<...<x_N=b.$$

$$D^j\left(f-s\right)\left({\left\{x\right\}}^{\text{in}}\right)=0,0\le j\le d-1,$$

Full text

Methods.

Results.

*Methods. *We have interpolated pressure as a function of temperature and density with quintic Hermite 2D-splines. 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 B-splines.

*Results. * The results were obtained with the grid of the SAHA-S equation-of-state (EOS) tables. The maximum $\lg P$ difference lies in the range from $10^{-9}$ to $10^{-4}$, and $\Gamma_{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\sim 1.5\cdot 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 equation-of-state formalisms, SAHA-S 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 solar-model points.

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 Msun. 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 $\Gamma_{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 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\left(T,\rho\right)$ and entropy $S(T,\rho)$, are required, each of which can be represented by a two-dimensional manifold of temperature $T$ and density $\rho$. 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\left(T,V\right)$. 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.

Each higher level of the structure of the potential in Fig. 1 is obtained as a result of differentiation with respect to temperature ${\partial}/{\partial T}\;$ (designated by right- and downward-pointing arrows in the figure) and to volume ${\partial}/{\partial 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 ${{\chi}_{T}}$, ${{\chi}_{\rho}}$, and the specific heat capacity $c_{V}$ on the second level:

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

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 $\Gamma_{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:

[TABLE]

where the dimensionless heat capacity is introduced by

[TABLE]

with the scaled pressure defined as

[TABLE]

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\left(T,\rho\right)$ 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$, ${{\Gamma}_{1}}$, ${{\chi}_{T}}$, ${{\chi}_{\rho}}$, ${{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}}=\rho\cdot{{\left({{T}_{6}}\right)}^{-2.25}}$, where $T_{6}=T/10^{6}$. 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 $\Pi$ is much more revealing (Fig. 2b). In the case of the ideal gas law, the scaled pressure $\Pi$ 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\cdot{{10}^{3}}$ K, the band where the hydrogen is half-ionized passes at $(10-20)\cdot{{10}^{3}}$K, and helium ionization regions lies at temperatures of $(20-100)\cdot{{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.

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 ${{\chi}_{T}}={{\chi}_{\rho}}=1$ are equivalent to the ideal-gas equation $P=R_{\rm g}T\rho/\mu$ with constant molecular weight $\mu$. $R_{\rm 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 ${{\chi}_{T}}$ and ${{\chi}_{\rho}}$ as well as corresponding effect of electron degeneracy.

Besides pressure $P$ as a function of $T$ and $\rho$, 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 $\Pi$ (see Eq. (6)). The behavior of ${{C}_{\Pi}}$ resembles that of ${{\chi}_{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 ${{\chi}_{T}}$. Second, electron degeneracy does not significantly affect the ${{C}_{\Pi}}$ surface.

The adiabatic exponent ${{\Gamma}_{1}}$ is defined by Eq. (5) and it depends both on the pressure derivatives and heat capacity. We used ${{\Gamma}_{1}}$ to estimate the accuracy of the interpolation. The ${{\Gamma}_{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.

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), § 16, equation (16.1)) may be rewritten by a transformation from volume $V$ to density $\rho$:

[TABLE]

Thus, the second pressure derivative with respect to temperature can be expressed via ${{\chi}_{T}}$ and the derivative of ${{c}_{V}}$ with respect to density:

[TABLE]

For the sake of simplicity in the spline differentiation, we rewrite (9) with dimensionless function ${{C}_{\Pi}}$ instead of ${{c}_{V}}$:

[TABLE]

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\left[{\hat{f}}\right]$ for this B-spline interpolation of a function $f$. The cap over the function $\hat{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 $\partial B\left[{\hat{f}}\right]$ differs from the true derivative of the function ${f}^{\prime}$. Moreover, the difference ${{\Delta}_{\partial B}}={f}^{\prime}-\partial B\left[{\hat{f}}\right]$ is not zero even at the nodes. In practice, we did not have an exact analytical estimation of ${f}^{\prime}$; but instead, the derivatives $\hat{{f}^{\prime}}$ are obtained by finite-difference method at the stage of calculating the EOS tables. Their spline interpolation is denoted by $B\left[{\hat{{f}^{\prime}}}\right]$. The difference $\Delta=B\left[{\hat{{f}^{\prime}}}\right]-\partial B\left[{\hat{f}}\right]$ can be calculated and used to estimate accuracy of spline differentiation $\partial B\left[{\hat{f}}\right]$.

We illustrate this with the example of SAHA-S tables. The tabulated pressure $\ln\hat{P}$ and its derivatives ${{\hat{\chi}}_{T}}$, ${{\hat{\chi}}_{\rho}}$ are interpolated with $B\left[\ln\hat{P}\right]$, $B\left[{{{\hat{\chi}}}_{T}}\right]$, and $B\left[{{{\hat{\chi}}}_{\rho}}\right]$ (Baturin et al. 2017). As a result, the spline derivative (with respect to temperature, for example) ${{\partial}_{T}}B\left[\ln\hat{P}\right]$ is not equal to $B\left[{{{\hat{\chi}}}_{T}}\right]$. The corresponding differences

[TABLE]

and

[TABLE]

are shown in Fig. 6.

The greatest differences of ${{\Delta}_{T}}$ are in the ionization and dissociation regions, since the functions $\ln\hat{P}$ and ${{\hat{\chi}}_{T}}$ 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 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\left(x,y\right)$ (smooth with second derivatives), the matching condition for mixed derivatives must be fulfilled:

[TABLE]

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

[TABLE]

The first is by differentiating $B\left[{{{\hat{\chi}}}_{\rho}}\right]$ with respect to temperature, and the second one by differentiating $B\left[{{{\hat{\chi}}}_{T}}\right]$ with respect to density. According to Eq. (13), the difference

[TABLE]

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).

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. (29) and (30) 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 (29) is equal to that of (30), 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

[TABLE]

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:

[TABLE]

The independent variables for 2D SAHA-S tables are

[TABLE]

The interpolated function is the logarithm of pressure

[TABLE]

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 ${{\hat{\chi}}_{T}}$ and ${{\hat{\chi}}_{\rho}}$:

[TABLE]

[TABLE]

The computation of mixed and higher derivatives in (16) is based on the spline differentiation of $B\left[{{{\hat{\chi}}}_{T}}\right]$ and $B\left[{{{\hat{\chi}}}_{\rho}}\right]$. The functions in the $\lg{{Q}_{s}}$ direction are smoother than along $\lg T$ (Fig. 3). So, differentiation with respect to $\lg{{Q}_{s}}$ is preferable whenever it is possible:

[TABLE]

[TABLE]

[TABLE]

Detailed formulae that explicitly relate to the higher derivatives of pressure with ${{\chi}_{T}}$ and ${{\chi}_{\rho}}$ are given in Appendix B.

The three derivatives in the third column of (16): ${{\partial}^{2}}H/\partial{{u}^{2}}$, ${{\partial}^{3}}H/\partial{{u}^{2}}\partial v$ and ${{\partial}^{4}}H/\partial{{u}^{2}}\partial{{v}^{2}}$ contain a double differentiation of $\lg P$ with respect to $\lg T$. One way to calculate them is to differentiate B-splines of ${{\hat{\chi}}_{T}}$ and ${{\hat{\chi}}_{\rho}}$ with respect to $\lg T$, but this leads to an inaccurate result, because ${{\hat{\chi}}_{T}}$ and ${{\hat{\chi}}_{\rho}}$ vary greatly with temperature. A better way is to avoid differentiation with respect to $\lg T$ altogether by using ${{c}_{V}}$ and ${{\chi}_{T}}$, which are derived from the thermodynamic relations in Sect. 2.3. Equation (10) allows one to estimate the second derivative

[TABLE]

as well as the higher derivatives, which are obtained from it by differentiation with respect to $\lg{{Q}_{s}}$:

[TABLE]

[TABLE]

The formulas for these derivatives in terms of $\left(T,{{Q}_{s}}\right)$ are given in Appendix B.

3.2 Comparison of B- and H-splines

Figure 8 shows the difference in $\lg P$ , 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}^{-4}}$. In the higher-temperature regions, they are five orders of magnitude lower. The red line shows the conditions $(T,\rho)$ of the solar model. The deviations with respect to these points are also shown in panel (b). They do not exceed $2.5\cdot{{10}^{-6}}$ in the low-temperature region while in the high-temperature region they are 3 orders of magnitude smaller.

Thus we notice that for pressure under solar conditions, the precision of B- and H- spline interpolation is better than $(2-3)\cdot{{10}^{-6}}$, 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 ${{\Gamma}_{1}}$ using two different methods. The first is with a B-spline that directly interpolates the discrete set of ${{\Gamma}_{1}}$. The other is based on Eq. (5), which requires the interpolation of the three values ${{\chi}_{T}}$, ${{\chi}_{\rho}}$, ${{C}_{\Pi}}$. The values ${{\chi}_{T}}$ and ${{\chi}_{\rho}}$ are calculated in accordance with Eq. (1) as derivatives of the Hermite spline of the pressure function:

[TABLE]

The function ${{C}_{\Pi}}$ 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}^{-3}}$ and it is attained in the temperature range of $(1.5-300)\cdot{{10}^{3}}\text{ K}$. The precision is much better for the points of the solar model. The maximum difference of the adiabatic exponent is $2\cdot{{10}^{-4}}$ in the low-temperature ionization region, but it is smaller by two orders of magnitude in the high-temperatures region of the radiative core.

Rapid oscillations of $\Delta{{\Gamma}_{1}}$ may be caused by different degrees of polynomials. In the case of B-splines, ${{\Gamma}_{1}}$ is represented by third-degree polynomials, whereas the derivatives ${{\chi}_{T}}$ and ${{\chi}_{\rho}}$ 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 $\left(T,\rho\right)$ 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\cdot{{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)\cdot{{10}^{5}}\text{ K}$ and it is generally coincides with profiles of Coulomb correction parameter inside the Sun.

Figure 11 shows the difference in the adiabatic exponent ${{\Gamma}_{1}}$. The largest values, up to 0.032, are in the outer low-temperature layer, at $\lg T<3.7$. This feature is an object of our next work. In the bulk of the solar model, the deviations do not exceed $2\cdot{{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.

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}^{5}}\text{ 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 Figure 12, which illustrates in both cases the difference between the adiabatic exponents ${{\Gamma}_{1}}$ 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.

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 $\lg P$ is up to ${{10}^{-4}}$, the difference in adiabatic exponent ${{\Gamma}_{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\sim(1.5-100)\cdot{{10}^{3}}\text{ K}$, and highest at ($T>100\cdot{{10}^{3}}\text{ 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 ${{\Gamma}_{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.

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 $\left\{{{x}_{i}}\right\}$ of the interval $\left[a,b\right]$:

[TABLE]

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,\text{ }m=2$) has continuous first and second derivatives while the third are discontinuous at the internal knots ${{\left\{x\right\}}^{\text{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 $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 ${{\left\{x\right\}}^{\text{in}}}=\left\{{{x}_{1}},...,{{x}_{N-1}}\right\}$ 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}_{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,\text{ }d$ on the mesh $\left\{{{x}_{i}}\right\}$ (Eq. (28)) forms a linear space with dimension $2m+d\left(N-1\right)$. 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 ${{\left\{x\right\}}^{\text{in}}}$ of the mesh:

[TABLE]

and $m$ values at the boundary points $a$ and $b$:

[TABLE]

The following are direct consequences of Eqs. (29) and (30). 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 $\left[{{x}_{i}},{{x}_{i+1}}\right]$. 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 Hermite-splines ${{H}_{5}}$ 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 ${{\left\{x,y\right\}}_{2D}}$ of the form ${{\left\{x,y\right\}}_{2D}}=\left\{{{x}_{i}}\right\}\otimes\left\{{{y}_{i}}\right\}$ , 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

[TABLE]

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 $f,{{{f}^{\prime}}_{x}},{{{f}^{\prime}}_{y}},$ and ${{{f}^{\prime\prime}}_{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 $\lg P(\lg T,\lg{{Q}_{s}})$, we used logarithmic derivatives ${{\chi}_{T}}$ and ${{\chi}_{\rho}}$ .

[TABLE]

[TABLE]

[TABLE]

The higher derivatives are obtained from Eqs. (33) and (34) by interpolation of tabulated functions ${{\hat{\chi}}_{T}}$ and ${{\hat{\chi}}_{\rho}}$ using B-splines and subsequent calculation of their derivatives

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

where

[TABLE]

The two remaining derivatives in Eq. (16) are calculated by differentiation of Eq. (38) with respect to $\lg Q_{s}$

[TABLE]

where

[TABLE]

Similarly, the last derivative is

[TABLE]

where

[TABLE]