© ESO, 2013

1. Introduction

Coulomb plasmas, i.e., the fully-ionized plasmas whose thermodynamics is strongly affected by electrostatic interactions, are encountered in many physical and astrophysical situations (e.g., Fortov 2009). Full ionization is reached either at high temperatures T and low densities ρ (thermal ionization) or at high densities ρ (pressure ionization). The latter case is typical of the interior conditions of low-mass stars, brown dwarfs, or giant planets (Chabrier & Baraffe 2000) as well as the interior and envelope conditions of white dwarfs and neutron stars. Coulomb interactions are crucial for the equation of state (EOS) under such conditions. In the interior or the envelope of compact objects such as white dwarfs and neutron stars, the electrons can be weakly or strongly degenerate, the plasma can be in the liquid or solid state, the electrons can have various degrees of degeneracy and relativism, and the quantum effects on ion motion can be substantial. Therefore, a wide-range EOS is needed for calculations of the structure and evolution of such stars.

In a previous work (Potekhin & Chabrier 2000, 2010, hereafter Papers I and II, respectively) we proposed a set of analytical expressions for the calculations of the EOSs of the Coulomb plasmas without magnetic fields and presented a code for thermodynamic functions based on the first, second, and mixed derivatives of the analytical Helmholtz free energy F with respect to density ρ and temperature T. This code has been employed in astrophysical modeling and adapted for the use in the Modules for Experiments in Stellar Astrophysics (mesa; Paxton et al. 2011).

The Bohr – van Leeuwen theorem states that an EOS of charged pointlike classical particles is not affected by a magnetic field (van Leeuwen 1921). However, a magnetic field can affect thermodynamic functions through intrinsic magnetic moments of particles and by the quantization of the motion of charged particles in Landau orbitals (Landau 1930; Landau & Lifshitz 1977). These effects can be important, for example, in magnetic white dwarfs whose magnetic fields B can reach 10⁷−10⁹ G (e.g., Wickramasinghe & Ferrario 2000, and references therein) and neutron stars with typical B ~ 10⁸−10¹⁴ G (e.g., Haensel et al. 2007, and references therein).

In this paper, we systematically consider analytical expressions for thermodynamic functions of magnetized Coulomb plasmas, discuss their validity range, and introduce some practical modifications. We take account of analytical and numerical results, currently available for various contributions to the Helmholtz free energy in quantizing magnetic fields. Taking advantage of recently published numerical results (Baiko 2009), we construct an analytical description of the thermodynamic functions of harmonic Coulomb crystals in quantizing magnetic fields.

In Sect. 2 we give definitions and simple estimates for the plasma parameters that determine different thermodynamic regimes. In Sect. 3 we outline the EOS of a nonmagnetized Coulomb plasma as the reference case. In Sect. 4 we consider the Boltzmann and Fermi gases in quantizing magnetic fields, present a general analytical description of their EOS, and simplify them for several limiting cases. In Sect. 5 we review nonideal contributions to the EOS of a Coulomb liquid in a strong magnetic field. In Sect. 6 we derive an analytical approximation for the EOS of a strongly magnetized Coulomb crystal. In Sect. 7 we present and discuss examples of thermodynamic functions for conditions typical of white dwarfs and neutron-star envelopes. The summary is given in Sect. 8. In the Appendices we give the explicit expressions for the thermodynamic functions used in Sect. 3.

2. Basic definitions

2.1. General parameters

Let n_e and n_i be the electron and ion number densities, A and Z the ion mass and charge numbers, respectively. In this paper we consider the neutral plasmas (therefore n_e = Zn_i) that contain a single type of ion and include neither positrons (they can be described using the same analytical functions as for the electrons; see, e.g., Blinnikov et al. 1996; Timmes & Arnett 1999), nor free neutrons (see Haensel et al. 2007 for a review).

The state of a free-electron gas is determined by the electron number density n_e and temperature T. Instead of n_e it is convenient to introduce the dimensionless density parameter r_s = a_e/a₀, where a₀ is the Bohr radius and ${\mathit{a}}_{\mathrm{e}}\mathrm{=}\mathrm{(}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{n}}_{\mathrm{e}}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}$. The parameter r_s can be quickly evaluated from the relations ${\mathit{r}}_{\mathrm{s}}\mathrm{=}\mathrm{1.1723}\hspace{0.17em}{\mathit{n}}_{\mathrm{24}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{=}\mathrm{(}{\mathit{\rho }}_{\mathrm{0}}\mathit{/}\mathit{\rho }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$ where n₂₄ ≡ n_e/10²⁴   cm^-3 and ρ₀ = 2.6752   (A/Z) g cm^-3. The analogous density parameter for the ions is R_S = a_im_i(Ze)²/ħ² = 1822.89AZ^7/3   r_s, where m_i is the ion mass and ${\mathit{a}}_{\mathrm{i}}\mathrm{\equiv }\mathrm{(}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }{\mathit{n}}_{\mathrm{i}}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}$ is the ion sphere radius.

At stellar densities it is convenient to use, instead of r_s, the nonmagnetic relativity parameter ${\mathit{x}}_{\mathrm{r}}\mathrm{=}{\mathit{p}}_{\mathrm{F}}\mathit{/}{\mathit{m}}_{\mathrm{e}}\mathit{c}\mathrm{=}\mathrm{1.00884}\hspace{0.17em}{\left({\mathit{\rho }}_{\mathrm{6}}\mathit{Z}\mathit{/}\mathit{A}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{=}\mathrm{0.014005}\hspace{0.17em}{\mathit{r}}_{\mathrm{s}}^{-1}\mathit{,}$(1)where p_F = ħ   (3π²n_e)^1/3 is the electron Fermi momentum in the absence of a magnetic field, and ρ₆ ≡ ρ/10⁶ g cm^-3. The Fermi energy (without the rest energy) for the electron gas is ${\mathit{ϵ}}_{\mathrm{F}}\mathrm{=}\mathit{c}\hspace{0.17em}\sqrt{\mathrm{(}{\mathit{m}}_{\mathrm{e}}\mathit{c}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{p}}_{\mathrm{F}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{-}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}\mathit{,}$ and the Fermi temperature T_F ≡ ϵ_F/k_B = T_r   (γ_r − 1), where T_r ≡ m_ec²/k_B = 5.93 × 10⁹   K, ${\mathit{\gamma }}_{\mathrm{r}}\mathrm{\equiv }\sqrt{\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathrm{r}}^{\mathrm{2}}}$, and k_B is the Boltzmann constant. A useful measure of electron degeneracy is θ = T/T_F. In the nonrelativistic limit (x_r ≪ 1), ${\mathit{T}}_{\mathrm{F}}\mathrm{\approx }\mathrm{1.163}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}\hspace{0.17em}{\mathit{r}}_{\mathrm{s}}^{-2}$ K, and $\mathit{\theta }\mathrm{=}\mathrm{0.543}\hspace{0.17em}{\mathit{r}}_{\mathrm{s}}\mathit{/}{\mathrm{\Gamma }}_{\mathrm{e}}\mathit{,}$(2)where ${\mathrm{\Gamma }}_{\mathrm{e}}\mathrm{\equiv }\frac{{\mathit{e}}^{\mathrm{2}}}{{\mathit{a}}_{\mathrm{e}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\mathrm{\approx }\frac{\mathrm{22.747}}{{\mathit{T}}_{\mathrm{6}}}{\left({\mathit{\rho }}_{\mathrm{6}}\frac{\mathit{Z}}{\mathit{A}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}$(3)In the opposite ultrarelativistic limit (x_r ≫ 1), θ ≈ (263   Γ_e)^-1. The strength of the Coulomb interaction of nonrelativistic ions is characterized by the Coulomb coupling parameter $\mathrm{\Gamma }\mathrm{=}\frac{\mathrm{(}\mathit{Ze}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{a}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{e}}{\mathit{Z}}^{\mathrm{5}\mathit{/}\mathrm{3}}\mathit{,}$(4)where T₆ ≡ T/10⁶ K.

Thermal de Broglie wavelengths of free ions and electrons are usually defined as ${\mathit{\lambda }}_{\mathrm{i}}\mathrm{=}{\left(\frac{\mathrm{2}\mathit{\pi }{ħ}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,} {\mathit{\lambda }}_{\mathrm{e}}\mathrm{=}{\left(\frac{\mathrm{2}\mathit{\pi }{ħ}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{e}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}\mathit{,}$(5)although in some publications these definitions differ by a numerical factor. The quantum effects on ion motion are important either at λ_i ≳ a_i or at T ≪ T_p, where T_p ≡ ħω_p/k_B is the ion plasma temperature and ω_p = (4πe²   n_iZ²/m_i(^1/2 is the ion plasma frequency. Since ${\mathit{\lambda }}_{\mathrm{i}}\mathit{/}{\mathit{a}}_{\mathrm{i}}\mathrm{=}\mathrm{(}{\mathit{T}}_{\mathrm{p}}\mathit{/}\mathit{T}\mathrm{)}\sqrt{\mathrm{2}\mathit{/}\mathrm{\Gamma }}$, the importance of the quantum effects in strongly coupled plasmas (i.e., at Γ ≫ 1) is determined by parameter $\mathit{\eta }\mathrm{\equiv }{\mathit{T}}_{\mathrm{p}}\mathit{/}\mathit{T}\mathrm{=}\mathrm{\Gamma }\sqrt{\mathrm{3}\mathit{/}{\mathit{R}}_{\mathrm{S}}}\mathrm{\approx }\mathrm{7.832}\hspace{0.17em}\mathrm{(}\mathit{Z}\mathit{/}\mathit{A}\mathrm{)}\sqrt{{\mathit{\rho }}_{\mathrm{6}}}\mathit{/}{\mathit{T}}_{\mathrm{6}}\mathit{.}$(6)

2.2. Magnetic-field parameters

In the nonrelativistic theory (Landau & Lifshitz 1977), the energy of an electron in magnetic field B equals $\mathit{n}ħ{\mathit{\omega }}_{\mathrm{c}}\mathrm{+}{\mathit{m}}_{\mathrm{e}}{\mathit{p}}_{\mathit{z}}^{\mathrm{2}}\mathit{/}\mathrm{2}$, where p_z is the momentum component along B, ω_c = eB/m_ec is the electron cyclotron frequency, $\mathit{n}\mathrm{=}{\mathit{n}}_{\mathrm{L}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\sigma }$ characterizes a Landau level, σ =  ±1 determines the spin projection on the field, and n_L is the non-negative integer Landau number related to the quantization of the kinetic motion transverse to the field. In the relativistic theory (Johnson & Lippmann 1949; Berestetskiĭ et al. 1982), the kinetic energy ϵ of an electron at the Landau level n and its longitudinal momentum p_z are inter-related as $\begin{array}{ccc}& & \mathit{ϵ}\mathrm{=}{\mathit{ϵ}}_{\mathit{n}}\mathrm{\left(}{\mathit{p}}_{\mathit{z}}\mathrm{\right)}\mathrm{=}\mathit{c}\hspace{0.17em}{\mathrm{(}{\mathit{m}}_{\mathrm{e}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}\mathrm{+}\mathrm{2}ħ{\mathit{\omega }}_{\mathrm{c}}{\mathit{m}}_{\mathrm{e}}\mathit{n}\mathrm{+}{{\mathit{p}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}\mathit{,}\\ & & \mathrm{\left|}{\mathit{p}}_{\mathit{z}}\mathrm{\right|}\mathrm{=}{\mathit{p}}_{\mathit{n}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}{\mathrm{[}\mathrm{(}{\mathit{m}}_{\mathrm{e}}\mathit{c}\mathrm{+}\mathit{ϵ}\mathit{/}\mathit{c}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{(}{\mathit{m}}_{\mathrm{e}}\mathit{c}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{m}}_{\mathrm{e}}ħ{\mathit{\omega }}_{\mathrm{c}}{\mathit{n}}^{\mathrm{]}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}\end{array}$The levels are double-degenerate with respect to σ. Their splitting due to the anomalous magnetic moment of the electron is  ≈ (m_ec²α_f/2π)   b at b ≪ 1 and  ~(m_ec²α_f/2π)    [lnb − 1.584] ² at b ≫ 1 (see Schwinger 1988; Suh & Mathews 2001), which is always much smaller than ħω_c and is negligible in the compact stars.

Convenient dimensionless parameters that characterize the magnetic field in a plasma are the ratios of the electron cyclotron energy ħω_c to the Hartree unit of energy, to the electron rest energy, and to k_BT: ${\mathit{\gamma }}_{\mathrm{m}}\mathrm{=}{ħ}^{\mathrm{3}}\mathit{B}\mathit{/}{\mathit{m}}_{\mathrm{e}}^{\mathrm{2}}\mathit{c}{\mathit{e}}^{\mathrm{3}}\mathrm{=}\mathit{B}\mathit{/}{\mathit{B}}_{\mathrm{0}}\mathit{,}$(9)where B₀ = 2.3505 × 10⁹   G, $\mathit{b}\mathrm{=}\frac{ħ{\mathit{\omega }}_{\mathrm{c}}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\mathrm{=}{\mathit{\alpha }}_{\mathrm{f}}^{\mathrm{2}}\hspace{0.17em}{\mathit{\gamma }}_{\mathrm{m}}\mathrm{=}\frac{\mathit{B}}{\mathrm{4.414}\mathrm{\times }{\mathrm{10}}^{\mathrm{13}} \mathrm{G}}\hspace{0.17em}\mathit{,}$(10)where α_f = e²/ħc is the fine-structure constant, and $\mathit{\zeta }\mathrm{=}ħ{\mathit{\omega }}_{\mathrm{c}}\mathit{/}{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathrm{\approx }\mathrm{134.34}{\mathit{B}}_{\mathrm{12}}\mathit{/}{\mathit{T}}_{\mathrm{6}}\mathit{,}$(11)where B₁₂ ≡ B/10¹² G. The magnetic length ${\mathit{a}}_{\mathrm{m}}\mathrm{=}\mathrm{(}ħ\mathit{c}\mathit{/}\mathit{eB}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}{\mathit{a}}_{\mathrm{0}}\mathit{/}\sqrt{{\mathit{\gamma }}_{\mathrm{m}}}$ gives a characteristic transverse scale of the electron wave function.

For the ions, the cyclotron energy is ħω_ci = Z   (m_e/m_i)   ħω_c, and the parameter analogous to ζ is ${\mathit{\zeta }}_{\mathrm{i}}\mathrm{=}ħ{\mathit{\omega }}_{\mathrm{ci}}\mathit{/}{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathrm{\approx }\mathrm{0.0737}\hspace{0.17em}\mathrm{(}\mathit{Z}\mathit{/}\mathit{A}\mathrm{)}{\mathit{B}}_{\mathrm{12}}\mathit{/}{\mathit{T}}_{\mathrm{6}}\mathit{.}$(12)Another important parameter is the ratio of the ion cyclotron frequency to the plasma frequency, $\mathit{\beta }\mathrm{=}{\mathit{\omega }}_{\mathrm{ci}}\mathit{/}{\mathit{\omega }}_{\mathrm{p}}\mathrm{=}{\mathit{\zeta }}_{\mathrm{i}}\mathit{/}\mathit{\eta }\mathrm{\approx }\mathrm{0.0094}\hspace{0.17em}{\mathit{B}}_{\mathrm{12}}\mathit{/}\sqrt{{\mathit{\rho }}_{\mathrm{6}}}\mathit{.}$(13)

2.3. Free energy and thermodynamic functions

The Helmholtz free energy F of a plasma can be conveniently written as $\mathit{F}\mathrm{=}{\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}\mathrm{+}{\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{+}{\mathit{F}}_{\mathrm{ee}}\mathrm{+}{\mathit{F}}_{\mathrm{ii}}\mathrm{+}{\mathit{F}}_{\mathrm{ie}}\mathit{,}$(14)where ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}$ and ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}$ denote the ideal free energy of the ions and the electrons, and the last three terms represent an excess free energy arising from the electron-electron, ion-ion, and ion-electron interactions, respectively. In the nonideal plasmas, correlations between any plasma particles depend on all interactions, therefore the separation in Eq. (14) is just a question of convenience.

An important reference case is the model of one-component plasma (OCP). In this model, the electrons are replaced by a rigid (nonpolarizable) background of the uniform charge distribution. It is convenient to define F_ii as the difference between F and ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}$ in the OCP model. Still stronger simplification is the ion-sphere model, in which the interaction energy in the OCP is evaluated as the electrostatic energy of a positive ion in the negatively charged sphere of radius a_i (Salpeter 1961). The electron exchange-correlation term is defined as ${\mathit{F}}_{\mathrm{ee}}\mathrm{=}\mathit{F}\mathrm{-}{\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}$ in the model of an electron gas without consideration of the ions, which are replaced by an uniform positive background to ensure the global charge neutrality. The ion-electron (electron polarization) contribution F_ie, then, is the difference between F and the other terms, when interactions between all types of particles are taken into account.

The pressure P, the internal energy U, and the entropy S of an ensemble of particles in volume V can be found from the thermodynamic relations P = −(∂F/∂V)_T,   S = −(∂F/∂T)_V, and U = F + TS. The second-order thermodynamic functions are derived by differentiating these first-order functions. The decomposition (Eq. (14)) induces analogous decompositions of P, U, S, the heat capacity C_V = (∂S/∂lnT)_V, and the logarithmic derivatives χ_T = (∂lnP/∂lnT)_V and χ_ρ = −(∂lnP/∂lnV)_T. Other second-order functions can be expressed through these functions by Maxwell relations (e.g., Landau & Lifshitz 1980).

3. EOS of nonmagnetized Coulomb plasmas

3.1. Ideal part of the free energy

The free energy of a gas of N_i = n_iV nonrelativistic classical ions is ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}\mathrm{=}{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathrm{\left[}\ln \mathrm{\right(}{\mathit{n}}_{\mathrm{i}}{\mathit{\lambda }}_{\mathrm{i}}^{\mathrm{3}}\mathit{/}{\mathit{M}}_{\mathrm{spin}}\mathrm{)}\mathrm{-}{\mathrm{1}}^{\mathrm{]}}\mathit{,}$(15)where M_spin is the spin multiplicity. Accordingly, ${\mathit{U}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\hspace{0.17em}{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}$, ${\mathit{P}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}\mathrm{=}{\mathit{n}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}$, ${\mathit{C}}_{\mathit{V,}\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}\hspace{0.17em}{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}$, and ${\mathit{\chi }}_{\mathit{T,}\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}\mathrm{=}{\mathit{\chi }}_{\mathit{\rho ,}\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}\mathrm{=}\mathrm{1}$. In the OCP, Eq. (15) can be written in terms of the dimensionless plasma parameters (Sect. 2) as $\frac{{\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}}{{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\mathrm{=}\mathrm{3}\ln \mathit{\eta }\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\ln \mathrm{\Gamma }\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\ln \frac{\mathrm{6}}{\mathit{\pi }}\mathrm{-}\ln {\mathit{M}}_{\mathrm{spin}}\mathrm{-}\mathrm{1.}$(16)The free energy of the electron gas is given by ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{=}{\mathit{\mu }}_{\mathrm{e}}{\mathit{N}}_{\mathrm{e}}\mathrm{-}{\mathit{P}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\hspace{0.17em}\mathit{V,}$(17)where N_e = n_eV is the number of electrons and μ_e is the electron chemical potential without the rest energy m_ec². The pressure and the number density are functions of μ_e and T: $\begin{array}{ccc}{\mathit{P}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}& \mathrm{=}& \frac{\mathrm{8}}{\mathrm{3}\sqrt{\mathit{\pi }}}\hspace{0.17em}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{{\mathit{\lambda }}_{\mathrm{e}}^{\mathrm{3}}}\left[{\mathit{I}}_{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{e}}\mathit{,\tau }\mathrm{\right)}\mathrm{+}\frac{\mathit{\tau }}{\mathrm{2}}{\mathit{I}}_{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{e}}\mathit{,\tau }\mathrm{\right)}\right]\mathit{,}\\ {\mathit{n}}_{\mathrm{e}}& \mathrm{=}& \frac{\mathrm{4}}{\sqrt{\mathit{\pi }}\hspace{0.17em}{\mathit{\lambda }}_{\mathrm{e}}^{\mathrm{3}}}\left[{\mathit{I}}_{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{e}}\mathit{,\tau }\mathrm{\right)}\mathrm{+}\mathit{\tau }{\mathit{I}}_{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{e}}\mathit{,\tau }\mathrm{\right)}\right]\mathit{,}\end{array}$where χ_e ≡ μ_e/k_BT, τ ≡ T/T_r, and ${\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{e}}\mathit{,\tau }\mathrm{\right)}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{{\mathit{x}}^{\mathit{\nu }}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\tau x}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\exp \mathrm{(}\mathit{x}\mathrm{-}{\mathit{\chi }}_{\mathrm{e}}\mathrm{)}\mathrm{+}\mathrm{1}}\hspace{0.17em}\mathrm{d}\mathit{x}$(20)is the Fermi-Dirac integral. The internal energy is ${\mathit{U}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{4}}{\sqrt{\mathit{\pi }}}\hspace{0.17em}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{TV}}{{\mathit{\lambda }}_{\mathrm{e}}^{\mathrm{3}}}\left[{\mathit{I}}_{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{e}}\mathit{,\tau }\mathrm{\right)}\mathrm{+}\mathit{\tau }\hspace{0.17em}{\mathit{I}}_{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{e}}\mathit{,\tau }\mathrm{\right)}\right]\mathit{.}$(21)Since we use V and T as independent variables, we need to find μ_e(V,T). This can be done either by inverting Eq. (19) numerically, or from the analytical approximation given in Chabrier & Potekhin (1998). Then the second-order thermodynamic functions are obtained using relations of the type $\begin{array}{ccc}{\left(\frac{\mathit{\partial f}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{e}}\mathit{,T}\mathrm{\right)}}{\mathit{\partial T}}\right)}_{\mathit{V}}& \mathrm{=}& {\left(\frac{\mathit{\partial f}}{\mathit{\partial T}}\right)}_{{\mathit{\chi }}_{\mathrm{e}}}\mathrm{-}{\left(\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{\chi }}_{\mathrm{e}}}\right)}_{\mathit{T}}\hspace{0.17em}\frac{\mathrm{(}\mathit{\partial }{\mathit{n}}_{\mathrm{e}}\mathit{/}\mathit{\partial T}{\mathrm{)}}_{{\mathit{\chi }}_{\mathrm{e}}}}{\mathrm{(}\mathit{\partial }{\mathit{n}}_{\mathrm{e}}\mathit{/}\mathit{\partial }{\mathit{\chi }}_{\mathrm{e}}{\mathrm{)}}_{\mathit{T}}}\mathit{,}\\ {\left(\frac{\mathit{\partial f}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{e}}\mathit{,T}\mathrm{\right)}}{\mathit{\partial V}}\right)}_{\mathit{T}}& \mathrm{=}& \mathrm{-}\frac{{\mathit{n}}_{\mathrm{e}}}{\mathit{V}}{\left(\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{n}}_{\mathrm{e}}}\right)}_{\mathit{T}}\mathrm{=}\mathrm{-}\frac{{\mathit{n}}_{\mathrm{e}}}{\mathit{V}}\hspace{0.17em}\frac{\mathrm{(}\mathit{\partial f}\mathit{/}\mathit{\partial }{\mathit{\chi }}_{\mathrm{e}}{\mathrm{)}}_{\mathit{T}}}{\mathrm{(}\mathit{\partial }{\mathit{n}}_{\mathrm{e}}\mathit{/}\mathit{\partial }{\mathit{\chi }}_{\mathrm{e}}{\mathrm{)}}_{\mathit{T}}}\mathrm{·}\end{array}$We use analytical approximations for I_ν(χ_e,τ) based on the fits of Blinnikov et al. (1996) and accurate typically to a few parts in 10⁴, with maximum error  <0.2% at τ ≤ 100 (Chabrier & Potekhin 1998). These approximations are given by different expressions in three ranges of χ_e: below, within, and above the interval 0.6 ≤ χ_e < 14. In particular, at large χ_e the Sommerfeld expansion (e.g., Chandrasekhar 1957; Girifalco 1973) yields¹${\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{e}}\mathit{,\tau }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}}\hspace{0.17em}{\mathit{\tau }}^{\mathit{\nu }\mathrm{+}\mathrm{1}}}\left({\mathrm{ℐ}}_{\mathit{\nu }}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{\mu }\end{array}\mathrm{\right)}\mathrm{+}\frac{{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{6}}\hspace{0.17em}{\mathit{\tau }}^{\mathrm{2}}{\mathrm{ℐ}}_{\mathit{\nu }}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{\mu }\end{array}\mathrm{\right)}\mathrm{+}\mathit{...}\right)\mathit{,}$(24)where $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\mu }\end{array}\mathrm{\equiv }{\mathit{\chi }}_{\mathrm{e}}\mathit{\tau }\mathrm{=}{\mathit{\mu }}_{\mathrm{e}}\mathit{/}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}$ is the electron chemical potential (without the rest energy) in the relativistic units, $\begin{array}{ccc}& & \\ & & {\mathrm{ℐ}}_{\mathrm{3}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{\mu }\end{array}\mathrm{\right)}\mathrm{=}\stackrel{}{{{\mathit{x}}^{\mathrm{3}}}_{\mathrm{˜}}}\mathit{/}\mathrm{3}\mathrm{-}{\mathrm{ℐ}}_{\mathrm{1}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{\mu }\end{array}\mathrm{\right)}\mathit{,}\\ & & \\ & & {\mathrm{ℐ}}_{\mathit{\nu }}^{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ \mathit{\mu }\end{array}\mathrm{\right)}\mathrm{=}\mathrm{d}{\mathrm{ℐ}}_{\mathit{\nu }}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathit{/}\mathrm{d}\begin{array}{}\mathrm{˜}\\ \mathit{\mu }\end{array}\mathit{.}\end{array}$Here, we have introduced notations $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{x}\end{array}\mathrm{\equiv }\sqrt{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\mu }\end{array}\mathrm{(}\mathrm{2}\mathrm{+}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\mu }\end{array}\mathrm{)}}$ and $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\gamma }\end{array}\mathrm{\equiv }\sqrt{\mathrm{1}\mathrm{+}\stackrel{}{{{\mathit{x}}^{\mathrm{2}}}_{\mathrm{˜}}}}\mathrm{=}\mathrm{1}\mathrm{+}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\mu }\end{array}$. At strong degeneracy, $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\mu }\end{array}\mathrm{\approx }{\mathit{ϵ}}_{\mathrm{F}}\mathit{/}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}\mathrm{-}\mathrm{1}$, $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{x}\end{array}\mathrm{\approx }{\mathit{x}}_{\mathrm{r}}$, and $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\gamma }\end{array}\mathrm{\approx }{\mathit{\gamma }}_{\mathrm{r}}$. In Paper II we also described an alternative expansion in powers of τ, which allows one to avoid numerical cancellations of close terms at small $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{\mu }\end{array}$ (we switch to this alternative expansion at $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\mu }\end{array}\mathit{<}\mathrm{0.1}$).

