Equation of state for solar near-surface convection

Numerical 3-D radiative hydrodynamical simulations are the main tool for the analysis of the interface between the solar convection zone and the photosphere. The equation of state is one of the necessary ingredients of these simulations. We compare two equations of state that are commonly used, one ideal and one nonideal, and quantify their differences. Using a numerical code we explore how these differences propagate with time in a 2-D convection simulation. We show that the runs with different equations of state (EOSs) and everything else identical relax to statistically steady states in which the mean temperature (in the range of the continuum optical depths typical for the solar photosphere) differs by less than 0.2 %. For most applications this difference may be considered insignificant.


Introduction
Realistic time-dependent 3-D numerical simulations (Nordlund, 1982;Stein andNordlund, 1998, 2000) radically improved our knowledge of the near-surface convection and the solar photosphere.In these simulations the equations of magnetohydrodynamics, the equation of state (EOS) and the radiative transfer equation (RTE) are solved ab initio with only few free parameters.Apart form the three-dimensionality, the key ingredient that makes these simulations realistic is the detailed description of the microphysical processes contributing to the EOS and the RTE.The realism of the simulations had been repeatedly confirmed by diverse numerical experiments comparing the actual observations of the sun with the spectra computed a posteriori from the models (As-plund et al., 2000b;Khomenko et al., 2005;Shelyag et al., 2007;Danilovic et al., 2008;Wedemeyer-Böhm and Rouppe van der Voort, 2009;Pereira et al., 2013).Beeck et al. (2012) compared the properties of the 3-D hydrodynamical (HD) simulations carried out by three codes that are the workhorses in the field: Stagger (Galsgaard and Nordlund, 1996), MPS/University of Chicago Radiative Magnetohydrodynamics (MURaM) (Vögler et al., 2005) and COnservative COde for the COmputation of COmpressible COnvection in a BOx of L Dimensions, L = 2, 3 (Co5Bold) (Wedemeyer et al., 2004).The codes are all primarily designed for simulations of near-surface convection, but with some intrinsic differences in the details of the microphysics, in the numerical algorithms, the boundary conditions, the parallelisation strategy, etc.Nevertheless, Beeck et al. (2012) demonstrated that the mean properties and the spatial distribution of the fundamental quantities in the HD simulations performed with these codes do not vary significantly.The difference in the results of these simulations is attributed to the difference in the input parameters (the opacities, the chemical composition).However, because the outcome of the code depends on these parameters in a highly non-linear way, it is impracticable to trace back the differences to the initial assumptions and approximations.An independent code was developed more recently by Tanner et al. (2012, the RHD code).Although the RTE in this code is solved in a simplified way, the authors showed that their results of the convection simulation are comparable to the results of Beeck et al. (2012).Tanner et al. (2013) employed the radiative hydrodynamics (RHD) code to study the variation in the simulated stellar convection with varying input metallicity.They concluded that the variation in the superadiabacity and the convection dynamics in their simulation remains relatively small for a large range of metallicities (0.01 ≤ Z ≤ 0.40).
Published by Copernicus Publications on behalf of the European Geosciences Union.
Figure 1.The boundaries in the T -ρ plane at which different nonideal effects start to be significant.The boundaries are estimated for the case of the pure hydrogen plasma (after Hansen et al., 2004).The T -ρ dependency of the 1-D models of the solar atmosphere (VALC model, red curve) and of the solar convection zone (after Spruit (1974), blue) is indicated.The shaded area shows the portion of the T -ρ plane occupied by a typical snapshot from a 3-D HD simulation.
While both studies clearly indicate that the results of the different realistic simulations provide a consistent qualitative and quantitative description of the solar convection, the differences between different simulation runs are still present.Moreover, these differences may become significant in certain experiments and measurements.A good example is the problem of the abundance measurement from the comparison of the synthetic and observed spectra.This problem is extremely model-dependent, and high accuracy is critical.
In this paper we isolate the influence of the EOS on the hydrodynamical convection simulations.In Sect.2, we quantitatively compare an ideal and a nonideal equilibrium EOS.They represent the main two types of EOS used in the realistic simulations.A simple numerical experiment designed to study the propagation of the differences caused by different EOSs through the simulated domain and in time is described in Sect.3.

