The following work was supported in part by NASA and NSF.
Equation of State of Hot, Dense Helium
Two first-principles simulation techniques, path integral Monte Carlo (PIMC) and density functional molecular dynamics (DFT-MD), are applied to study hot, dense helium in the density-temperature range of 0.387 ... 5.35 gcm-3 and 500 K ... 1.28*108 K. One coherent equation of state (EOS) is derived by combining DFT-MD data at lower temperatures with PIMC results at higher temperatures. Good agreement between both techniques is found in an intermediate temperature range. For the highest temperatures, the PIMC results converge to the Debye-Huckel limiting law. In order to derive the entropy, a thermodynamically consistent free energy fit is introduced that reproduces the internal energies and pressure derived from the first-principles simulations. The equation of state is presented in the form of a table as well as a fit.
This code derives the equation of hot, dense helium. Using a 2D spline interpolation as function of temperature and density, the code computes pressure, P, internal energy, E, Helmholtz free energy, F, and entropy, S.
This is all done free of charge but we recommended two citations:
B. Militzer,
"Path Integral Monte Carlo and Density Functional Molecular Dynamics Simulations of Hot, Dense Helium", Phys. Rev. B 79 (2009) 155105, cond-mat/08050317.
B. Militzer, "First Principles Calculations of Shock Compressed Fluid Helium",
Phys. Rev. Lett. 97 (2006) 175501.
(1) Download the code
(2) Installation
To install our helium EOS code, cd into your preferred directory, download and untar the latest version, and compile our simple code with your favorite Fortran 77 compile like GNU's gfortran for example:
my_machine:~/Downloads/helium_EOS> tar -xzf eoshe_11-04-22.tar.gz
my_machine:~/Downloads/helium_EOS> gfortran eoshe.f log55.f spline.f -o eoshe
In case the compilation fails, try switching to a different compiler or tinker with the compiler flag. (I do not use Fortran for my research but this code was written for my collaborators in 2009 who still used it at the time.)
(3) Files included
The EOS tables are stored in files
| eoshe.f | Main code that provides function call for P, E, F, and S with arguments T and rho. Prints a table as a test. |
| log55.f | 2D free energy interpolation. |
| spline.f | Simple 1D spline interpolation. |
| eoshe.out | Example output file for comparison. |
(4) Execution and tests
my_machine:~/Downloads/helium_EOS> eoshe > my_eoshe.out
my_machine:~/Downloads/helium_EOS> diff eoshe.out my_eoshe.out
If your output agrees with our example output, the diff should not print anything.
(4) Units and conventions
Here is what our example interpolation code prints:
rs= 2.400000 rhoHe(g/cc)= 0.387289 T(K)= 500.000 E(Ha/e)= -1.431865 P(GPa)= 1.406642 [F](Ha/e)= -1.437469 [S](Ha/kb/e)= 3.539252
rs= 2.400000 rhoHe(g/cc)= 0.387289 T(K)= 641.409 E(Ha/e)= -1.431474 P(GPa)= 1.594792 [F](Ha/e)= -1.439103 [S](Ha/kb/e)= 3.755843
rs= 2.400000 rhoHe(g/cc)= 0.387289 T(K)= 822.810 E(Ha/e)= -1.430913 P(GPa)= 1.835963 [F](Ha/e)= -1.441332 [S](Ha/kb/e)= 3.998785
...
rs= 0.800000 rhoHe(g/cc)= 10.456811 T(K)= 77782226.582 E(Ha/e)= 564.811592 P(GPa)= 5071000.761780 [F](Ha/e)= -4581.997299 [S](Ha/kb/e)= 20.894643
rs= 0.800000 rhoHe(g/cc)= 10.456811 T(K)= 99780383.856 E(Ha/e)= 715.528920 P(GPa)= 6510475.331458 [F](Ha/e)= -6058.023284 [S](Ha/kb/e)= 21.436239
rs= 0.800000 rhoHe(g/cc)= 10.456811 T(K)= 128000000.000 E(Ha/e)= 902.363130 P(GPa)= 8367649.190793 [F](Ha/e)= -7997.343159 [S](Ha/kb/e)= 21.955482
As we have done for our hydrogen-helium EOS table, this code uses atomic units:
Wigner-Seitz radius, rs, is given Bohr radii = 5.29177208607*10-11 meters.
4/3*pi*rs3=V/Ne where Ne is total number of electrons (bound or free) in volume V.
The number density of electrons is n=Ne / V. Every helium atom contributes two electrons regardless of its ionization state.
The Helmholtz free energy, F, and the internal energy, E, are given in Hartree/electron = 4.35974380267*10-18 Joules
The entory, S, is given in kb/electron where kb is Boltzmann's constant.
(5) Limitations
We provide thermodynamic quantities on a rectangular grid in temperature and density.
The temperature range is very large from 500 to 1.28*108 K. The density range spans from 0.387289 to 5.35 g/cm3.
So conditions in the deep interior of a giant planet are well covered but the low-density outer layers are not included and one needs to combine our results with another EOS table.
If you have feedback, comments, or a specific scientific problem you would like to solve please send email to militzer[at]berkeley[dot]edu.
Best wishes and happy EOS computation!
Burkhard Militzer
November, 4, 2022.