The discontinuities of the Blinnikov et al. (1996) approximations for I_ν(χ_e,τ) at χ_e = 0.6 and χ_e = 14 are typically a few parts in 10⁴ at τ ≲ 10², but they may reach  ≈ 1% for the second derivatives. This accuracy is sufficient for many applications. Nevertheless, the jumps may produce problems, e.g., when higher derivatives are evaluated numerically in a stellar evolution code. In our calculations of white-dwarf evolution (to be published elsewhere), we smoothly interpolate between the two analytical approximations for the adjacent intervals near the boundary at the cost of a slight violation of the thermodynamic consistency in the interpolation regions (this version of the EOS code is now also available at our web site).

If a higher accuracy is needed, one can numerically calculate tables of I_ν(χ_e,τ) (e.g., Timmes & Arnett 1999) and interpolate in them with an algorithm that preserves thermodynamic consistency (Timmes & Swesty 2000) and is available at mesa (Paxton et al. 2011).

3.2. Nonideal contributions

3.2.1. Electron and ion liquids

The contribution to the free energy due to the electron-electron interactions has been studied by many authors. For the reasons explained in Paper II, we adopt the fit to F_ee derived by Ichimaru et al. (1987A.1) (see Appendix A).

The ion-ion interactions are described using the OCP model. In the liquid regime, the numerical results obtained for the OCP of nonrelativistic pointlike charged particles in different intervals of the Coulomb coupling parameter from Γ = 0 to Γ ~ 200 by different numerical and analytical methods are reproduced by a simple expression given in Appendix B.1. The accurate fit for classical OCP is supplemented by the Wigner-Kirkwood correction, which extends the applicability range of our approximations to lower temperatures T ~ T_p. In spite of the significant progress in numerical ab initio modeling of quantum ion liquids, available results do not currently allow us to establish an analytical extension to still lower temperatures T ≪ T_p (see Chabrier et al. 2002, for references and discussion).

3.2.2. Coulomb crystal

At T < T_m, where T_m is the melting temperature, the ions in thermodynamic equilibrium are arranged in the body-centered cubic (bcc) Coulomb lattice. In the harmonic approximation (e.g., Kittel 1963), the free energy of the lattice is ${\mathit{F}}_{\mathrm{lat}}\mathrm{=}{\mathit{U}}_{\mathrm{0}}\mathrm{+}{\mathit{U}}_{\mathrm{q}}\mathrm{+}{\mathit{F}}_{\mathrm{th}}\mathit{,}$(29)where U₀ = N_iC₀(Ze)²/a_i is the classical static-lattice energy, C₀ ≈ −0.9 is the Madelung constant, ${\mathit{U}}_{\mathrm{q}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{N}}_{\mathrm{i}}ħ{\mathit{\omega }}_{\mathrm{p}}{\mathit{u}}_{\mathrm{1}}$(30)accounts for zero-point quantum vibrations, u₁ =  ⟨ ω_kα ⟩ _ph/ω_p ≈ 0.5 is the reduced first moment of phonon frequencies, ${\mathit{F}}_{\mathrm{th}}\mathrm{=}\mathrm{3}{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}{⟨\ln \mathrm{[}\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}ħ{\mathit{\omega }}_{{k}\mathit{\alpha }}\mathit{/}{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathrm{\left)}\mathrm{\right]}⟩}_{\mathrm{ph}}$(31)is the thermal contribution, ω_kα are phonon frequencies, and  ⟨ ... ⟩ _ph denotes the averaging over phonon polarizations α and wave vectors k in the first Brillouin zone. Here we do not separate the classical-gas free energy, therefore F_lat replaces ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}\mathrm{+}{\mathit{F}}_{\mathrm{ii}}^{}$ in Eq. (14).