Ideal versus nonideal EOS
Two types of EOS are used in the codes for the solar convection simulations: (1) the nonideal EOS including the effects of the pressure ionisation, Coulomb interaction and electron degeneracy and (2) the ideal EOS for a mixture including the partial ionisation effects.MHD1 (Hummer and Mihalas, 1988;Mihalas et al., 1988) and OPAL (Rogers et al., 1996;Rogers and Nayfonov, 2002) are commonly used nonideal EOSs.The former is used in Stagger, the latter in MURaM, RHD, in the ANTARES code (Grimm-Strele et al., 2015) and in the MANCHA code (Khomenko et al., 2014).The ideal EOS of Wolf (1983) is implemented in Co5Bold.MU-RaM can work with a simple ideal EOS described in Vögler et al. (2005), while MANCHA incorporates an algorithm for the evaluation of the ideal EOS of a partially ionised mixture described by Vardya (1965) and improved by Mihalas (1967) and Wittmann (1974) (we shall refer to this EOS as VMW).The VMW algorithm is implemented in several radiative transfer codes to provide the electron pressure for known temperature and gas pressure (e.g.Bellot Rubio et al., 1999;Socas-Navarro, 2011).
The range of validity of the ideal EOS may be estimated for the trivial case of partially ionised pure hydrogen.In Fig. 1, we reproduce Fig. 3.9 of Hansen et al. (2004), showing the T -ρ plane with the boundaries for the nonideal effects and the curve where the ionisation fraction of hydrogen is 0.99.To Fig. 1 we add plots of two 1-D stratifications (the chromospheric VALC model of Vernazza et al. (1981) and the convection zone model of Spruit, 1974) and the area covered by a typical 3-D HD simulation (shaded area).All three models are within the limits of where the ideal EOS for pure hydrogen is valid.The deep end of the convection-zone model is close to the pressure-ionisation boundary, while, on the other end, the ionisation of hydrogen in the chromosphere requires a non-equilibrium solution of hydrogen ionisation (Carlsson and Stein, 2002) that cannot be represented by a simple curve in this plot.Nevertheless, the near-surface convection and the photosphere fall in the region of the T -ρ plane, where the ideal EOS is a safe assumption.
To compute nonideal EOS is extremely time-consuming and impossible to do on the fly in codes for 3-D numerical simulations.Instead EOS is precomputed and results are stored in lookup tables suitable for fast interpolation.The computational cost of the ideal EOS is much lower; however, when the molecules are taken into account, it is necessary to solve the equations iteratively, and, therefore, the use of the precompiled lookup tables may save considerable computing time as long as the table grid is sufficiently fine to limit the interpolation errors.

Ideal EOS for partially ionised mixture
The VMW EOS accounts for the partial ionisation of the mixture composed of atomic and molecular hydrogen (H, H−, H+, H 2 and H + 2 ), neutral, singly and doubly ionised atoms of helium and other elements.The ionisation equilibrium is computed using Saha's equation, while for the dissociation of the molecules, the instantaneous chemical equilibrium is assumed.We implemented the algorithm in a stand-alone code with several improvements, mainly related to the atomic data.For the atomic partition functions, we adopted polynomial fits of Irwin (1981).The molecular partition function data for H 2 and H + 2 come from the survey of Sauval and Tatum (1984).The chemical equilibrium constants and the energies of the rotational-vibrational bound states are evaluated using these partition functions.Internal energy is defined up to an arbitrary constant.We follow the choice made in the OPAL and MHD EOS and set the potential component of the internal energy to 0 when all hydrogen atoms are in the ground state of the H 2 molecule and all other atoms are neutral and in their atomic ground states.For the chemical composition we adopt the high-metallicity abundances of Anders and Grevesse (1989).Since our prime interest here is differential analysis, the choice of the abundance set is not critical.On the other hand, the set of Anders and Grevesse is built into the OPAL EOS and cannot be altered once the tables are produced.We do not compare the VMW directly to other ideal EOSs as they are based on the very similar approach and the differences among them are minor as long as the same chemical composition is selected.
We evaluated the VMW EOS for a large T -p grid covering 17 orders of magnitude in pressure and nearly 6 in temperature.This range of values matches the range of the OPAL EOS tables.Figure 2 shows the ionisation fraction of hydrogen, x H = n H /n tot H , with n H and n tot H being the number density of the neutral hydrogen and the total number density of the hydrogen atoms and molecules.Figure 3 shows the distribution in the log T -log p plane of the mass density, the electron number density and the specific internal energy per mass 2 .

Nonideal EOS: OPAL
The OPAL EOS is computed in the physical picture through an activity expansion of the grand canonical partition function of the plasma.The results are distributed as lookup ta-2 Hereafter we use "internal energy" for the specific internal energy per mass.bles3 of the internal energy, the pressure, the electron density and other derived thermodynamical quantities as functions of the temperature and the density of several chemical compositions.We reversely interpolate the tables (for the mass fraction X = 0.7062 and Z = 0.0197, corresponding to the abundances of Anders and Grevesse) to obtain the internal energy, the density and the electron pressure in the T -p coordinates.The interpolation is performed using the method of the overlapping biquadratics that is included in the distribution of the OPAL EOS.The results for the internal energy and the density shown in the same range of the temperature and pressure as in Fig. 3 are visually nearly identical to the results of VMW.The electron density distribution shows a number of artefacts in the low-temperature region where ad hoc electron density is introduced (Rogers and Nayfonov, 2002), while it is similar to the VMW distribution elsewhere.Trampedach et al. (2006) examined the differences between the OPAL EOS and the MHD EOS.