Beyond the harmonic-lattice approximation, the total reduced free energy f_lat ≡ F_lat/N_ik_BT can be written as ${\mathit{f}}_{\mathrm{lat}}\mathrm{=}{\mathit{C}}_{\mathrm{0}}\mathrm{\Gamma }\mathrm{+}\mathrm{1.5}{\mathit{u}}_{\mathrm{1}}\mathit{\eta }\mathrm{+}{\mathit{f}}_{\mathrm{th}}\mathrm{+}{\mathit{f}}_{\mathrm{ah}}\mathit{.}$(32)Here, the first three terms correspond to the three terms in Eq. (29), and f_ah is the anharmonic correction. The most accurate values of the constants C₀ and u₁ were calculated by Baiko (2000) (see Appendix B.2). For f_th = F_th/N_ik_BT, we use the highly precise fit of Baiko et al. (2001) (Appendix B.2). In the classical limit η ≪ 1, it reduces to f_th ≃ 3lnη − 2.49389−1.5u₁η + η²/24, where the term with u₁ cancels that in Eq. (32) and the last term represents the Wigner-Kirkwood quantum correction ${\mathit{f}}_{\mathrm{ii}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ (Eq. (B.3)), which is the same in the liquid and solid phases (Pollock & Hansen 1973). In the opposite limit T ≪ T_p, we have f_th ≃ −209.3323   η^-3 (here the constant is given for the bcc crystal; for other lattice types, see Baiko et al. 2001).

Anharmonic corrections for Coulomb lattices were studied by many authors in the limits η → 0 and η → ∞, but only a few numerical results of low precision are available at finite η values (see Paper II for references and discussion). In Paper II we constructed an analytical interpolation between these limits, which is applicable at arbitrary η and is consistent with the available numerical results within accuracy of the latter ones (Appendix B.2). It should be replaced by a more accurate function in the future when accurate finite-temperature anharmonic quantum corrections become available.

3.2.3. Electron polarization

Electron polarization in Coulomb liquids was studied by perturbation (Galam & Hansen 1976; Yakovlev & Shalybkov 1989) and hypernetted-chain (HNC) techniques (Chabrier & Ashcroft 1990; Chabrier & Potekhin 1998; Paper I). The results have been reproduced by an analytical expression (Appendix C.1), which exactly recovers the Debye-Hückel limit for the weakly coupled (Γ ≪ 1) electron-ion plasmas and the Thomas-Fermi limit for the strongly coupled (Γ ≫ 1) Coulomb liquids at Z ≫ 1.

For classical ions, the simplest screening model consists in replacing the Coulomb potential by the Yukawa potential. Molecular-dynamics and path-integral Monte Carlo simulations of classical liquid and solid Yukawa systems were performed in several works (e.g., Hamaguchi et al. 1997; Militzer & Graham 2006). However, the Yukawa interaction reflects only the small-wavenumber asymptote of the electron dielectric function (Jancovici 1962; Galam & Hansen 1976). The first-order perturbation approximation for the dynamical matrix of a classical Coulomb solid with the polarization corrections was developed by Pollock & Hansen (1973). The phonon spectrum in such a quantum crystal has been calculated only in the harmonic approximation (Baiko 2002), which has a restricted applicability to this problem (for example, it is obviously incapable of reproducing the polarization contribution to the heat capacity in the classical limit η → 0, where it gives C_V = 3N_ik_B independent of the polarization).

In Paper I we calculated F_ie using the semiclassical perturbation theory of Galam & Hansen (1976) with a model structure factor, and fit the results by an analytical function of x_r and η. In Paper II we improved the η-dependence of this function to completely eliminate the screening contribution in the strong quantum limit η ≪ 1, because the employed model of the structure factor failed at η ≲ 1. The latter approximation is reproduced in Appendix C.2. It can be improved in the future, when the polarization corrections for the quantum Coulomb crystal at η ≲ 1 have been accurately evaluated.

3.2.4. Ion mixtures

In Sects. 3.2.1–3.2.3 we have considered plasmas containing identical ions. In the case where several ion types are present in a strongly coupled Coulomb plasma, a common approximation is the linear mixing rule (LMR), ${\mathit{f}}_{\mathrm{ex}}^{\mathrm{LM}}\mathrm{\approx }\sum _{\mathit{j}}{\mathit{x}}_{\mathit{j}}{\mathit{f}}_{\mathrm{ex}}\mathrm{(}{\mathrm{\Gamma }}_{\mathit{j}}\mathit{,}{\mathit{x}}_{\mathit{j}}\mathrm{=}\mathrm{1}\mathrm{)}\mathit{,}$(33)where x_j are the number fractions of ions with charge numbers Z_j and ${\mathrm{\Gamma }}_{\mathit{j}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{e}}\hspace{0.17em}{\mathit{Z}}_{\mathit{j}}^{\mathrm{5}\mathit{/}\mathrm{3}}$. In Eq. (33), f_ex ≡ F_ex/N_ik_BT is the reduced nonideal part of the free energy, F_ex is the excess free energy, which is equal to F_ii in the case of the rigid charge-neutralizing electron background and to F_ii + F_ie + F_ee in the case of the polarizable background. The high accuracy of Eq. (33) for binary ionic mixtures in the rigid background was first demonstrated by calculations in the HNC approximation (Hansen et al. 1977; Brami et al. 1979) and confirmed later by Monte Carlo simulations (DeWitt et al. 1996; Rosenfeld 1996; DeWitt & Slattery 2003). The validity of the LMR in the case of an ionic mixture immersed in a polarizable finite-temperature electron background has been examined by Hansen et al. (1977) in the first-order thermodynamic perturbation approximation and by Chabrier & Ashcroft (1990) by solving the HNC equations with effective screened potentials. These authors found that the LMR remains accurate when the electron response is taken into account in the inter-ionic potential, as long as the Coulomb coupling is strong (Γ_j > 1, ∀j).

However, the LMR is not exact, and Eq. (33) should be replaced by the Debye – Hückel formula in the limit of weak coupling (Γ_j ≪ 1, ∀j). Even in the strong-coupling regime, the small deviations from the LMR are important for establishing phase equilibria (see Medin & Cumming 2010). The deviations from the LMR were studied by Brami et al. (1979); Chabrier & Ashcroft (1990); DeWitt et al. (1996); DeWitt & Slattery (2003); Potekhin et al. (2009a,b) for strongly coupled Coulomb liquids and by Ogata et al. (1993) and DeWitt & Slattery (2003) for Coulomb solids.

The analytical expression that describes deviations from the LMR, Δf ≡ f − f_LM, in Coulomb liquids for arbitrary coupling parameters Γ_j reads (Potekhin et al. 2009b) $\mathrm{\Delta }{\mathit{f}}_{\mathrm{liq}}\mathrm{=}\frac{{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{⟨}{\mathit{Z}}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{⟩}}{\sqrt{\mathrm{3}}}\hspace{0.17em}\frac{\mathit{\delta }}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{a}\hspace{0.17em}\mathrm{⟨}\mathrm{\Gamma }{\mathrm{⟩}}^{\mathit{\alpha }}\mathrm{)}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{b}\hspace{0.17em}\mathrm{⟨}\mathrm{\Gamma }{\mathrm{⟩}}^{\mathit{\alpha }}{\mathrm{)}}^{\mathit{\beta }}}\mathit{,}$(34)where $\mathrm{⟨}{\mathit{Z}}^{\mathit{k}}\mathrm{⟩}\mathrm{\equiv }{\sum }_{\mathit{j}}{\mathit{x}}_{\mathit{j}}{\mathit{Z}}_{\mathit{j}}^{\mathit{k}}$,  ⟨ Γ ⟩  = Γ_e ⟨ Z^5/3 ⟩ , δ is defined either as $\mathit{\delta }\mathrm{=}\mathrm{1}\mathrm{-}\frac{\mathrm{⟨}{\mathit{Z}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{⟨}\mathit{Z}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}\mathrm{⟨}{\mathit{Z}}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{⟩}}$(35)for rigid electron background model, or as $\mathit{\delta }\mathrm{=}\frac{\mathrm{⟨}\mathit{Z}\hspace{0.17em}\mathrm{(}\mathit{Z}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{⟩}}{\mathrm{⟨}{\mathit{Z}}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{⟩}}\mathrm{-}\frac{\mathrm{(}\mathrm{⟨}{\mathit{Z}}^{\mathrm{2}}\mathrm{⟩}\mathrm{+}\mathrm{⟨}\mathit{Z}\mathrm{⟩}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{⟨}\mathit{Z}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}\mathrm{⟨}{\mathit{Z}}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{⟩}}$(36)for polarizable background, and parameters a, b, α, and β depend on the plasma composition as follows: $\begin{array}{ccc}& & \mathit{a}\mathrm{=}\frac{\mathrm{2.6}\hspace{0.17em}\mathit{\delta }\mathrm{+}\mathrm{14}\hspace{0.17em}{\mathit{\delta }}^{\mathrm{3}}}{\mathrm{1}\mathrm{-}\mathit{\alpha }}\mathit{,} \mathit{\alpha }\mathrm{=}\frac{\mathrm{⟨}\mathit{Z}{\mathrm{⟩}}^{\mathrm{2}\mathit{/}\mathrm{5}}}{\mathrm{⟨}{\mathit{Z}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{1}\mathit{/}\mathrm{5}}}\mathit{,}\\ & & \mathit{b}\mathrm{=}\mathrm{0.0117}\hspace{0.17em}{\left(\frac{\mathrm{⟨}{\mathit{Z}}^{\mathrm{2}}\mathrm{⟩}}{\mathrm{⟨}\mathit{Z}{\mathrm{⟩}}^{\mathrm{2}}}\right)}^{\mathrm{2}}\mathit{a,} \mathit{\beta }\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}\mathit{\alpha }}\mathrm{-}\mathrm{1.}\end{array}$For Coulomb solids, one should distinguish regular crystals containing different ion types and disordered solid mixtures, where different ions are randomly distributed in regular lattice sites (Ogata et al. 1993). Each regular lattice type corresponds to a fixed composition, whereas random lattices allow variable fractions of different ion types. The free energy correction Δf mainly arises from the difference in the Madelung energies. It is generally larger for regular crystals than for “random” crystals with the same composition. Ogata et al. (1993) performed Monte Carlo simulations of solid ionic mixtures and fitted the calculated deviation, Δf_sol, from linear-mixing prediction for the reduced free energy in a random binary ion crystal. Medin & Cumming (2010) and Hughto et al. (2012) used this fit to study the phase separation and solidification of ion mixtures in the interiors of white dwarfs. We note, however, that the fit of Ogata et al. (1993) exhibits nonphysical features: for example, it is nonmonotonic as a function of R_Z = Z₂/Z₁ at a fixed number fraction x₂ = 1 − x₁ for a binary ion mixture with Z₂/Z₁ > 2 and x₂ < 0.5. A much simpler fit, which does not exhibit unphysical behavior, was suggested by DeWitt & Slattery (2003). It can be written as $\mathrm{\Delta }{\mathit{f}}_{\mathrm{sol}}\mathrm{=}\mathrm{0.00326}\hspace{0.17em}{\mathit{x}}_{\mathrm{1}}{\mathit{x}}_{\mathrm{2}}{\mathit{R}}_{\mathit{Z}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{⟨}\mathrm{\Gamma }\mathrm{⟩}\mathit{.}$ However, the latter fit is valid only for relatively small charge ratios R_Z ≲ 3/2. We replace it by the expression $\mathrm{\Delta }{\mathit{f}}_{\mathrm{sol}}\mathrm{=}{\mathit{x}}_{\mathrm{1}}{\mathit{x}}_{\mathrm{2}}\hspace{0.17em}{\mathrm{\Gamma }}_{\mathrm{1}}\hspace{0.17em}\mathrm{\Delta }\mathit{g}\mathrm{\left(}{\mathit{x}}_{\mathrm{2}}\mathit{,}{\mathit{R}}_{\mathit{Z}}\mathrm{\right)}\mathit{,}$(39)where $\mathrm{\Delta }\mathit{g}\mathrm{\left(}{\mathit{x}}_{\mathrm{2}}\mathit{,}{\mathit{R}}_{\mathit{Z}}\mathrm{\right)}\mathrm{=}\mathrm{0.012}\hspace{0.17em}\frac{\mathit{x}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{x}\mathrm{)}}{{\mathit{x}}_{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{x}}_{\mathrm{2}}\mathrm{)}}\mathrm{(}\mathrm{1}\mathrm{-}{{\mathit{R}}_{\mathit{Z}}^{-2}}^{\mathrm{)}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{x}}_{\mathrm{2}}\mathrm{+}{\mathit{x}}_{\mathrm{2}}{{\mathit{R}}_{\mathit{Z}}^{\mathrm{5}\mathit{/}\mathrm{3}}}^{\mathrm{)}}$(40)and $\mathit{x}\mathrm{=}{\mathit{x}}_{\mathrm{2}}\mathit{/}{\mathit{R}}_{\mathit{Z}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathit{Z}}^{-1}\mathrm{)}\hspace{0.17em}{\mathit{x}}_{\mathrm{2}}^{{\mathit{R}}_{\mathit{Z}}}\mathit{.}$ The approximation in Eq. (40) reproduces reasonably well the results of both Ogata et al. (1993) and DeWitt & Slattery (2003) for random two-component ionic bcc lattices. For a multicomponent ion crystal, Medin & Cumming (2010) proposed the extrapolation from the two-component plasma case $\mathrm{\Delta }{\mathit{f}}_{\mathrm{sol}}\mathrm{=}\sum _{\mathit{i}}\sum _{\mathit{j}\hspace{0.17em}\mathit{>}\hspace{0.17em}\mathit{i}}{\mathit{x}}_{\mathit{i}}{\mathit{x}}_{\mathit{j}}\hspace{0.17em}{\mathrm{\Gamma }}_{\mathit{i}}\hspace{0.17em}\mathrm{\Delta }\mathit{g}\left(\frac{{\mathit{x}}_{\mathit{j}}}{{\mathit{x}}_{\mathit{i}}\mathrm{+}{\mathit{x}}_{\mathit{j}}}\mathit{,}\frac{{\mathit{Z}}_{\mathit{j}}}{{\mathit{Z}}_{\mathit{i}}}\right)\mathit{,}$(41)where the indices are arranged so that Z_j < Z_j + 1.

4. EOS of a fully ionized magnetized gas

4.1. Ions

We consider only nondegenerate and nonrelativistic ions (for a discussion of the EOS of degenerate nuclear matter in strong magnetic fields see, e.g., Broderick et al. 2000; Suh & Mathews 2001). In this case (cf. Potekhin et al. 1999) $\frac{{\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}}{{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\mathrm{=}\ln \left(\mathrm{2}\mathit{\pi }\frac{{\mathit{n}}_{\mathrm{i}}{\mathit{\lambda }}_{\mathrm{i}}{\mathit{a}}_{\mathrm{m}}^{\mathrm{2}}}{\mathit{Z}}\right)\mathrm{+}\ln \mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\zeta }}_{\mathrm{i}}}\mathrm{)}\mathrm{-}\mathrm{1}\mathrm{+}\frac{{\mathit{\zeta }}_{\mathrm{i}}}{\mathrm{2}}\mathrm{+}\frac{{\mathit{F}}_{\mathrm{spin}}}{{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\mathrm{·}$(42)The last term arises from the energy of the magnetic moments of the ions in the magnetic field, ${\mathit{F}}_{\mathrm{spin}}\mathrm{=}\mathrm{-}{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}\ln [\frac{\sinh \mathrm{(}{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}{\mathit{M}}_{\mathrm{spin}}\mathit{/}\mathrm{4}\mathrm{)}}{\sinh \mathrm{(}{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}\mathit{/}\mathrm{4}\mathrm{)}}]\mathit{,}$(43)where M_spin is the ion spin multiplicity, and g_i is the g-factor (g_i = 5.5857 for protons). For ions with spin one-half (M_spin = 2), the expression in the square brackets in Eq. (43) simplifies to  [2cosh(g_i   ζ_i/4)] . For zero-spin ions, such as ⁴He, ⁵⁶Fe, and other even-even nuclei in the ground state, F_spin = 0.

The ion pressure obeys the nonmagnetic ideal-gas relation ${\mathit{P}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}\mathrm{=}{\mathit{n}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}$, but expressions for the internal energy and heat capacity are different: $\begin{array}{ccc}\frac{{\mathit{U}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}}{{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\frac{{\mathit{\zeta }}_{\mathrm{i}}}{{\mathit{e}}^{{\mathit{\zeta }}_{\mathrm{i}}}\mathrm{-}\mathrm{1}}\mathrm{+}\frac{{\mathit{\zeta }}_{\mathrm{i}}}{\mathrm{2}}\mathrm{+}{\mathit{u}}_{\mathrm{spin}}\mathit{,}\\ \frac{{\mathit{C}}_{\mathit{V,}\mathrm{id}}^{\mathrm{\left(}\mathrm{i}\mathrm{\right)}}}{{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}{\left(\frac{{\mathit{\zeta }}_{\mathrm{i}}}{{\mathit{e}}^{{\mathit{\zeta }}_{\mathrm{i}}}\mathrm{-}\mathrm{1}}\right)}^{\mathrm{2}}\mathrm{+}{\mathit{c}}_{\mathrm{spin}}\mathit{.}\end{array}$Here, the terms u_spin and c_spin arise from F_spin, $\begin{array}{ccc}{\mathit{u}}_{\mathrm{spin}}& \mathrm{=}& \frac{{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}\mathit{/}\mathrm{4}}{\tanh \mathrm{(}{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}\mathit{/}\mathrm{4}\mathrm{)}}\mathrm{-}\frac{{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}{\mathit{M}}_{\mathrm{spin}}\mathit{/}\mathrm{4}}{\tanh \mathrm{(}{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}{\mathit{M}}_{\mathrm{spin}}\mathit{/}\mathrm{4}\mathrm{)}}\mathit{,}\\ {\mathit{c}}_{\mathrm{spin}}& \mathrm{=}& {\left(\frac{{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}\mathit{/}\mathrm{4}}{\sinh \mathrm{(}{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}\mathit{/}\mathrm{4}\mathrm{)}}\right)}^{\mathrm{2}}\mathrm{-}{\left(\frac{{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}{\mathit{M}}_{\mathrm{spin}}\mathit{/}\mathrm{4}}{\sinh \mathrm{(}{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}{\mathit{M}}_{\mathrm{spin}}\mathit{/}\mathrm{4}\mathrm{)}}\right)}^{\mathrm{2}}\mathrm{·}\end{array}$They simplify at M_spin = 2: ${\mathit{u}}_{\mathrm{spin}}\mathrm{=}\mathrm{-}\frac{{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}}{\mathrm{4}}\tanh (\frac{{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}}{\mathrm{4}})\mathit{,} {\mathit{c}}_{\mathrm{spin}}\mathrm{=}{\left(\frac{{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}\mathit{/}\mathrm{4}}{\cosh \mathrm{(}{\mathit{g}}_{\mathrm{i}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{i}}\mathit{/}\mathrm{4}\mathrm{)}}\right)}^{\mathrm{2}}\mathrm{·}$(48)

4.2. Electrons

4.2.1. General case

Thermodynamic functions of the electron gas in a magnetic field are easily derived from first principles (Landau & Lifshitz 1980). The number of quantum states per longitudinal momentum interval Δp_z for an electron with given B-projections of the spin and the orbital moment and with a fixed Landau number n in a volume V equals $\mathit{V}\mathrm{\Delta }{\mathit{p}}_{\mathit{z}}\mathit{/}\mathrm{\left(}\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}{\mathit{a}}_{\mathrm{m}}^{\mathrm{2}}ħ\mathrm{\right)}$ (Landau & Lifshitz 1977). Thus one can express the electron number density n_e and the grand thermodynamic potential ${\mathrm{\Omega }}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{=}\mathrm{-}{\mathit{P}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathit{V}$ as $\begin{array}{ccc}{\mathit{n}}_{\mathrm{e}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}{\mathit{a}}_{\mathrm{m}}^{\mathrm{2}}ħ}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{\sigma }}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\frac{\mathrm{d}{\mathit{p}}_{\mathit{z}}}{\exp \mathrm{\left[}\mathrm{\right(}{\mathit{ϵ}}_{\mathit{n}}\mathrm{\left(}{\mathit{p}}_{\mathit{z}}\mathrm{\right)}\mathrm{-}{\mathit{\mu }}_{\mathrm{e}}\mathrm{)}\mathit{/}{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathrm{)}\mathrm{]}\mathrm{+}\mathrm{1}}\mathit{,}\\ {\mathrm{\Omega }}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}& \mathrm{=}& \mathrm{-}\frac{\mathit{V}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}{\mathit{a}}_{\mathrm{m}}^{\mathrm{2}}ħ}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{\sigma }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\ln \left(\mathrm{1}\mathrm{+}\exp \left[\frac{{\mathit{\mu }}_{\mathrm{e}}\mathrm{-}{\mathit{ϵ}}_{\mathit{n}}\mathrm{\left(}{\mathit{p}}_{\mathit{z}}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{B}}\mathit{T}}\right]\right)\mathrm{d}{\mathit{p}}_{\mathit{z}}\mathit{,}\end{array}$where ϵ_n(p_z) is given by Eq. (7) and  ∑ _σ denotes the sum over spin projections, which amounts to the factor 2 for n ≥ 1 since we neglect the anomalous magnetic moment of electrons. This derivation equally holds in the relativistic and nonrelativistic theories. Equations (49) and (50) can be rewritten, using integration by parts, as $\begin{array}{ccc}{\mathit{n}}_{\mathrm{e}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\pi }}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{a}}_{\mathrm{m}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{e}}}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{\sigma }}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathit{bn}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{4}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{I}}_{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\left(}{\mathit{\chi }}_{\mathit{n}}\mathit{,}{\mathit{\tau }}_{\mathit{n}}\mathrm{\right)}}{\mathit{\partial }{\mathit{\chi }}_{\mathit{n}}}\mathit{,}& & \\ {\mathit{P}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{{\mathit{\pi }}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{a}}_{\mathrm{m}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{e}}}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{\sigma }}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathit{bn}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{4}}\hspace{0.17em}{\mathit{I}}_{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}{\mathit{\chi }}_{\mathit{n}}\mathit{,}{\mathit{\tau }}_{\mathit{n}}\mathrm{)}\mathit{,}& & \end{array}$where ${\mathit{\tau }}_{\mathit{n}}\mathrm{=}\mathit{\tau }\mathit{/}\sqrt{\mathrm{1}\mathrm{+}\mathrm{2}\mathit{bn}}$ and ${\mathit{\chi }}_{\mathit{n}}\mathrm{=}{\mathit{\chi }}_{\mathrm{e}}\mathrm{+}{\mathit{\tau }}^{-1}\mathrm{-}{\mathit{\tau }}_{\mathit{n}}^{-1}\mathit{.}$ The free energy ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}$ is given by Eqs. (17), (51), and (52).

The calculation of n_e, ${\mathit{P}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}$, and their derivatives at given χ_e and τ can be performed using Eqs. (51, (52) and the same analytical approximations to the Fermi-Dirac integrals as for the nonmagnetized electron gas. The reduced electron chemical potential χ_e at constant n_e and T is found by numerical inversion of Eq. (51). Then the derivatives over T at constant V and over V at constant T are given by Eqs. (22) and (23). We use this approach in the current research, but we should note that for quantizing magnetic fields it is less precise than at B = 0. As mentioned in Sect. 3.1, the inaccuracy of the employed approximations for I_ν(χ_e,τ) is within a fraction of percent, but it grows for the derivatives. Since the first derivatives are already employed in Eq. (51), evaluation of the second-order thermodynamic functions such as χ_T or C_V involves third derivatives. Therefore, the error in the evaluation of these functions may rise to several percent. This level of accuracy may be sufficient for astrophysical applications, but otherwise one should resort to a thermodynamically consistent interpolation in numerical tables of the Fermi-Dirac integrals (Timmes & Arnett 1999; Timmes & Swesty 2000). Equations (51) and (52) can be simplified in several limiting cases considered below.

4.2.2. Strongly quantizing and nonquantizing magnetic fields

The field is called strongly quantizing if most of the electrons reside on the ground Landau level. The electron Fermi momentum, then, equals ${\mathit{p}}_{\mathrm{F}}\mathrm{=}\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}{\mathit{a}}_{\mathrm{m}}^{\mathrm{2}}ħ{\mathit{n}}_{\mathrm{e}}\mathrm{=}\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}\hspace{0.17em}\mathrm{(}{\mathit{a}}_{\mathrm{m}}\mathit{/}{\mathit{a}}_{\mathrm{e}}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}{\mathit{p}}_{\mathrm{F}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathit{,}$(53)where ${\mathit{p}}_{\mathrm{F}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ is the zero-field Fermi momentum at the given density. Equation (53) can be written as p_F = m_ecx_B, where ${\mathit{x}}_{\mathit{B}}\mathrm{=}\frac{\mathrm{2}{\mathit{x}}_{\mathrm{r}}^{\mathrm{3}}}{\mathrm{3}\mathit{b}}\mathrm{\approx }\mathrm{30.2}\hspace{0.17em}\frac{\mathit{Z}}{\mathit{A}}\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{6}}}{{\mathit{B}}_{\mathrm{12}}}$(54)is the relativity parameter modified by the field and ${\mathit{x}}_{\mathrm{r}}\mathrm{=}{\mathit{p}}_{\mathrm{F}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathit{/}{\mathit{m}}_{\mathrm{e}}\mathit{c}$ (Sect. 2.1). With increasing n_e at constant B and zero temperature, the electrons start to populate the first excited Landau level when n_e reaches ${\mathit{n}}_{\mathit{B}}\mathrm{=}\mathrm{(}{\mathit{\pi }}^{\mathrm{2}}\sqrt{\mathrm{2}}\hspace{0.17em}{\mathit{a}}_{\mathrm{m}}^{\mathrm{3}}{\mathrm{)}}^{-1}$. Therefore, the field is strongly quantizing at T ≪ T_cycl and ρ < ρ_B, where T_cycl = ħω_c/k_B ≈ 1.3434 × 10⁸B₁₂ K and ${\mathit{\rho }}_{\mathit{B}}\mathrm{=}{\mathit{m}}_{\mathrm{i}}{\mathit{n}}_{\mathit{B}}\mathit{/}\mathit{Z}\mathrm{\approx }\mathrm{7045}\hspace{0.17em}\mathrm{(}\mathit{A}\mathit{/}\mathit{Z}\mathrm{)}\hspace{0.17em}{\mathit{B}}_{\mathrm{12}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{g} \mathrm{cm}-3\mathit{.}$(55)The condition n_e < n_B can be written as ${\mathit{a}}_{\mathrm{m}}\mathit{/}{\mathit{a}}_{\mathrm{e}}\mathit{<}\sqrt{\mathrm{2}}\hspace{0.17em}\mathrm{(}\mathrm{3}\mathit{\pi }{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}$. Then Eq. (53) shows that in a strongly quantizing field ${\mathit{p}}_{\mathrm{F}}\mathit{<}{\mathit{p}}_{\mathrm{F}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$, except for densities n_e close to the threshold n_B. Thus T_F is reduced, compared to its nonmagnetic value ${\mathit{T}}_{\mathrm{F}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$, by factor ${\mathit{T}}_{\mathrm{F}}\mathit{/}{\mathit{T}}_{\mathrm{F}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{=}\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{]}\mathit{/}\mathrm{(}{\mathit{\gamma }}_{\mathrm{r}}\mathrm{-}\mathrm{1}\mathrm{)}$. In the nonrelativistic limit, ${\mathit{T}}_{\mathrm{F}}\mathit{/}{\mathit{T}}_{\mathrm{F}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{=}\mathrm{(}{\mathit{p}}_{\mathrm{F}}\mathit{/}{\mathit{p}}_{\mathrm{F}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}{\mathrm{)}}^{\mathrm{2}}$, and the parameter θ = T/T_F becomes ${\mathit{\theta }}_{\mathit{B}}\mathrm{=}\mathrm{8}{\mathit{\gamma }}_{\mathrm{m}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{s}}^{\mathrm{5}}\mathit{/}\mathrm{\left(}\mathrm{9}{\mathit{\pi }}^{\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{e}}\mathrm{\right)}\mathrm{\approx }\mathrm{0.166}\hspace{0.17em}{\mathit{\theta }}_{\mathrm{0}}\hspace{0.17em}{\mathit{\gamma }}_{\mathrm{m}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{s}}^{\mathrm{4}}\mathit{,}$(56)where θ₀ is the nonmagnetic value given by Eq. (2).

The opposite case of a nonquantizing magnetic field occurs at T ≫ T_B, where T_B is the temperature at which the thermal kinetic energy of the electrons becomes sufficient to smear their distribution over many Landau levels. It can be estimated as T_B = T_cycl, if ρ ≤ ρ_B and T_B = T_cycl/γ_r, if ρ > ρ_B (a more sophisticated but qualitatively similar definition of T_B was introduced by Lai 2001). Then we can approximately replace the sum over Landau level numbers n by the integral over a continuous variable n. Integrating over n by parts, we can reduce Eq. (52) to Eq. (18) and Eq. (51) to Eq. (19), i.e., to recover the zero-field thermodynamics. At ρ ≫ ρ_B, the electrons also fill many Landau levels and the magnetic field can be approximately treated as nonquantizing.

In the intermediate region, where the magnetic field is neither strongly quantizing nor nonquantizing, the summation over n manifests itself in quantum oscillations of the thermodynamic functions with changing B and/or ρ, similar to the de Haas – van Alphen oscillations of magnetic susceptibility (e.g., Landau & Lifshitz 1980). The oscillations are smoothed by the thermal broadening of the Fermi distribution function and by the quantum broadening of the Landau levels (particularly, owing to electron collisions; see Yakovlev & Kaminker 1994, for references). Some examples of such oscillations will be given in Sect. 7.

Figure 1 presents the ρ – T diagram of outer neutron-star envelopes at two magnetic field strengths, B = 10¹² G and 10¹⁵ G, assuming fully-ionized iron (this assumption may be crude in the lower left part of the diagram). In the strongly-quantizing magnetic-field domain, bounded by ρ_B and T_cycl, the dependence T_F(ρ) is steeper than at B = 0, in agreement with Eq. (56). The line T_m(ρ) separates Coulomb liquid from Coulomb crystal. Near the lower right corner of the figure, where T ≪ T_p, the quantum effects on the ions become important (i.e., the ions cannot be treated as classical pointlike particles). In the lower left corner, at ρ < ρ_s, the plasma is unstable to the phase separation into the gaseous and condensed phases (this phase transition will be discussed in Sect. 7.3).

4.2.3. Strongly degenerate electrons

If the electrons are strongly degenerate, then one can apply the Sommerfeld expansion (Sect. 3.1) and obtain ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{\approx }{\mathit{F}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{+}\mathrm{\Delta }\mathit{F}$, where ${\mathit{F}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{=}{\mathit{ϵ}}_{\mathrm{F}}{\mathit{N}}_{\mathrm{e}}\mathrm{-}{\mathit{P}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathit{V}$ is the zero-temperature value and ΔF is a thermal correction. According to Eqs. (24) and (52), the zero-temperature pressure P₀ is ${\mathit{P}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{r}}\hspace{0.17em}\mathit{b}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{{\mathit{n}}_{\max }}\sum _{\mathit{\sigma }}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathit{bn}\mathrm{)}\hspace{0.17em}{\mathrm{ℐ}}_{\mathrm{1}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ {\mathit{ϵ}}_{\mathit{n}}\end{array}\mathrm{\right)}\mathit{,}$(57)where ${\mathit{P}}_{\mathrm{r}}\mathrm{\equiv }{\mathit{m}}_{\mathrm{e}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{3}}\mathit{/}{ħ}^{\mathrm{3}}\mathrm{=}\mathrm{1.4218}\mathrm{\times }{\mathrm{10}}^{\mathrm{25}}{ \mathrm{dyn} \mathrm{cm}}^{-2}$ is the relativistic unit of pressure, n_max is the maximum integer n for which ${\mathit{p}}_{\mathit{n}}^{\mathrm{2}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathit{>}\mathrm{0}$, and $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{ϵ}}_{\mathit{n}}\end{array}\mathrm{\equiv }{\mathit{ϵ}}_{\mathrm{F}}\mathit{/}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{-}\sqrt{\mathrm{1}\mathrm{+}\mathrm{2}\mathit{bn}}$. According to Eqs. (24) and (51), the Fermi energy ϵ_F is determined by the condition ${\mathit{n}}_{\mathrm{e}}\mathrm{=}{\left(\frac{{\mathit{m}}_{\mathrm{e}}\mathit{c}}{ħ}\right)}^{\mathrm{3}}\frac{\mathit{b}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{{\mathit{n}}_{\max }}\sum _{\mathit{\sigma }}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathit{bn}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathrm{ℐ}}_{\mathrm{1}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{(}\begin{array}{}\mathrm{˜}\\ {\mathit{ϵ}}_{\mathit{n}}\end{array}\mathrm{)}\mathit{.}$(58)In order to obtain the chemical potential μ_e = ϵ_F + Δϵ with fractional accuracy  ~ ${\mathit{\chi }}_{\mathrm{e}}^{-2}$, we retain two terms in Eq. (24), insert it into Eq. (51), approximate ${\mathrm{ℐ}}_{\mathit{\nu }}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{\mu }}_{\mathit{n}}\end{array}\mathrm{\right)}$ in the vicinity of $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{\mu }}_{\mathit{n}}\end{array}\mathrm{=}{\mathit{ϵ}}_{\mathit{n}}$ by ${\mathrm{ℐ}}_{\mathit{\nu }}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{n}}\end{array}\mathrm{\right)}\mathrm{\approx }{\mathrm{ℐ}}_{\mathit{\nu }}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ {\mathit{ϵ}}_{\mathit{n}}\end{array}\mathrm{\right)}\mathrm{+}{\mathrm{ℐ}}_{\mathit{\nu }}^{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ {\mathit{ϵ}}_{\mathit{n}}\end{array}\mathrm{\right)}\hspace{0.17em}\mathrm{\Delta }\begin{array}{}\mathrm{˜}\\ \mathit{ϵ}\end{array}\mathit{,}$(59)where $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{\mu }}_{\mathit{n}}\end{array}\mathrm{\equiv }{\mathit{\chi }}_{\mathit{n}}\mathit{\tau }$ and $\mathrm{\Delta }\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{ϵ}\end{array}\mathrm{\equiv }\mathrm{\Delta }\mathit{ϵ}\mathit{/}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}$, and drop the higher-order terms containing $\mathrm{\left(}{\mathit{\tau }}^{\mathrm{2}}\mathrm{\Delta }\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{ϵ}\end{array}\mathrm{\right)}$. Then $\mathrm{\Delta }\begin{array}{}\mathrm{˜}\\ \mathit{ϵ}\end{array}\mathrm{\approx }\mathrm{-}\frac{{\mathit{\pi }}^{\mathrm{2}}{\mathit{\tau }}^{\mathrm{2}}}{\mathrm{6}}\frac{{\sum }_{\mathit{n}\mathrm{=}\mathrm{0}}^{{\mathit{n}}_{\max }}{\sum }_{\mathit{\sigma }}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathit{bn}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathrm{ℐ}}_{\mathrm{1}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}\mathrm{(}\begin{array}{}\mathrm{˜}\\ {\mathit{ϵ}}_{\mathit{n}}\end{array}\mathrm{)}}{{\sum }_{\mathit{n}\mathrm{=}\mathrm{0}}^{{\mathit{n}}_{\max }}{\sum }_{\mathit{\sigma }}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathit{bn}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathrm{ℐ}}_{\mathrm{1}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{(}\begin{array}{}\mathrm{˜}\\ {\mathit{ϵ}}_{\mathit{n}}\end{array}\mathrm{)}}\mathrm{·}$(60)The thermal correction to the pressure equals $\mathrm{\Delta }\mathit{P}\mathrm{=}{\mathit{P}}_{\mathrm{r}}\hspace{0.17em}\mathit{b}\sum _{\mathit{n}\mathrm{=}\mathrm{0}}^{{\mathit{n}}_{\max }}\sum _{\mathit{\sigma }}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathit{bn}\mathrm{)}\left(\frac{{\mathit{\tau }}^{\mathrm{2}}}{\mathrm{12}}\hspace{0.17em}{\mathrm{ℐ}}_{\mathrm{1}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ {\mathit{ϵ}}_{\mathit{n}}\end{array}\mathrm{\right)}\mathrm{+}\frac{\mathrm{\Delta }\begin{array}{}\mathrm{˜}\\ \mathit{ϵ}\end{array}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}\hspace{0.17em}{\mathrm{ℐ}}_{\mathrm{1}\mathit{/}\mathrm{2}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{\left(}\begin{array}{}\mathrm{˜}\\ {\mathit{ϵ}}_{\mathit{n}}\end{array}\mathrm{\right)}\right)\mathit{,}$(61)and the thermal correction to the free energy and internal energy $\mathrm{\Delta }\mathit{F}\mathrm{=}\mathrm{-}\mathrm{\Delta }\mathit{U}\mathrm{=}{\mathit{N}}_{\mathrm{e}}\hspace{0.17em}\mathrm{\Delta }\mathit{ϵ}\mathrm{-}\mathit{V}\mathrm{\Delta }\mathit{P}\mathit{.}$(62)The leading contribution to the heat capacity is ${\mathit{C}}_{\mathit{V}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{=}\mathrm{2}\mathrm{\Delta }\mathit{U}\mathit{/}\mathit{T}$. As in the nonmagnetic case, ${\mathit{C}}_{\mathit{V}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}$ is proportional to T at T → 0, but with a different proportionality coefficient.

4.2.4. Strongly degenerate electrons in a strongly quantizing magnetic field

If the magnetic field is strongly quantizing and the electrons are strongly degenerate (which corresponds to the triangular domains in Fig. 1 defined by conditions ρ < ρ_B and T < T_F), then $\begin{array}{ccc}{\mathit{F}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}& \mathrm{=}& \mathrm{[}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{]}}{\mathit{N}}_{\mathrm{e}}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}\mathrm{-}{\mathit{P}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathit{V,}\\ {\mathit{P}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}& \mathrm{=}& \frac{{\mathit{P}}_{\mathrm{r}}\hspace{0.17em}\mathit{b}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{[}{\mathit{x}}_{\mathit{B}}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}\ln {\mathrm{(}{\mathit{x}}_{\mathit{B}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}}^{\mathrm{]}}\mathit{.}\end{array}$In the nonrelativistic (x_B ≪ 1) and ultrarelativistic (x_B ≫ 1) limits, we have ${\mathit{P}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{\simeq }{\mathit{P}}_{\mathrm{r}}\hspace{0.17em}\mathit{b}\hspace{0.17em}{\mathit{x}}_{\mathit{B}}^{\mathrm{3}}\mathit{/}\mathrm{6}{\mathit{\pi }}^{\mathrm{2}}\mathrm{\propto }{\mathit{n}}_{\mathrm{e}}^{\mathrm{3}}$ and ${\mathit{P}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{\simeq }{\mathit{P}}_{\mathrm{r}}\hspace{0.17em}\mathit{b}\hspace{0.17em}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}\mathit{/}\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}\mathrm{\propto }{\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}$, respectively. Compared with the nonmagnetic case (Papers I and II), the dependence ${\mathit{P}}_{\mathrm{0}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}\mathrm{\left(}{\mathit{n}}_{\mathrm{e}}\mathrm{\right)}$ is steeper, but P is lower everywhere except for n_e ≈ n_B. Thus, a strongly quantizing magnetic field softens the EOS of degenerate electrons.

The thermal corrections (Eqs. (60)–(62)) simplify to $\begin{array}{ccc}& & \mathrm{\Delta }\begin{array}{}\mathrm{˜}\\ \mathit{ϵ}\end{array}\mathrm{\approx }\mathrm{-}\frac{{\mathit{\pi }}^{\mathrm{2}}{\mathit{\tau }}^{\mathrm{2}}}{\mathrm{6}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}}\mathit{,}\\ & & \mathrm{\Delta }\mathit{P}\mathrm{=}{\mathit{P}}_{\mathrm{r}}\hspace{0.17em}\frac{\mathit{b}{\mathit{\tau }}^{\mathrm{2}}}{\mathrm{12}}\hspace{0.17em}\frac{\mathrm{2}\mathrm{+}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{x}}_{\mathit{B}}^{}}\mathit{,}\\ & & \frac{\mathrm{\Delta }\mathit{U}}{\mathit{V}}\mathrm{=}\mathrm{-}\frac{\mathrm{\Delta }\mathit{F}}{\mathit{V}}\mathrm{=}{\mathit{P}}_{\mathrm{r}}\hspace{0.17em}\frac{\mathit{b}{\mathit{\tau }}^{\mathrm{2}}}{\mathrm{12}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{{\mathit{x}}_{\mathit{B}}}\mathit{,}\\ & & \frac{{\mathit{C}}_{\mathit{V}}^{\mathrm{\left(}\mathrm{e}\mathrm{\right)}}}{{\mathit{N}}_{\mathrm{e}}{\mathit{k}}_{\mathrm{B}}}\mathrm{=}\frac{{\mathit{\pi }}^{\mathrm{2}}\mathit{\tau }}{\mathrm{3}}\hspace{0.17em}\frac{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{{\mathit{x}}_{\mathit{B}}^{\mathrm{2}}}\mathrm{·}\end{array}$The last equation differs from the nonmagnetic case (Paper II) in that x_B replaces x_r and π²/3 replaces π².

4.2.5. Nonrelativistic limit

In the nonrelativistic limit (p_F ≪ m_ec and T ≪ T_r), Eqs. (51) and (52) simplify to ${\mathit{n}}_{\mathrm{e}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{a}}_{\mathrm{m}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{e}}}\sum _{\mathit{n,\sigma }}{\mathit{I}}_{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\left(}{\mathit{\chi }}_{\mathit{n}}\mathrm{\right)}\mathit{,} {\mathit{P}}_{\mathrm{e}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{{\mathit{\pi }}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{a}}_{\mathrm{m}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{e}}}\sum _{\mathit{n,\sigma }}{\mathit{I}}_{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\left(}{\mathit{\chi }}_{\mathit{n}}\mathrm{\right)}\mathit{.}$(69)In the nondegenerate regime (T ≫ T_F), one has I_ν(χ) ≈ e^χ   Γ(ν + 1), where Γ(ν + 1) is the gamma-function. Then Eq. (69) yields ${\mathit{P}}_{\mathrm{id}}^{\mathrm{e}}\mathrm{=}{\mathit{n}}_{\mathrm{e}}{\mathit{k}}_{\mathrm{B}}\mathit{T}$ and ${\mathit{\chi }}_{\mathrm{e}}\mathrm{=}\ln \mathrm{(}{\mathit{n}}_{\mathrm{e}}{\mathit{\lambda }}_{\mathrm{e}}^{\mathrm{3}}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{-}\ln \mathrm{(}\mathit{\zeta }\mathit{/}\mathrm{4}\mathrm{)}\mathrm{+}\ln \mathrm{\left[}\tanh \mathrm{\right(}\mathit{\zeta }\mathit{/}\mathrm{4}\mathrm{\left)}\mathrm{\right]}$(70)which provides the free energy ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathit{e}\mathrm{\right)}}\mathrm{=}\mathrm{(}{\mathit{\chi }}_{\mathrm{e}}\mathrm{-}\mathrm{1}\mathrm{)}\hspace{0.17em}{\mathit{N}}_{\mathrm{e}}{\mathit{k}}_{\mathrm{B}}\mathit{T}$.