Comparison
The relative differences between the mass density, the electron number density and the internal energy in the VMW and the OPAL EOS in the log T -log p plane are shown in Fig. 4. The positive values indicate that the VMW quantity is larger.The largest differences in the density and the internal energy occur in the area with a temperature below 10 6 K for pressure higher than 10 6 dyn cm −2 and densities higher than 10 −6 g cm −3 .This is the area where hydrogen is still neutral in the VMW EOS (cf.Fig. 2).The OPAL tables do not include the partial pressures of different particles needed for a direct comparison.However, the negative difference (the red area) in the density at constant pressure together with correspondingly higher internal energy and the electron density indicate that the VMW hydrogen ionisation fraction is higher than the one from OPAL.The opposite applies for the blue area.While the latter may be interpreted as the effect of the pressure ionisation (cf.Fig. 1), it is difficult to attribute the former to any particular effect without knowing the partial pressures consistent with the OPAL.Nevertheless, the area where the differences appear is not covered by the near-surface convection.At low pressure and a temperature below ≈ 2 × 104 K, the relative differences are below 1 % (except for the electron densities at low T , where the artefacts dominate).Figure 5 shows a blow-up of Fig. 4 in that region.Vertical stripes in the density plot and the horizontal stripes in the electron pressure are due to the interpolation errors (coming from the reverse interpolation of the OPAL table to the T -ρ domain) 4 .The vertical bands below 3500 K in the relative difference of the electron density correspond  to the area where the OPAL electron densities are unreliable.

The convection simulation
To study propagation with time of differences caused by EOS choice in a near-surface convection simulation, we initiate and run two simulations with all parameters identical except for the EOS lookup tables.For this experiment we use the version of the MURaM code described in Vögler (2003) with several minor adaptations needed to run the code with an alternative set of EOS lookup tables.The code requires two tables, one with the temperature and the density tabulated as functions of the internal energy and the pressure and another with the temperature and the pressure as functions of the internal energy and the density.The required tables are produced from the VMW and the OPAL EOS described above.Since we are primarily interested in the relative difference between the two runs and not in the physical realism of the result, we set up a simple HD convection run in 2-D with grey opacities.The physical size of the domain is 6.0 Mm in the horizontal direction and 1.4 Mm in the vertical, with 288 and 100 grid points respectively.The bottom boundary is open for the mass flow, the top boundary is closed and the vertical boundaries are set as periodic.The mass fluctuations through the bottom boundary are dynamically adjusted by the boundary pressure of the upflowing material p up (see Vögler, 2003, Eq. 3.56, p. 32).
The initial model is prepared by replicating the temperature and the density of a 1-D convectively stable model in two dimensions.The magnetic field is equal to 0 throughout the simulation runs.The internal energy is computed from these two quantities using both VMW and OPAL, so that we have two initial snapshots, one consistent with each EOS.The specific energy in the snapshots is then perturbed with a multiplicative factor of 1 ± 0.05, with identical randomisation over the simulation box.We start the runs with a damping factor for the radiative exchange term, gradually releasing it until the convection cells fully develop.Once the simulation is in the statistically steady state, we let it run for about 16 h of solar time while taking a snapshot every 20 s.

Results
In the first iteration step, the relative difference between the temperature in the two runs is within ±0.2%, which is in agreement with Fig. 5.As the flows develop and accelerate in the initial phase, the differences between the runs grow and the point-by-point comparison between their parameters becomes meaningless.Figure 6 shows the velocity for the two runs at constant height in the box (850 km above the bottom of the box) versus time.The flows in the two cases are indistinguishably similar over the first 13 min; after that initial phase they take different evolutionary paths toward the statistically steady solutions.
In MURaM the flux of the emergent radiation F out is controlled by the specific energy of the upflowing material at the bottom boundary so that its mean value matches the value of the real sun (see Vögler, 2003, Eq. 3.50, p. 30).The relative difference of the mean F out for the two runs is thus small, as expected: 0.3 %.The temperature profiles of the two runs averaged over time and over the horizontal direction are very similar.The difference in the geometrical height scale (left panel) and in the scale of the continuum optical depth at 500 nm (right) is shown in Fig. 7.The largest absolute difference between the two in the geometrical height scale is at 0.8 Mm above the bottom of the box and amounts to nearly 200 K (the relative difference is 2.3 %).This is the region where the radiative cooling begins to dominate producing a steep gradient of temperature with height.In the continuum optical depth scale at 500 nm, in the photosphere (τ 500 ≤ 1), the difference between the two temperature profiles is up to 10 K (less than 0.2 % relative to the local temperature).This difference causes the 0.3 % difference in the mean flux of the emergent intensity, which is negligible for most purposes.