As follows from Eq. (70), the reduced internal energy ${\mathit{U}}_{\mathrm{id}}^{\mathrm{e}}\mathit{/}{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}$ and heat capacity ${\mathit{C}}_{\mathit{V,}\mathrm{id}}^{\mathrm{e}}\mathit{/}{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}$ of the Boltzmann gas decrease with increasing ζ. In a strongly quantizing magnetic field (ζ ≫ 1), they tend to 1/2 instead of 3/2 because the gas becomes effectively one-dimensional. The only kinetic degree is along the magnetic field.

In the nonquantizing field (ζ ≪ 1), the two last terms in Eq. (70) cancel out, so that the standard expression ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathit{e}\mathrm{\right)}} \mathrm{=} {\mathit{N}}_{\mathrm{e}}{\mathit{k}}_{\mathrm{B}}\mathit{T}\hspace{0.17em}\mathrm{\left[}\ln \mathrm{\right(}{\mathit{n}}_{\mathrm{e}}{\mathit{\lambda }}_{\mathrm{e}}^{\mathrm{3}}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{-}\mathrm{1}\mathrm{]}$ is recovered. In the strongly quantizing, nondegenerate regime (ρ < ρ_B and T_F ≪ T ≪ T_cycl), the last term of Eq. (70) vanishes, which yields ${\mathit{F}}_{\mathrm{id}}^{\mathrm{\left(}\mathit{e}\mathrm{\right)}}\mathrm{=}{\mathit{N}}_{\mathrm{e}}{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathrm{\left[}\ln \mathrm{\right(}\mathrm{2}\mathit{\pi }{\mathit{a}}_{\mathrm{m}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{e}}{\mathit{n}}_{\mathrm{e}}\mathrm{)}\mathrm{-}{\mathrm{1}}^{\mathrm{]}}\mathit{.}$(71)

4.3. Thermodynamic and kinetic pressures

The above expressions for pressure of a magnetized gas of charged particles are based on the principles of thermodynamics (Landau & Lifshitz 1980), according to which P = −(∂F/∂V)_T,B. Alternatively, the pressure can be calculated from the microscopic dynamics as the sum of the changes of kinetic momenta of all particles reflected off a unit surface per unit time. The result of the latter calculation, the kinetic pressure P^kin, depends on the orientation of the surface relative to B (Canuto & Chiu 1968). If the surface is perpendicular to B, then one gets the kinetic pressure ${\mathit{P}}_{\mathrm{\parallel }}^{\mathrm{kin}}\mathrm{=}\mathit{P}$, which acts along the field lines, the longitudinal pressure. If the surface is parallel to the field, one gets a different (transverse) kinetic pressure, which can be expressed (Blandford & Hernquist 1982) as ${\mathit{P}}_{\mathrm{\perp }}^{\mathrm{kin}}\mathrm{=}\mathrm{-}\mathrm{\Omega }\mathit{/}\mathit{V}\mathrm{+}\mathit{B}\hspace{0.17em}\mathrm{(}\mathit{\partial }\mathrm{\Omega }\mathit{/}\mathit{\partial B}{\mathrm{)}}_{\mathit{V,T}}\mathrm{=}\mathit{P}\mathrm{-}\mathit{MB,}$(72)where M = −∂Ω/∂B is the magnetization.

In order to resolve the apparent paradox, one should take the magnetization current density j_m = c   ∇ × M into account; when boundaries are present, this volume current should be supplemented by the surface current cM × B/B (see, e.g., Griffith 1999). As argued by Blandford & Hernquist (1982), if we compress the electron gas perpendicular to B then we must do work against the Lorentz force density j_m × B/c, which gives an additional contribution to the total transverse pressure and makes it equal to P. Because this point still causes confusion in some publications, let us illustrate it with a simple example.

If the pressure were anisotropic, then one might expect an anisotropic density gradient in a strongly magnetized star. Let us consider a small volume element in the star, assuming that we can treat B, T, and gravitational acceleration g as constants within this volume, and z axis is directed along g. Hydrostatic equilibrium implies that the density of gravitational force, ρg, is balanced by the density of forces created by plasma particles. The crucial point is that the magnetization contributes to this balance.

Let us compare the cases where B is parallel and perpendicular to g. In the first case, the z-component of the Lorentz force is absent, and we get the standard equation of hydrostatic equilibrium: $\mathit{\rho g}\mathrm{=}\mathrm{d}{\mathit{P}}_{\mathrm{\parallel }}^{\mathrm{kin}}\mathit{/}\mathrm{d}\mathit{z}\mathrm{=}\mathrm{d}\mathit{P}\mathit{/}\mathrm{d}\mathit{z}$. In the second case, the gradients of the kinetic pressure $\mathrm{d}{\mathit{P}}_{\mathrm{\perp }}^{\mathrm{kin}}\mathit{/}\mathrm{d}\mathit{z}$ and of the Lorentz force density B   dM/dz act in parallel. In the constant and uniform magnetic field, dM/dz is not zero, but is related to the density gradient: $\frac{\mathrm{d}\mathit{M}}{\mathrm{d}\mathit{z}}\mathrm{=}\frac{\mathit{\partial M}\mathrm{\left(}\mathit{\rho ,T,B}\mathrm{\right)}}{\mathit{\partial \rho }}\hspace{0.17em}\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{z}}\mathrm{=}\mathrm{-}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Omega }\mathrm{\left(}\mathit{\rho ,T,B}\mathrm{\right)}}{\mathit{V\partial \rho }\hspace{0.17em}\mathit{\partial B}}\hspace{0.17em}\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{z}}\mathrm{·}$(73)Then the equilibrium condition takes the form $\mathit{\rho g}\mathrm{=}\mathrm{d}{\mathit{P}}_{\mathrm{\perp }}^{\mathrm{kin}}\mathit{/}\mathrm{d}\mathit{z}\mathrm{+}\mathit{B}\hspace{0.17em}\mathrm{d}\mathit{M}\mathit{/}\mathrm{d}\mathit{z}\mathrm{=}\mathrm{d}\mathit{P}\mathit{/}\mathrm{d}\mathit{z}\mathit{,}$ which is the same as in the first case. Furthermore, one can express P through ρ using some EOS. In the considered example, dP/dz =  [∂P(ρ,T,B)/∂ρ]    dρ/dz, so that dρ/dz is also the same in both cases. Thus, the stellar hydrostatic profile is determined by the isotropic thermodynamic pressure P, which automatically includes magnetization.

5. Magnetic effects on the EOS of a Coulomb liquid

5.1. Electron exchange and correlation

The effects of a magnetic field on the contribution to the free energy due to electron exchange and correlation were studied either in the regime of strong degeneracy and strongly quantizing magnetic fields (Danz & Glasser 1971; Banerjee et al. 1974; Fushiki et al. 1989; see also Morbec & Capelle 2008 for an instructive discussion of the previous results and the inclusion of the second Landau level contribution), or at low densities (Alastuey & Jancovici 1980; Cornu 1998; Steinberg et al. 1998, 2000). In a previous work (Potekhin et al. 1999) we suggested a modification of the field-free expression for F_ee, which matches available exact limiting expressions, including the cases of nonquantizing, strongly quantizing degenerate, and strongly quantizing nondegenerate plasmas. The modification consists in replacing F_ee(θ,Γ_e) by F_ee(θ ^∗ ,Γ_e), where ${\mathit{\theta }}^{\mathrm{\ast }}\mathrm{=}\frac{{\mathit{\theta }}_{\mathrm{0}}\mathrm{+}{\mathit{\theta }}_{\mathit{B}}}{\mathrm{1}\mathrm{+}\frac{{\mathit{\theta }}_{\mathit{B}}}{{\mathit{\theta }}_{\mathrm{0}}}\hspace{0.17em}\exp \mathrm{(}\mathrm{-}{\mathit{\theta }}_{\mathit{B}}^{-1}\mathrm{)}\hspace{0.17em}\frac{\cosh \mathrm{(}\mathit{\zeta }\mathit{/}\mathrm{2}\mathrm{)}}{\mathrm{\left[}\cosh \mathrm{\right(}\mathit{\zeta }\mathit{/}\mathrm{4}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}\hspace{0.17em}\frac{\tanh \mathrm{(}\mathit{\zeta }\mathit{/}\mathrm{4}\mathrm{)}}{\mathit{\zeta }\mathit{/}\mathrm{4}}\hspace{0.17em}\frac{\mathrm{arctanh}\hspace{0.17em}\mathit{\xi }}{\mathit{\xi }}}\hspace{0.17em}\mathit{,}$(74)ξ =  [1−(4/ζ)tanh(ζ/4)] ^1/2, and θ₀ and θ_B are given by Eqs. (2) and (56), respectively, at fixed n_e and T.

5.2. Wigner-Kirkwood term

For the same reason as in Sect. 3.2.1, the treatment of the quantum effects in the ion liquid is restricted by the Wigner-Kirkwood term. Its expression in an arbitrary magnetic field was obtained by Alastuey & Jancovici (1980): ${\mathit{f}}_{\mathrm{ii}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{\eta }}^{\mathrm{2}}}{\mathrm{24}}\hspace{0.17em}\left[\frac{\mathrm{4}}{{\mathit{\zeta }}_{\mathrm{i}}\hspace{0.17em}\tanh \mathrm{(}{\mathit{\zeta }}_{\mathrm{i}}\mathit{/}\mathrm{2}\mathrm{)}}\mathrm{-}\frac{\mathrm{8}}{{\mathit{\zeta }}_{\mathrm{i}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}}\right]\mathrm{·}$(75)The function in the square brackets monotonously varies from 1 at ζ_i → 0 to 1/3 at ζ_i → ∞, reflecting the effective reduction of the degrees of freedom of the classical ion motion from d = 3 at B = 0 to d = 1 for a strongly quantizing field. At small ζ_i, ${\mathit{f}}_{\mathrm{ii}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\approx }\mathrm{(}{\mathit{\eta }}^{\mathrm{2}}\mathit{/}\mathrm{24}\mathrm{)}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\zeta }}_{\mathrm{i}}^{\mathrm{2}}\mathit{/}\mathrm{90}\mathrm{)}$.

5.3. Electron-ion correlations

Using the linear response theory in the Thomas-Fermi limit, Fushiki et al. (1989) evaluated the electron polarization energy for a dense plasma in a strongly quantizing magnetic field at zero temperature, assuming that the ions remain classical (unaffected by the field). A comparison with the analogous zero-field result shows that the strongly quantizing magnetic field (${\mathit{\gamma }}_{\mathrm{m}}{\mathit{r}}_{\mathrm{s}}^{\mathrm{2}}\mathit{>}\mathrm{2.23}$) increases the polarization energy at high densities (r_s ≪ 1) by a factor of $\mathrm{0.8846}\hspace{0.17em}{\mathit{\gamma }}_{\mathrm{m}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{s}}^{\mathrm{4}}$ (Potekhin et al. 1999).

Recently, Sharma & Reddy (2011) calculated the screening of the ion-ion potential due to electrons in a large magnetic field B at T = 0, using the one-loop representation of the polarization function. Their results for the strongly quantizing magnetic field show that the screening is anisotropic, and the screened ion potential exhibits Friedel oscillations with period πħ/p_F in a cylinder of a radius  ~πħ/p_F along the magnetic field line that passes through the Coulomb center. Sharma & Reddy suggest that this long-range oscillatory behavior can affect the ion lattice structure. However, finite temperature should damp these oscillations, so that they are pronounced only at T ≪ T_F, i.e., deep within the triangular domains formed by the lines T_F and ρ_B in Fig. 1. At the typical pulsar magnetic fields B ~ 10¹² G, this requires an unusually low temperature of the neutron-star crust. On the other hand, the conditions T ≪ T_F and ρ < ρ_B can be easily fulfilled in the outer crust of magnetars at B ~ 10¹⁵ G (cf. Fig. 1 and the top panel of Fig. 8), but in this case the Friedel oscillations are strongly suppressed because the electrons are ultrarelativistic.

To the best of our knowledge, the magnetic effects on the electron polarization energy have not been calculated at finite temperatures or in the case where the field is not strongly quantizing. In view of the limited scope and limited applicability of the available results on the magnetic effects, we use the nonmagnetic expression for F_ie in our code (Appendix C).

6. Harmonic Coulomb crystals in the magnetic field

The magnetic effects on Coulomb crystals have been studied only in the harmonic approximation. Nagai & Fukuyama (1982, 1983) calculated phonon spectra of body-centered cubic (bcc), face-centered cubic (fcc), and hexagonal closely-packed (hcp) OCP lattices. They compared the energies of zero-point vibrations at different values of parameters β and R_S and found conditions of stability of every lattice type. However, Baiko (2000, 2009) noticed that their choice of the magnetic-field direction did not provide the minimum of the total energy.

Usov et al. (1980) obtained the equations for oscillation modes of a harmonic OCP crystal and studied its phonon spectrum in a quantizing magnetic field in several limiting cases. These authors discovered a “soft” phonon mode with dispersion relation ω_kα ∝ k² near the center of the Brillouin zone, which leads to the unusual dependence of the heat capacity of the lattice C_V,lat ∝ T^3/2 at T → 0 instead of the Debye law C_V,lat ∝ T³. Usov et al. (1980) argued that a strong magnetic field should increase stability of the crystal.

Baiko (2000, 2009) studied the magnetic effects on the phonon spectrum of the harmonic Coulomb crystals and calculated its energy, entropy, and heat capacity. We have found that his results can be approximately reproduced by the analytical expressions presented below.

6.1. Thermal phonon contributions

Without a magnetic field, the thermal phonon contribution f_th to the reduced free energy of a Coulomb crystal F/N_ik_BT is a function of a single argument η, described by a simple analytical expression (Baiko et al. 2001). The magnetic field introduces the second independent dimensionless argument β. The functional dependence of thermodynamic functions on η and β is not simple. Baiko (2009) identified five characteristic sectors of the η – β plane:

• 1.

  η < 1 and β < η^-1 – weakly magnetized classical crystal,

• 2.

  η > 1 and β < η^-1 – weakly magnetized quantum crystal,

• 3.

  η < 1 and β > η^-1 – strongly magnetized classical crystal,

• 4.

  η > β > η^-1 – strongly magnetized quantum crystal,

• 5.

  β > η > 1 – very strongly magnetized quantum crystal.

Note that the condition β > η^-1 is equivalent to ζ_i > 1. Thus, the magnetic field strongly affects the thermodynamic functions of a Coulomb crystal when ħω_ci > k_BT, that is the same condition as for the gas and liquid phases.

For astrophysical applications, we have constructed an analytical representation of the EOS of the magnetized Coulomb crystal, which is asymptotically exact in each of the five sectors far from their boundaries, exactly recovers the nonmagnetic fit of Baiko et al. (2001) in the limit β → 0, and reaches a reasonable compromise between simplicity and accuracy.

The term f_th in Eq. (32) can be rewritten as f_th = u_th − s_th, where u_th = U_th/N_ik_BT and s_th = S_th/N_ik_B are the thermal contributions to the reduced internal energy and entropy. We approximately represent u_th by the function $\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}\mathrm{=}\frac{{\mathit{u}}_{\mathrm{th}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}{\mathrm{1}\mathrm{+}\mathrm{0.5}{\left[\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{3}\mathit{/}{\mathit{\zeta }}_{\mathrm{i}}{\mathrm{)}}^{\mathrm{5}\mathit{/}\mathrm{2}}\right]}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{5}}}\mathrm{+}\frac{\mathit{\psi }\hspace{0.17em}\mathit{/}\sqrt{\mathrm{1}\mathrm{+}\mathrm{24}\mathit{/}{\mathit{\eta }}^{\mathrm{2}}}}{\mathrm{1}\mathrm{+}\mathrm{15}\hspace{0.17em}\mathit{\beta }\mathit{/}\mathit{\eta }\mathrm{+}\mathit{\psi }}\mathit{,}$(76)where ψ = 12.5   (β/η)^3/2 + 119   (β/η)² and ζ_i = βη, and we represent s_th by the function $\begin{array}{}\mathrm{˜}\\ \mathit{s}\end{array}\mathrm{=}{\mathit{s}}_{\mathrm{th}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}\ln \left(\mathrm{1}\mathrm{+}\frac{\mathrm{20.8}\hspace{0.17em}\mathrm{(}\mathit{\beta }\mathit{/}\mathit{\eta }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\mathrm{122.5}\hspace{0.17em}\mathrm{(}\mathit{\beta }\mathit{/}\mathit{\eta }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{7.9}\hspace{0.17em}\mathit{\beta }\mathit{/}\mathit{\eta }\mathrm{)}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{42}\mathit{/}{\mathit{\eta }}^{\mathrm{2}}\mathrm{)}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\zeta }}_{\mathrm{i}}^{-1}{\mathrm{)}}^{\mathrm{4}}}\right)\mathrm{·}$(77)In these equations, ${\mathit{u}}_{\mathrm{th}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ and ${\mathit{s}}_{\mathrm{th}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ are the values of u_th and s_th at β = 0. Equations (76) and (77) exactly reproduce the known asymptotic limits: u_th = 3 in the classical nonmagnetic limit (η ≪ β ≪ 1), u_th = 2 in the classical magnetic limit (η ≪ β^-1 ≪ 1), u_th = 1 in the case where β ≫ η ≫ 1, and u_th = 0.6s_th ∝ η^−3/2, if η → ∞ at β = constant.

The functions $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}\mathrm{\left(}\mathit{\eta ,\beta }\mathrm{\right)}$ and $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{s}\end{array}\mathrm{\left(}\mathit{\eta ,\beta }\mathrm{\right)}$ are displayed in Figs. 2 and 3. Their accuracy is seen from a comparison with the numerical results (Baiko 2009), also shown in the figures. However, if the complete consistency of different thermodynamic functions is required, Eqs. (76) and (77) should not be used directly, but should be first combined into ${\mathit{f}}_{\mathrm{th}}\mathrm{=}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}\mathrm{-}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{s}\end{array}\mathit{.}$ Then one can calculate thermodynamic functions by differentiating the function f_th(η,β). In this way we obtain, for example, $\begin{array}{ccc}& & {\mathit{u}}_{\mathrm{th}}\mathrm{=}\frac{{\mathit{U}}_{\mathrm{th}}}{{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\mathrm{=}\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}\mathrm{+}\mathrm{\Delta }\mathit{u,}\\ & & {\mathit{s}}_{\mathrm{th}}\mathrm{=}{\mathit{S}}_{\mathrm{th}}\mathit{/}{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathrm{=}\begin{array}{}\mathrm{˜}\\ \mathit{s}\end{array}\mathrm{+}\mathrm{\Delta }\mathit{u,}\\ & & \\ & & \mathrm{where} \mathrm{\Delta }\mathit{u}\mathrm{=}{\left(\frac{\mathit{\partial }\begin{array}{}\mathrm{˜}\\ \mathit{s}\end{array}}{\mathit{\partial }\ln \mathit{\eta }}\right)}_{\mathit{\beta }}\mathrm{-}{\left(\frac{\mathit{\partial }\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}{\mathit{\partial }\ln \mathit{\eta }}\right)}_{\mathit{\beta }}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}\mathit{.}\end{array}$Note that the relation between the phonon contributions to pressure and internal energy, 2P_thV = U_th, which is standard for a nonmagnetized harmonic crystal, is invalid in the strongly magnetized crystal because both dimensionless arguments η and β depend on density.

Approximations in Eqs. (78)–(80) are shown in Figs. 2–4. Their reasonable behavior beyond the range of available numerical data is demonstrated by plotting them also at larger β = 10³ and 10⁴.

6.2. Zero-point vibrations

Because motions of the ions are confined by the magnetic field in the transverse direction, they exhibit quantum oscillations in the ground state (Landau 1930; Landau & Lifshitz 1977). The energy of these oscillations is ħω_ci/2 for every ion, which gives the term ζ_i/2 in Eq. (42). In a Coulomb crystal, the motion of an ion is confined in an effective potential well, centered at its equilibrium lattice site. The total energy of the zero-point quantum lattice oscillations U_q is given by Eq. (30).

In the case where the crystal is placed in a magnetic field, U_q includes contributions due to both magnetic and lattice confinements of the ion motion. However, since the magnetic contribution N_iħω_ci/2 is common in all phase states, we take it as the zero energy point in our code and separate it from the lattice contribution that is specific to the solid phase. Then Eq. (30) becomes ${\mathit{U}}_{\mathrm{q}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{N}}_{\mathrm{i}}ħ{\mathit{\omega }}_{\mathrm{p}}{\mathit{u}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{N}}_{\mathrm{i}}ħ{\mathit{\omega }}_{\mathrm{ci}}\mathit{,}$(82)and in Eq. (32) we have $\mathrm{1.5}{\mathit{u}}_{\mathrm{1}}\mathit{\eta }\mathrm{=}\mathrm{1.5}{\mathit{u}}_{\mathrm{1}}^{\mathrm{\prime }}\mathit{\eta }\mathrm{+}{\mathit{\zeta }}_{\mathrm{i}}\mathit{/}\mathrm{2}$, where we have defined ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{=}{\mathit{u}}_{\mathrm{1}}\mathrm{-}\mathit{\beta }\mathit{/}\mathrm{3.}$(83)The reduced frequency moment ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\prime }}$ still depends on β, because the character of ion vibrations is affected by the magnetic field (they become essentially one-dimensional if B is extremely large), but the latter dependence is relatively weak. Having extracted ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\prime }}$ from the available numerical results for u₁ (Baiko 2009; Baiko & Yakovlev 2012), we can represent it by the simple interpolation ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}\frac{{\mathit{u}}_{\mathrm{1}}^{\mathrm{0}}\mathrm{+}\mathrm{1.27}{\mathit{\beta }}^{\mathrm{9}\mathit{/}\mathrm{8}}{\mathit{u}}_{\mathrm{1}}^{\mathrm{\infty }}}{\mathrm{1}\mathrm{+}\mathrm{1.27}{\mathit{\beta }}^{\mathrm{9}\mathit{/}\mathrm{8}}}\mathit{,}$(84)where ${\mathit{u}}_{\mathrm{1}}^{\mathrm{0}}$ is the zero-field value and ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\infty }}$ is the limit of ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}$ at β → ∞. Only one of the three phonon branches contributes to ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\prime }}$ in the latter limit, therefore ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\infty }}\mathrm{~}{\mathit{u}}_{\mathrm{1}}^{\mathrm{0}}\mathit{/}\mathrm{3}$. For the bcc crystal, ${\mathit{u}}_{\mathrm{1}}^{\mathrm{0}}\mathrm{=}\mathrm{0.5113875}$, whereas ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\infty }}$ varies between 0.18 and 0.19 depending on the orientation of the lattice in the magnetic field. In Fig. 5 we show ${\mathit{u}}_{\mathrm{1}}^{\mathrm{\prime }}$ and the logarithm (base 10) of u₁ versus β.

Usov et al. (1980) noticed that the energy of a crystal depends on its orientation in a strong magnetic field. However, numerical calculations (Baiko 2000, 2009) show that this dependence is very weak. For example, the difference Δu₁ between the values of u₁ for two orientations, where the field lines connect an ion with its nearest neighbor in the first case and with a next-order nearest neighbor in the second case, is approximately $\mathrm{\Delta }{\mathit{u}}_{\mathrm{1}}\mathrm{=}\frac{\hspace{0.17em}{\mathit{\beta }}^{\mathrm{15}\mathit{/}\mathrm{4}}}{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\beta }}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{)}}^{\mathrm{5}\mathit{/}\mathrm{2}}}\hspace{0.17em}\mathrm{(}\mathrm{\Delta }{\mathit{u}}_{\mathrm{1}}{\mathrm{)}}_{\max }$(85)with saturation level (Δu₁)_max = 0.0064 for the bcc lattice.

6.3. Comparison with the Baiko-Yakovlev fit

After the present work was completed, we became aware of an independent study by Baiko & Yakovlev (2012, priv. comm.), who developed another set of approximations for the free energy of the harmonic Coulomb crystal in a magnetic field. They presented the free energy as a sum of three terms corresponding to the contributions from each of the three phonon modes in the bcc crystal. Thus each of these terms has a clear physics meaning, while our fitting expressions give only the total contribution, which cannot be easily decomposed to three parts corresponding to the separate phonon modes.

Unlike our fit, the fit of Baiko & Yakovlev (2012) does not exactly reproduce the very accurate results of Baiko et al. (2001) in the limit β → 0. At finite β, both sets of fitting expressions accurately reproduce the asymptotes at T ≪ T_p and T ≫ T_p and have similar accuracies within several percent points in the intermediate range 0.1T_p/β ≲ T ≲ 10T_p. Meanwhile, our approximation is simpler: the Baiko-Yakovlev approximation contains 27 independent numerical fitting parameters, whereas our fits (76) and (77) contain together only 9 such parameters.

7. Examples and discussion

7.1. Thermodynamic functions

Characteristic features of the EOS can be seen in Fig. 6. Here, we have chosen the plasma parameters that are typical for outer envelopes of isolated neutron stars: we consider fully-ionized iron (Z = 26, A = 56) at T = 10⁷ K and B = 10¹² G (for illustration, the density range is extended to ρ ≲ 10⁵ neglecting the bound states that can be important in this ρ – T domain). We plot the normalized pressure p = P/n_ik_BT, entropy S/N_ik_B, heat capacity c_V = C_V/N_ik_B, and logarithmic derivatives of pressure χ_ρ and χ_T as functions of density. Dashed lines show these functions in the absence of quantizing magnetic field. The vertical dotted lines marked by numbers separate different characteristic domains, consecutively entered with increasing density: onset of electron degeneracy at B = 0 (${\mathit{T}}_{\mathrm{F}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{=}\mathit{T}$) and at B = 10¹² G (T_F = T), population of excited Landau levels (ρ = ρ_B), melting point with formation of a classical Coulomb crystal (T_m = T), and quantum effects in the crystal (T_p = T).

At low densities, the ideal-gas values are approached: p = 1 + Z, χ_ρ = χ_T = 1, c_V = (3 + 3Z)/2 at B = 0, and c_V = (3 + Z)/2 at B = 10¹² G. The latter difference is because at ρ < ρ_B the electron gas is effectively one-dimensional due to the strong magnetic quantization.

With increasing density, the reduced pressure p first decreases below its ideal-gas value due to the Coulomb nonideality and then increases due to the electron degeneracy. The increase occurs earlier at zero field than in the strong magnetic field, because of the delayed onset of the degeneracy (Sect. 4). When ρ > ρ_B ≈ 1.5 × 10⁴ g cm^-3, the thermodynamic functions approach their zero-field values. The gradually decreasing oscillations correspond to consecutive filling of the electron Landau levels. The magnetic field B = 10¹² G does not affect the ion contributions in Fig. 6, because it is nonquantizing for the iron nuclei at T = 10⁷ K (ζ_i = 0.00342).

The liquid-solid phase transition occurs in Fig. 6 at ρ ≈ 8.25 × 10⁴ g cm^-3, where we adopt the classical OCP melting condition Γ = 175.2 (Paper I). With further increase in density (ρ ≳ 10⁶) the degeneracy becomes so strong that the energy and pressure are nearly independent of T and χ_T strongly decreases. The normalized heat capacity gradually tends to its value c_V = 3 characteristic of the classical simple crystal. At still higher density the ion motions become quantized (T_p ≫ T) which leads to the further decrease in the heat capacity and the entropy.

7.2. Melting

The electron polarization, ion quantum effects, and quantizing magnetic field can shift the melting temperature. The lower panel of Fig. 7 shows the Coulomb coupling parameter Γ at the melting (that is, the value Γ_m at which the free energies of the two phases are equal to each other); the upper panel displays the difference between the internal energies in the liquid and solid phases at the melting point (the latent heat Q_m = U_liq − U_sol), divided by N_ik_BT. We plot the data for fully ionized ¹²C and ⁵⁶Fe at B = 0 and 10¹³ G for carbon, B = 0, 10¹⁴ G, and 10¹⁵ G for iron. The density range shown in the figure is typical for the outer envelopes of a neutron star and is also relevant for white dwarfs.

The position of the melting point is very sensitive to the accuracy of the free energies of the Coulomb liquid and crystal (see, e.g., Paper I). The polarization and quantum corrections to the classical OCP free energy are not known sufficiently well for finding the position of the melting point in the whole interval of densities shown in Fig. 7. The dot-dashed curves in this figure correspond to the domain where T < 0.5   T_p. Here, the perturbation theory for the quantum effects in the liquid phase becomes progressively inaccurate. The dashed curves correspond to the domain where the binding energy of the hydrogen-like ion exceeds 10% of the Fermi energy, Z²   Ry > 0.1   ϵ_F, which corresponds to a_e ≳ 0.6a₀/Z. Here, the perturbation treatment of the electron polarization in the crystal starts to be inaccurate. In addition, the position of the melting point cannot be traced to Γ ≲ 100–120, because the available results for the anharmonic corrections to the free energy of a Coulomb crystal (Appendix B.2) are accurate only at larger Γ. Nevertheless, we can evaluate Γ_m in a certain interval of densities for each type of ions (the solid segments of the curves in the figure).

The values of Q_m in Fig. 7 roughly (within a factor of two) agree with the OCP value ${\mathit{Q}}_{\mathrm{m}}^{\mathrm{OCP}}\mathrm{=}\mathrm{0.746}{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}$ at Γ_m = 175.2 and with the values used in theoretical models of white dwarf cooling (e.g., Hansen 2004 and references therein). Most of the neutron star cooling models currently ignore the release of the latent heat at the crystallization of the neutron star envelopes (e.g., Yakovlev & Pethick 2004, and references therein).

Figure 7 shows that the strong magnetic fields tend to decrease Γ_m and thus stabilize the Coulomb crystal, in qualitative agreement with previous conjectures (Ruderman 1971; Kaplan & Glasser 1972; Usov et al. 1980; Lai 2001). At densities ρ ~ 10⁷–10⁸ g cm^-3 corresponding to the “sensitivity strip” in the neutron-star cooling theory (Yakovlev & Pethick 2004), the stabilization proves to be significant at the magnetar field strengths B = 10¹⁴–10¹⁵ G. The results for ⁵⁶Fe in this B interval, shown in the lower panel of Fig. 7, can be roughly (within 10%) described by the formula Γ_m(B) ≈ Γ_m(0)/(1 + 0.2   β). At the typical pulsar field strengths, B = 10¹²–10¹³ G, the effect is noticeable at lower densities. However, these conclusions remain preliminary in the absence of an evaluation of the magnetic-field effects on the anharmonic corrections. In view of the limited applicability and incompleteness of the evaluation of Γ_m with account of the quantum, polarization, and magnetic effects, in applications we use the classical OCP value Γ_m = 175.2 as the fiducial melting criterion.

7.3. Magnetic condensation

Ruderman (1971) suggested that the strong magnetic field may stabilize molecular chains (polymers) aligned with the magnetic field and eventually turn the surface of a neutron star into the metallic solid state. Later studies have provided support for this conjecture, although the critical temperature T_crit, below which this condensation occurs, remains very uncertain. Condensed surface density ρ_s is usually estimated as ${\mathit{\rho }}_{\mathrm{s}\mathit{,\xi }}\mathrm{=}\mathrm{561}\hspace{0.17em}\mathit{\xi }\hspace{0.17em}\mathit{A}{\mathit{Z}}^{-0.6}\hspace{0.17em}{\mathit{B}}_{\mathrm{12}}^{\mathrm{1.2}} \mathrm{g} \mathrm{cm}-3\mathit{,}$(86)where ξ ~ 1 is an unknown numerical factor, which absorbs the theoretical uncertainty (Lai 2001; Medin & Lai 2006). The value ξ = 1 corresponds to the EOS provided by the ion-sphere model (Salpeter 1961), which is close to the uniform model of Fushiki et al. (1989). For comparison, the results of the zero-temperature Thomas-Fermi model for ⁵⁶Fe at  G (Rögnvaldsson et al. 1993) can be approximated (within 4%) by ρ_s,ξ with $\mathit{\xi }\mathrm{\approx }\mathrm{0.2}\mathrm{+}\mathrm{0.01}\mathit{/}{\mathit{B}}_{\mathrm{12}}^{\mathrm{0.56}}$, whereas the finite-temperature Thomas-Fermi model of Thorolfsson et al. (1998) does not predict magnetic condensation at all. Our EOS for partially ionized hydrogen plasmas in strong magnetic fields (Potekhin et al. 1999; Potekhin & Chabrier 2004) exhibits a phase transition with ${\mathit{T}}_{\mathrm{crit}}\mathrm{\approx }\mathrm{3}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}\hspace{0.17em}{\mathit{B}}_{\mathrm{12}}^{\mathrm{0.39}}$ K and critical density ${\mathit{\rho }}_{\mathrm{crit}}\mathrm{\approx }\mathrm{143}\hspace{0.17em}{\mathit{B}}_{\mathrm{12}}^{\mathrm{1.18}} \mathrm{g} \mathrm{cm}{}^{-3}\mathrm{~}{\mathit{\rho }}_{\mathrm{s}\mathit{,}\mathrm{0.25}}$ at 1 ≲ B₁₂ ≲ 10³. According to another study (Lai & Salpeter 1997; Lai 2001), T_crit for hydrogen is several times smaller.