Conclusions
The tables of thermodynamical quantities (the mass density, the specific internal energy and the electron density) computed using two equations of state -one ideal (VMW, after Vardya, 1965;Mihalas, 1967;Wittmann, 1974) and one nonideal (OPAL, Rogers et al., 1996;Rogers and Nayfonov, 2002) and both widely used in solar physics -are compared.The relative differences between the quantities in the region of the T -p plane usually considered in the realistic nearsurface convection simulations (approximately 1 Mm below the level where the mean continuum optical depth is equal to 1 and 0.5 Mm above it) are below 1 %, except for the electron density at the low temperatures, where the values of the OPAL tables are not reliable.Since we have access only to the tables precomputed with the OPAL code but not to the code itself, we cannot determine the origin of these differences.Possible causes are the differences in the atomic and molecular data, in the approach and the interpolation errors.
To check how this difference may affect the result of the hydrodynamical simulations, we simulated the 2-D solar convection with the two EOSs with an otherwise identical set-up.The two simulation runs separate in the parameter space as soon as the convective cells are developed.Nevertheless, they evolve separately to a statistically nearly identical steady state, with the temperature difference in the optical depth scale below 10 K throughout the photosphere.This experiment demonstrates that the two EOSs are interchangeable for simulations of near-surface solar convection.Moreover, the electron density in the OPAL tables, corrupted at low temperatures, may be recomputed for a large portion of the T -ρ grid of OPAL with an accuracy better than 1 % using the VMW ideal EOS (cf.Fig. 4).To estimate what the effect of the EOS choice in a more realistic set-up is (3-D, with magnetic field, non-grey radiative transfer including scattering, etc.), it would be necessary to repeat this experiment in such conditions.The change of the mean temperature between 2-D and 3-D solar convection simulations has been described by Asplund et al. (2000a) (see lower panel of their Fig.9).Without analysing the origin of the differences in their study, it may be noted that these differences are of the same order as the differences between the highly realistic 3-D HD simulations performed by different codes (Beeck et al., 2012).In any case, we do not expect significant changes in our results as long as the temperature and the density in a simulation remain within nearly the same range of the T -ρ plane as in our experiment.
This result is in agreement with the studies of Beeck et al. (2012) and Tanner et al. (2012), demonstrating robustness of the near-surface simulations.Regarding the deep convection, this result is consistent with the conclusions of Bahcall et al. (2004) that the depth of the convection zone does not depend significantly on the uncertainties in the EOS.It is also consistent with the recent study of Lord et al. (2014), who compared the effect of the OPAL EOS and a simple Saha-based EOS on the horizontal velocity spectrum of a solar deep convection simulation and found that the two produced nearly identical results, especially in the higher portion of the convection zone.

Figure 2 .
Figure 2. The ionisation fraction of hydrogen in the T -p grid computed with the VMW EOS.Those portions of the domain that are dominated by the neutral H, H 2 molecule and H + are indicated.

Figure 3 .
Figure 3.The mass density, the electron number density and the specific internal energy per mass computed on the temperature-pressure grid using the VMW EOS.

Figure 4 .
Figure 4.The relative difference between the VMW EOS and the OPAL EOS in the mass density, the electron number density and the internal energy ((VMW-OPAL)/VMW).The VMW EOS is computed on the T -p grid and the OPAL EOS computed on the T -ρ grid and then interpolated to the T -p grid.The two EOSs are computed for the same mass fractions.The contour lines are copied from Fig. 3 for comparison.

Figure 5 .
Figure 5. Blow-up of Fig. 4 in the low-pressure and low-temperature region relevant for the near-surface convection.The fragments of the contour lines from Fig. 4 are added.

Figure 6 .
Figure 6.Vertical velocity at fixed height (850 km above the bottom of the simulation box) over the first 24 min of the simulation with the OPAL EOS (left) and with the VMW EOS (right).The dashed line approximately marks the instance when the two simulation runs separate.

Figure 7 .
Figure7.The difference between the mean temperature of the simulation runs with OPAL and VMW in the geometrical height scale (left) and in the continuum optical depth scale at 500 nm (right).Note different scales in the two panels.