Medin & Lai (2006) performed density-functional calculations of the cohesive energy Q_s of the condensed phases of H, He, C, and Fe in strong magnetic fields. A comparison with previous density-functional calculations of other authors prompts that Q_s may vary within a factor of two at B₁₂ ≳ 1, depending on the approximations (see Medin & Lai 2006 for references and discussion). In a subsequent study, Medin & Lai (2007) calculated the equilibrium densities of saturated vapors of He, C, and Fe atoms and polymers above the condensed surfaces, and obtained T_crit at several values of B by equating the vapor density to ρ_s. Unlike the previous authors, Medin & Lai (2006, 2007) have taken the electronic band structure of the condensed matter into account self-consistently, but they did not allow for the atomic motion across the magnetic field and mostly neglected the contributions of the excited atomic and molecular states in the gaseous phase. Medin & Lai obtained ρ_s assuming that the linear molecular chains form a rectangular array with sides 2R in the plane perpendicular to B, and that the distance a between the nuclei along B in the condensed matter remains the same as in the gaseous phase, so that ρ_s = m_i/4aR² (Medin 2012, priv. comm.). Using Tables 3–5 of Medin & Lai (2006) for a and R, we can describe their results for the surface density of ¹²C at  and ⁵⁶Fe at by Eq. (86) with $\mathit{\xi }\mathrm{=}\mathrm{0.517}\mathrm{+}\mathrm{0.24}\mathit{/}{\mathit{B}}_{\mathrm{12}}^{\mathrm{1}\mathit{/}\mathrm{5}}\mathrm{\pm }\mathrm{0.011}$ and ξ = 0.55 ± 0.11, respectively.

Medin & Lai (2007) found that the critical temperature is T_crit ≈ 0.08Q_s/k_B. Their numerical results for He, C, and Fe can be roughly (within a factor of 1.5) described as ${\mathit{T}}_{\mathrm{crit}}\mathrm{~}\mathrm{5}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\hspace{0.17em}{\mathit{Z}}^{\mathrm{1}\mathit{/}\mathrm{4}}\hspace{0.17em}{\mathit{B}}_{\mathrm{12}}^{\mathrm{3}\mathit{/}\mathrm{4}}$ K at 1 ≲ B₁₂ ≲ 1000. For comparison, the results of Lai & Salpeter (1997) for H at 10 ≲ B₁₂ ≲ 500 suggest ${\mathit{T}}_{\mathrm{crit}}\mathrm{~}\mathrm{1.6}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\hspace{0.17em}{\mathit{B}}_{\mathrm{12}}^{\mathrm{0.7}}$ K. The discrepancies between different estimates of ρ_s and T_crit reflect the current theoretical uncertainties.

The filled dots on the curves for magnetized carbon in both panels of Fig. 7 correspond to the condensed surface position in the fully-ionized plasma model. In this model ${\mathit{T}}_{\mathrm{crit}}\mathrm{\approx }\mathrm{2.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}\hspace{0.17em}{\mathit{Z}}^{\mathrm{0.9}}\hspace{0.17em}{\mathit{B}}_{\mathrm{12}}^{\mathrm{0.4}}$ K and ρ_crit ≈ ρ_s,0.47. With decreasing temperature T below T_crit, the surface density increases and tends to the limit ρ_s,1 given using Eq. (86) as ρ_s ≈ ρ_s,1/ [1 + 1.1   (T/T_crit)⁵]  (cf. Fig. 1). At smaller densities, there is a thermodynamically unstable region in this model, therefore the curves in Fig. 7 are not continued to the left beyond this point. For the magnetized iron model, the melting curve does not cross the surface because T_m > T_crit. In this case, the open circles mark the density that the condensed surface would have at much smaller temperatures T ~ 10⁶   K ≪ T_crit. The parts of the curves to the left of the open circles cannot be reached in a stationary stellar envelope.

7.4. Thermal structure of a magnetar envelope

The results presented above have a direct application to the calculations of the thermal and mechanical structure of neutron-star envelopes with strong magnetic fields. Figure 8 illustrates the structure of a typical magnetar envelope with the ground-state nuclear composition. For illustration we have assumed that the magnetar has mass 1.4   M_⊙ and radius 12 km, and the considered patch of the stellar surface has effective temperature 10^6.5 K and magnetic field B = 10¹⁵ G inclined at 45°. The top panel shows the thermal structure of the envelope, which has been calculated by numerical solution of the system of heat balance equations, taking the general relativity effects and neutrino emission into account (Potekhin et al. 2007). The middle panel presents the ion charge Z as function of ρ (Rüster et al. 2006). In the bottom panel (analogous to Fig. 6) we plot several reduced thermodynamic functions of ρ and T along the thermal profile (i.e., taking T from the top panel), starting at the condensed solid surface.

The temperature quickly grows at the solid surface and reaches the melting point at the depth z ≈ 7 cm. Thus, at the given conditions, the liquid ocean of a magnetar turns out to be covered by a thin layer of “ice” (solid substance). We treat the solid crust as immobile, but the liquid layer below the “ice” is convective up to the depth z ~ 1 m. We treat the convective heat transport through this layer in the adiabatic approximation (Schwarzschild 1958). The change of the heat-transport mechanism from conduction to convection causes the break of the temperature profile at the melting point. We underline that this treatment is only an approximation. In reality, the superadiabatic growth of temperature can lead to a hydrostatic instability of the shell of “ice” and eventually to its cracking and fragmentation into turning-up “ice floes”. This can result in transient enhancements of the thermal luminosity of magnetars.

The temperature profile flattens with density increase, and the Coulomb plasma freezes again at the interface between the layers of ⁶⁶Ni and ⁸⁶Kr at ρ = 1.5 × 10⁹ g cm^-3 (z = 73.8 m). These phase transitions do not cause any substantial breaks in χ_ρ, χ_T, or C_V/N_ik_B, because the Coulomb plasmas have similar structure factors in the liquid and crystalline phases in the melting region (cf. Baiko et al. 1998).

At the boundaries between layers composed of different chemical elements, the reduced thermodynamic functions do not exhibit substantial discontinuities, except for the abrupt increases in P/n_ik_BT at the interfaces ⁶⁶Ni/⁸⁶Kr (ρ = 1.5 × 10⁹ g cm^-3) and ⁷⁸Ni/¹²⁴Mo (ρ = 1.32 × 10¹¹ g cm^-3), which are caused by the decreases in n_i with the large jumps in A (of course, the non-normalized pressure P is continuous). The specific heat per ion C_V/N_ik_B is almost continuous at these interfaces, which means that heat capacity of unit volume abruptly decreases. The drop in χ_T at the ⁷⁸Ni/¹²⁸Mo interface is due to the same decrease in n_i, which leads to the decrease in the ionic contribution that mostly determines ∂P/∂T at the strong degeneracy.

The oscillations of the reduced thermodynamic functions (most noticeable for χ_ρ and χ_T) correspond to consecutive population of excited Landau levels by degenerate electrons with density increase, analogous to the oscillations in Fig. 6.

The magnetic effects on the nonideal part of the plasma thermodynamic functions have almost no influence on the temperature profile in the magnetar envelope, as illustrated in the upper panel of Fig. 8 where the corresponding solid and dashed lines virtually coincide. For comparison, the dotted line in the upper panel shows the result of a calculation totally neglecting the Coulomb nonideality. In this case, the profile is quite different at low densities, where there is no longer a solid surface. However, even in this case the thermal profile is almost the same at large ρ. This means that the Coulomb nonideality has a minor impact on the relation between the internal and effective temperatures and therefore on the cooling curves (Yakovlev & Pethick 2004), but it can be important for the shape of the thermal spectrum (cf., e.g., Potekhin et al. 2012).

8. Conclusions

We have systematically reviewed analytical approximations for the EOS of fully-ionized electron-ion plasmas in magnetic fields and described several improvements to the previously published approximations, taking nonideality attributable to ion-ion, electron-electron, and electron-ion interactions into account. The presented formulae are applicable in a wide range of plasma parameters, including the domains of nondegenerate and degenerate, nonrelativistic and relativistic electrons, weakly and strongly coupled Coulomb liquids, classical and quantum Coulomb crystals. As an application, we have calculated and discussed the behavior of thermodynamic functions, melting, and latent heat at crystallization of strongly coupled Coulomb plasmas with the parameters appropriate for cooling white dwarfs and envelopes of nonmagnetized and strongly magnetized neutron stars. We have also shown that a typical outer envelope of a magnetar can have a liquid layer beneath the solid surface.

---

Acknowledgments

We are grateful to José Pons for pointing out a bug in a previous version of one of our subroutines, to D. A. Baiko and D. G. Yakovlev for making their results available to us prior to publication and for useful discussions, and to D. G. Yakovlev for valuable remarks on a preliminary version of this paper. The work of A.Y.P. was partially supported by the Ministry of Education and Science of the Russian Federation (Agreement No. 8409, 2012), the Russian Foundation for Basic Research (RFBR grant 11-02-00253-a), and the Russian Leading Scientific Schools program (grant NSh-3769.2010.2).

References

Appendix A: Nonideal part of the free energy of electrons

Tanaka et al. (1985) calculated the interaction energy of the electron fluid at finite T and presented a fitting formula that reproduced their numerical results as well as the results of other authors in various limits. Subsequently the behavior of the fit at T ≪ T_F was improved by Ichimaru et al. (1987A.1). The result reads $\begin{array}{ccc}& & \frac{{\mathit{F}}_{\mathit{ee}}}{{\mathit{k}}_{\mathrm{B}}\mathit{T}{\mathit{N}}_{\mathrm{e}}}\mathrm{=}\mathrm{-}\mathit{g}{\mathrm{\Gamma }}_{\mathrm{e}}\mathrm{-}\frac{\mathrm{2}\mathit{A}}{\mathit{f}}\hspace{0.17em}\sqrt{{\mathrm{\Gamma }}_{\mathrm{e}}}\\ & & \mathrm{+}\frac{\mathrm{2}\mathrm{(}\mathit{dB}\mathrm{+}\mathit{CA}\mathrm{)}}{\mathit{fD}}\left[\arctan \left(\frac{\mathrm{2}\mathit{f}\hspace{0.17em}\sqrt{{\mathrm{\Gamma }}_{\mathrm{e}}}\mathrm{+}\mathit{d}}{\mathit{D}}\right)\mathrm{-}\arctan \left(\frac{\mathit{d}}{\mathit{D}}\right)\right]\\ & & \end{array}$(A.1)where A = b − gd, B = a − g, C = 2 − d²/f, $\mathit{D}\mathrm{=}\sqrt{\mathrm{4}\mathit{f}\mathrm{-}{\mathit{d}}^{\mathrm{2}}}$, and a, b, d, f, and g are the following functions of θ = T/T_F (at B = 0): $\begin{array}{ccc}\mathit{a}& \mathrm{=}& {\left(\frac{\mathrm{9}}{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\tanh \frac{\mathrm{1}}{\mathit{\theta }}\\ & & \mathrm{\times }\frac{\mathrm{0.75}\mathrm{+}\mathrm{3.04363}{\mathit{\theta }}^{\mathrm{2}}\mathrm{-}\mathrm{0.09227}{\mathit{\theta }}^{\mathrm{3}}\mathrm{+}\mathrm{1.7035}{\mathit{\theta }}^{\mathrm{4}}}{\mathrm{1}\mathrm{+}\mathrm{8.31051}{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}\mathrm{5.1105}{\mathit{\theta }}^{\mathrm{4}}}\hspace{0.17em}\mathit{,}\\ \mathit{b}& \mathrm{=}& \frac{\mathrm{0.341308}\mathrm{+}\mathrm{12.0708}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}\mathrm{1.148889}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{4}}}{\mathrm{1}\mathrm{+}\mathrm{10.495346}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}\mathrm{1.326623}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{4}}}\sqrt{\mathit{\theta }}\hspace{0.17em}\tanh \frac{\mathrm{1}}{\sqrt{\mathit{\theta }}}\hspace{0.17em}\mathit{,}\\ \mathit{d}& \mathrm{=}& \frac{\mathrm{0.614925}\mathrm{+}\mathrm{16.996055}{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}\mathrm{1.489056}{\mathit{\theta }}^{\mathrm{4}}}{\mathrm{1}\mathrm{+}\mathrm{10.10935}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}\mathrm{1.22184}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{4}}}\sqrt{\mathit{\theta }}\hspace{0.17em}\tanh \frac{\mathrm{1}}{\sqrt{\mathit{\theta }}}\hspace{0.17em}\mathit{,}\\ \mathit{f}& \mathrm{=}& \frac{\mathrm{0.539409}\mathrm{+}\mathrm{2.522206}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}\mathrm{0.178484}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{4}}}{\mathrm{1}\mathrm{+}\mathrm{2.555501}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{2}}\mathrm{+}\mathrm{0.146319}\hspace{0.17em}{\mathit{\theta }}^{\mathrm{4}}}\hspace{0.17em}\mathit{\theta }\hspace{0.17em}\tanh \frac{\mathrm{1}}{\hspace{0.17em}\mathit{\theta }}\hspace{0.17em}\mathit{,}\\ \mathit{g}& \mathrm{=}& \mathrm{0.872496}\mathrm{+}\mathrm{0.025248}\hspace{0.17em}\exp \mathrm{(}\mathrm{-}\mathrm{1}\mathit{/}\mathit{\theta }\mathrm{)}\mathit{.}\end{array}$The accuracy of Eq. (A.1) is 1%.

In a quantizing magnetic field, we replace the argument θ in these expressions by the quantity θ ^∗  defined by Eq. (74), as explained in Sect. 5.1.

Appendix B: Nonideal part of the free energy of ions in the rigid background

B.1. Coulomb liquid

For the reduced free energy f_ii ≡ F_ii/N_ik_BT of the classic OCP, we have the following analytical formula (Paper I): $\begin{array}{ccc}{\mathit{f}}_{\mathrm{ii}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}& \mathrm{=}& {\mathit{A}}_{\mathrm{1}}[\sqrt{\mathrm{\Gamma }\mathrm{(}{\mathit{A}}_{\mathrm{2}}\mathrm{+}\mathrm{\Gamma }\mathrm{)}}\mathrm{-}{\mathit{A}}_{\mathrm{2}}\ln \left(\sqrt{\mathrm{\Gamma }\mathit{/}{\mathit{A}}_{\mathrm{2}}}\mathrm{+}\sqrt{\mathrm{1}\mathrm{+}\mathrm{\Gamma }\mathit{/}{\mathit{A}}_{\mathrm{2}}}\right)]\\ & & \mathrm{+}\mathrm{2}{\mathit{A}}_{\mathrm{3}}\left[\sqrt{\mathrm{\Gamma }}\mathrm{-}\arctan \sqrt{\mathrm{\Gamma }}\right]\\ & \mathrm{+}& {\mathit{B}}_{\mathrm{1}}\left[\mathrm{\Gamma }\mathrm{-}{\mathit{B}}_{\mathrm{2}}\ln \left(\mathrm{1}\mathrm{+}\frac{\mathrm{\Gamma }}{{\mathit{B}}_{\mathrm{2}}}\right)\right]\mathrm{+}\frac{{\mathit{B}}_{\mathrm{3}}}{\mathrm{2}}\hspace{0.17em}\ln \left(\mathrm{1}\mathrm{+}\frac{{\mathrm{\Gamma }}^{\mathrm{2}}}{{\mathit{B}}_{\mathrm{4}}}\right)\mathit{,}\end{array}$(B.1)where A₁ = −0.907347, A₂ = 0.62849, B₁ = 0.0045, B₂ = 170, B₃ = −8.4 × 10^-5, B₄ = 0.0037, and ${\mathit{A}}_{\mathrm{3}}\mathrm{=}\mathrm{-}\sqrt{\mathrm{3}}\mathit{/}\mathrm{2}\mathrm{-}{\mathit{A}}_{\mathrm{1}}\mathit{/}\sqrt{{\mathit{A}}_{\mathrm{2}}}$. The derivative $\frac{\mathit{\partial }{\mathit{f}}_{\mathrm{ii}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}}{\mathit{\partial }\ln \mathrm{\Gamma }}\mathrm{=}{\mathrm{\Gamma }}^{\mathrm{3}\mathit{/}\mathrm{2}}\left[\frac{{\mathit{A}}_{\mathrm{1}}}{\sqrt{\mathrm{\Gamma }\mathrm{+}{\mathit{A}}_{\mathrm{2}}}}\mathrm{+}\frac{{\mathit{A}}_{\mathrm{3}}}{\mathrm{\Gamma }\mathrm{+}\mathrm{1}}\right]\mathrm{+}\frac{{\mathit{B}}_{\mathrm{1}}\hspace{0.17em}{\mathrm{\Gamma }}^{\mathrm{2}}}{\mathrm{\Gamma }\mathrm{+}{\mathit{B}}_{\mathrm{2}}}\mathrm{+}\frac{{\mathit{B}}_{\mathrm{3}}\hspace{0.17em}{\mathrm{\Gamma }}^{\mathrm{2}}}{{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathrm{4}}}\mathit{,}$(B.2)reproduces Monte Carlo calculations of the reduced internal energy U_ii/N_ik_BT at 1 ≤ Γ ≤ 190 (Caillol 1999) within the accuracy of these calculations,  ≲10^-3. For any values of the coupling parameter in a liquid OCP, , the fractional error of the approximation (B.1) does not exceed 2 × 10^-4.

The classical treatment of ion motion is justified at T ≫ T_p. One can extend the applicability range of the analytical EOS to T ~ T_p by the Wigner-Kirkwood quantum corrections (Wigner 1932; Landau & Lifshitz 1980). The lowest-order correction to the reduced free energy is ${\mathit{f}}_{\mathrm{ii}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}{\mathit{\eta }}^{\mathrm{2}}\mathit{/}\mathrm{24.}$(B.3)The next-order correction  ∝ ħ⁴ was obtained by Hansen & Vieillefosse (1975). It can be written as $\begin{array}{ccc}{\mathit{f}}_{\mathrm{ii}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}& \mathrm{\approx }& (\mathrm{-}\mathrm{2.085}\mathrm{\times }{\mathrm{10}}^{-4}\mathrm{-}\frac{\mathrm{2.411}\mathrm{\times }{\mathrm{10}}^{-4}}{{\mathrm{\Gamma }}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{-}\frac{\mathrm{0.001288}}{\mathrm{\Gamma }}\\ & & \end{array}$(B.4)This expression, unlike Eq. (B.3), is not exact. Both corrections have limited applicability, because as soon as η becomes large, the Wigner expansion diverges and the plasma forms a quantum liquid, whose free energy is not known in an analytical form. Therefore we use only the lowest-order correction (B.3), i.e., ${\mathit{f}}_{\mathrm{ii}}\mathrm{\approx }{\mathit{f}}_{\mathrm{ii}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{+}{\mathit{f}}_{\mathrm{ii}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$. In a magnetic field, Eq. (B.3) is replaced by Eq. (75) (Sect. 5.2).

B.2. Coulomb crystal

The reduced free energy of an OCP in the crystalline phase is given by Eq. (32), where the first three terms describe the harmonic lattice model (Baiko et al. 2001). For the bcc crystal, we have C₀ = −0.895   929   255   68 and u₁ = 0.511   3875, and for f_th the following fitting formula can be used: ${\mathit{f}}_{\mathrm{th}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{3}}\ln \left(\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\alpha }}_{\mathit{k}}\mathit{\eta }}\right)\mathrm{-}\frac{\mathit{A}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}}{\mathit{B}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}}\mathit{,}$(B.5)where α₁ = 0.932446, α₂ = 0.334547, α₃ = 0.265764, $\begin{array}{ccc}\mathit{A}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}& \mathrm{=}& \mathrm{1}\mathrm{+}\mathrm{0.1839}\hspace{0.17em}\mathit{\eta }\mathrm{+}\mathrm{0.593}\hspace{0.17em}\mathrm{586}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{2}}\mathrm{+}\mathrm{0.005}\hspace{0.17em}\mathrm{4814}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{3}}\\ & & \mathrm{+}\mathrm{5.018}\hspace{0.17em}\mathrm{13}\mathrm{\times }{\mathrm{10}}^{-4}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{4}}\mathrm{+}\mathrm{3.9247}\mathrm{\times }{\mathrm{10}}^{-7}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{6}}\\ & & \mathrm{+}\mathrm{5.8356}\mathrm{\times }{\mathrm{10}}^{-11}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{8}}\mathit{,}\\ \mathit{B}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}& \mathrm{=}& \mathrm{261.66}\mathrm{+}\mathrm{7.079}\hspace{0.17em}\mathrm{97}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{2}}\mathrm{+}\mathrm{0.0409}\hspace{0.17em}\mathrm{484}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{4}}\\ & & \mathrm{+}\mathrm{3.973}\hspace{0.17em}\mathrm{55}\mathrm{\times }{\mathrm{10}}^{-4}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{5}}\mathrm{+}\mathrm{5.111}\hspace{0.17em}\mathrm{48}\mathrm{\times }{\mathrm{10}}^{-5}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{6}}\\ & & \mathrm{+}\mathrm{2.197}\hspace{0.17em}\mathrm{49}\mathrm{\times }{\mathrm{10}}^{-6}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{7}}\mathrm{+}\mathrm{1.866}\hspace{0.17em}\mathrm{985}\mathrm{\times }{\mathrm{10}}^{-9}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{9}}\\ & & \mathrm{+}\mathrm{2.787}\hspace{0.17em}\mathrm{72}\mathrm{\times }{\mathrm{10}}^{-13}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{11}}\mathit{.}\end{array}$The Taylor expansion of Eq. (B.5) at small η is consistent with the Wigner correction (B.3). However, the next Taylor term  ~η³ is absent in the Wigner expansion, and therefore Eq. (B.5) does not reproduce higher-order Wigner corrections. Nevertheless, approximation (B.5) is very accurate: it reproduces the numerical results in Baiko et al. (2001) with fractional deviations within 5 × 10^-6, and its first and second derivatives reproduce the calculated contributions to the internal energy and heat capacity with deviations up to several parts in 10⁵. Other types of simple lattices are described by the same expressions with slightly different parameters (see Baiko et al. 2001).

Anharmonic corrections for Coulomb lattices were studied in a number of works (see Papers I and II for references). In the classical regime η → 0, we have chosen one of the 11 parametrizations proposed by Farouki & Hamaguchi (1993): ${\mathit{f}}_{\mathrm{ah}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathrm{\Gamma }\mathrm{\right)}\mathrm{=}\mathrm{-}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathrm{3}}\frac{{\mathit{a}}_{\mathit{k}}}{\mathit{k}{\mathrm{\Gamma }}^{\mathit{k}}}\mathit{,}$(B.6)where a₁ = 10.9, a₂ = 247, and a₃ = 1.765 × 10⁵. A continuation to arbitrary η, which is consistent with available analytical and numerical results for quantum crystals, reads (Paper II) ${\mathit{f}}_{\mathrm{ah}}\mathrm{=}{\mathit{f}}_{\mathrm{ah}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathrm{\Gamma }\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{0.0112}{\mathit{\eta }}^{\mathrm{2}}}\mathrm{-}\mathrm{0.12}\hspace{0.17em}{\mathit{\eta }}^{\mathrm{2}}\mathit{/}\mathrm{\Gamma }\mathit{.}$(B.7)Superstrong magnetic fields can significantly change these expressions under the conditions ζ_i ≳ 1 and η ≳ 1. Analytical approximations for the free energy of a harmonic Coulomb crystal in quantizing magnetic fields are derived in Sect. 6. Analogous results for the anharmonic corrections are currently unavailable.

Appendix C: Electron polarization corrections

C.1. Coulomb liquid

The screening contribution to the reduced free energy of the Coulomb liquid at 0 < Γ ≲ 300 has been calculated by the HNC technique and fitted by the expression (Paper I) ${\mathit{f}}_{\mathrm{ie}}\mathrm{\equiv }\frac{{\mathit{F}}_{\mathrm{ie}}}{{\mathit{N}}_{\mathrm{i}}{\mathit{k}}_{\mathrm{B}}\mathit{T}}\mathrm{=}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{e}}\hspace{0.17em}\frac{{\mathit{c}}_{\mathrm{DH}}\sqrt{{\mathrm{\Gamma }}_{\mathrm{e}}}\mathrm{+}{\mathit{c}}_{\mathrm{TF}}\mathit{a}{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathit{\nu }}{\mathit{g}}_{\mathrm{1}}{\mathit{h}}_{\mathrm{1}}}{\mathrm{1}\mathrm{+}\mathrm{[}\mathit{b}\hspace{0.17em}\sqrt{{\mathrm{\Gamma }}_{\mathrm{e}}}\mathrm{+}\mathit{a}{\mathit{g}}_{\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathit{\nu }}\mathit{/}{\mathit{r}}_{\mathrm{s}}\mathrm{]}{\mathit{\gamma }}_{\mathrm{r}}^{-1}}\mathit{,}$(C.1)where ${\mathit{c}}_{\mathrm{DH}}\mathrm{=}\mathrm{(}\mathit{Z}\mathit{/}\sqrt{\mathrm{3}}\mathrm{)}{}^{\mathrm{[}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{Z}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{-}{\mathit{Z}}^{\mathrm{3}\mathit{/}\mathrm{2}}{}^{\mathrm{]}}$ ensures exact transition to the Debye-Hückel limit at Γ → 0, c_TF = (18/175)   (12/π)^2/3Z^7/3(1 − Z ^− 1/3 + 0.2   Z^−1/2( fits the numerical data at large Γ and reproduces the Thomas-Fermi limit (Salpeter 1961) at Z → ∞, the parameters a = 1.11   Z^0.475, b = 0.2 + 0.078   (lnZ)², and ν = 1.16 + 0.08lnZ provide a low-order approximation to F_ie for intermediate r_s and Γ. The functions $\begin{array}{ccc}{\mathit{g}}_{\mathrm{1}}& \mathrm{=}& \mathrm{1}\mathrm{+}\mathrm{0.78}{\mathrm{[}\mathrm{21}\mathrm{+}{\mathrm{\Gamma }}_{\mathrm{e}}\mathrm{(}\mathit{Z}\mathit{/}{\mathit{r}}_{\mathrm{s}}{\mathrm{)}}^{\mathrm{3}}\mathrm{]}}^{-1}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{e}}\mathit{/}\mathit{Z}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}\\ {\mathit{g}}_{\mathrm{2}}& \mathrm{=}& \mathrm{1}\mathrm{+}\frac{\mathit{Z}\mathrm{-}\mathrm{1}}{\mathrm{9}}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{1}}{\mathrm{0.001}\hspace{0.17em}{\mathit{Z}}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathrm{\Gamma }}_{\mathrm{e}}}\right)\frac{{\mathit{r}}_{\mathrm{s}}^{\mathrm{3}}}{\mathrm{1}\mathrm{+}\mathrm{6}\hspace{0.17em}{\mathit{r}}_{\mathrm{s}}^{\mathrm{2}}}\hspace{0.17em}\mathit{,}\end{array}$improve the fit at relatively large r_s. Finally, the function ${\mathit{h}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{1}\mathrm{+}{\mathit{x}}_{\mathrm{r}}^{\mathrm{2}}\mathit{/}\mathrm{5}}{\mathrm{1}\mathrm{+}\mathrm{0.18}\hspace{0.17em}{\mathit{Z}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{4}}{\mathit{x}}_{\mathrm{r}}\mathrm{+}\mathrm{0.37}\hspace{0.17em}{\mathit{Z}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{x}}_{\mathrm{r}}^{\mathrm{2}}\mathrm{+}{\mathit{x}}_{\mathrm{r}}^{\mathrm{2}}\mathit{/}\mathrm{5}}$is the relativistic correction, as is ${\mathit{\gamma }}_{\mathrm{r}}^{-1}$ in the denominator.

C.2. Coulomb crystal

The screening contribution to the reduced free energy of the Coulomb crystals was evaluated using the semiclassical perturbation approach with an effective structure factor (Paper I) and fitted by the expression (Paper II) ${\mathit{f}}_{\mathrm{ie}}\mathrm{=}\mathrm{-}{\mathit{f}}_{\mathrm{\infty }}\mathrm{\left(}{\mathit{x}}_{\mathrm{r}}\mathrm{\right)}\hspace{0.17em}\mathrm{\Gamma }\mathrm{\{}\mathrm{1}\mathrm{+}\mathit{A}\mathrm{(}{\mathit{x}}_{\mathrm{r}}\mathrm{)}\hspace{0.17em}[\mathit{Q}\mathrm{(}\mathit{\eta }\mathrm{)}\mathit{/}\mathrm{\Gamma }{]}^{\mathit{s}}\mathrm{\}}\mathit{,}$(C.2)where $\begin{array}{ccc}{\mathit{f}}_{\mathrm{\infty }}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& {\mathit{a}}_{\mathrm{TF}}{\mathit{Z}}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{b}}_{\mathrm{1}}\hspace{0.17em}\sqrt{\mathrm{1}\mathrm{+}{\mathit{b}}_{\mathrm{2}}\mathit{/}{\mathit{x}}^{\mathrm{2}}}\mathit{,}\\ \mathit{A}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{b}}_{\mathrm{3}}\mathrm{+}\mathrm{17.9}{\mathit{x}}^{\mathrm{2}}}{\mathrm{1}\mathrm{+}{\mathit{b}}_{\mathrm{4}}{\mathit{x}}^{\mathrm{2}}}\mathit{,}\\ \mathit{Q}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}& \mathrm{=}& {\left[\ln \left(\mathrm{1}\mathrm{+}{\mathit{e}}^{\mathrm{(}\mathrm{0.205}\mathit{\eta }{\mathrm{)}}^{\mathrm{2}}}\right)\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\left[\ln \left(\mathit{e}\mathrm{-}\mathrm{(}\mathit{e}\mathrm{-}\mathrm{2}\mathrm{)}{\mathit{e}}^{\mathrm{-}\mathrm{(}\mathrm{0.205}\mathit{\eta }{\mathrm{)}}^{\mathrm{2}}}\right)\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}\end{array}$the parameter a_TF = 0.00352 is related to c_TF in Eq. (C.1), and parameters s and b₁ – b₄ depend on Z: $\begin{array}{ccc}\mathit{s}& \mathrm{=}& {\mathrm{[}\mathrm{1}\mathrm{+}\mathrm{0.01}\hspace{0.17em}\mathrm{(}\ln \mathit{Z}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{+}\mathrm{0.097}\hspace{0.17em}{\mathit{Z}}^{-2}\mathrm{]}}^{-1}\mathit{,}\\ {\mathit{b}}_{\mathrm{1}}& \mathrm{=}& \mathrm{1}\mathrm{-}\mathrm{1.1866}\hspace{0.17em}{\mathit{Z}}^{-0.267}\mathrm{+}\mathrm{0.27}\hspace{0.17em}{\mathit{Z}}^{-1}\mathit{,}\\ {\mathit{b}}_{\mathrm{2}}& \mathrm{=}& \mathrm{1}\mathrm{+}\frac{\mathrm{2.25}}{{\mathit{Z}}^{\mathrm{1}\mathit{/}\mathrm{3}}}\hspace{0.17em}\frac{\mathrm{1}\mathrm{+}\mathrm{0.684}\hspace{0.17em}{\mathit{Z}}^{\mathrm{5}}\mathrm{+}\mathrm{0.222}\hspace{0.17em}{\mathit{Z}}^{\mathrm{6}}}{\mathrm{1}\mathrm{+}\mathrm{0.222}\hspace{0.17em}{\mathit{Z}}^{\mathrm{6}}}\mathit{,}\\ {\mathit{b}}_{\mathrm{3}}& \mathrm{=}& \mathrm{41.5}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\ln \mathit{Z}\mathrm{)}\mathit{,}\\ {\mathit{b}}_{\mathrm{4}}& \mathrm{=}& \mathrm{0.395}\ln \mathit{Z}\mathrm{+}\mathrm{0.347}\hspace{0.17em}{\mathit{Z}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathit{.}\end{array}$Here, the numerical parameters are given for the bcc crystal; their values for the fcc lattice are slightly different (Paper I).

All Figures