Equation of state for magnetized Coulomb plasmas^{⋆}
^{1}
CRAL (UMR CNRS 5574), École Normale Supérieure de Lyon,
69364
Lyon Cedex 07,
France
email: palex@astro.ioffe.ru, chabrier@enslyon.fr
^{2}
Ioffe PhysicalTechnical Institute, Politekhnicheskaya 26, 194021
St. Petersburg,
Russia
^{3}
Isaac Newton Institute of Chile, St. Petersburg Branch,
Russia
^{4}
School of Physics, University of Exeter,
Exeter,
EX4 4QL,
UK
Received:
23
July
2012
Accepted:
5
December
2012
We have developed an analytical equation of state (EOS) for magnetized fullyionized plasmas that cover a wide range of temperatures and densities, from lowdensity classical plasmas to relativistic, quantum plasma conditions. This EOS directly applies to calculations of structure and evolution of strongly magnetized white dwarfs and neutron stars. We review available analytical and numerical results for thermodynamic functions of the nonmagnetized and magnetized Coulomb gases, liquids, and solids. We propose a new analytical expression for the free energy of solid Coulomb mixtures. Based on recent numerical results, we have constructed analytical approximations for the thermodynamic functions of harmonic Coulomb crystals in quantizing magnetic fields. The analytical description ensures a consistent evaluation of all astrophysically important thermodynamic functions based on the first, second, and mixed derivatives of the free energy. Our numerical code for calculation of thermodynamic functions based on these approximations has been made publicly available. Using this code, we calculate and discuss the effects of electron screening and magnetic quantization on the position of the melting point in a range of densities and magnetic fields relevant to white dwarfs and outer envelopes of neutron stars. We consider also the thermal and mechanical structure of a magnetar envelope and argue that it can have a frozen surface which covers the liquid ocean above the solid crust.
Key words: dense matter / equation of state / magnetic fields / stars: neutron / white dwarfs
The Fortran code that realizes the analytical approximations described in this paper is available at http://www.ioffe.ru/astro/EIP/ and at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/550/A43
© ESO, 2013
1. Introduction
Coulomb plasmas, i.e., the fullyionized 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 lowmass 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 widerange 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^{7}−10^{9} G (e.g., Wickramasinghe & Ferrario 2000, and references therein) and neutron stars with typical B ~ 10^{8}−10^{14} 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 neutronstar 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 freeelectron 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_{0}, where a_{0} is the Bohr radius and . The parameter r_{s} can be quickly evaluated from the relations where n_{24} ≡ n_{e}/10^{24} cm^{3} and ρ_{0} = 2.6752 (A/Z) g cm^{3}. The analogous density parameter for the ions is R_{S} = a_{i}m_{i}(Ze)^{2}/ħ^{2} = 1822.89AZ^{7/3} r_{s}, where m_{i} is the ion mass and is the ion sphere radius.
At stellar densities it is convenient to use, instead of r_{s}, the nonmagnetic relativity parameter (1)where p_{F} = ħ (3π^{2}n_{e})^{1/3} is the electron Fermi momentum in the absence of a magnetic field, and ρ_{6} ≡ ρ/10^{6} g cm^{3}. The Fermi energy (without the rest energy) for the electron gas is and the Fermi temperature T_{F} ≡ ϵ_{F}/k_{B} = T_{r} (γ_{r} − 1), where T_{r} ≡ m_{e}c^{2}/k_{B} = 5.93 × 10^{9} K, , and k_{B} is the Boltzmann constant. A useful measure of electron degeneracy is θ = T/T_{F}. In the nonrelativistic limit (x_{r} ≪ 1), K, and (2)where (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 (4)where T_{6} ≡ T/10^{6} K.
Thermal de Broglie wavelengths of free ions and electrons are usually defined as (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^{2} n_{i}Z^{2}/m_{i}(^{1/2} is the ion plasma frequency. Since , the importance of the quantum effects in strongly coupled plasmas (i.e., at Γ ≫ 1) is determined by parameter (6)
2.2. Magneticfield parameters
In the nonrelativistic theory (Landau & Lifshitz 1977), the energy of an electron in magnetic field B equals , where p_{z} is the momentum component along B, ω_{c} = eB/m_{e}c is the electron cyclotron frequency, characterizes a Landau level, σ = ±1 determines the spin projection on the field, and n_{L} is the nonnegative 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 interrelated as The levels are doubledegenerate with respect to σ. Their splitting due to the anomalous magnetic moment of the electron is ≈ (m_{e}c^{2}α_{f}/2π) b at b ≪ 1 and ~(m_{e}c^{2}α_{f}/2π) [lnb − 1.584] ^{2} 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_{B}T: (9)where B_{0} = 2.3505 × 10^{9} G, (10)where α_{f} = e^{2}/ħc is the finestructure constant, and (11)where B_{12} ≡ B/10^{12} G. The magnetic length 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 (12)Another important parameter is the ratio of the ion cyclotron frequency to the plasma frequency, (13)
2.3. Free energy and thermodynamic functions
The Helmholtz free energy F of a plasma can be conveniently written as (14)where and denote the ideal free energy of the ions and the electrons, and the last three terms represent an excess free energy arising from the electronelectron, ionion, and ionelectron 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 onecomponent 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 in the OCP model. Still stronger simplification is the ionsphere 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 exchangecorrelation term is defined as 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 ionelectron (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 secondorder thermodynamic functions are derived by differentiating these firstorder 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 secondorder 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_{i}V nonrelativistic classical ions is (15)where M_{spin} is the spin multiplicity. Accordingly, , , , and . In the OCP, Eq. (15) can be written in terms of the dimensionless plasma parameters (Sect. 2) as (16)The free energy of the electron gas is given by (17)where N_{e} = n_{e}V is the number of electrons and μ_{e} is the electron chemical potential without the rest energy m_{e}c^{2}. The pressure and the number density are functions of μ_{e} and T: where χ_{e} ≡ μ_{e}/k_{B}T, τ ≡ T/T_{r}, and (20)is the FermiDirac integral. The internal energy is (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 secondorder thermodynamic functions are obtained using relations of the type 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^{4}, 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^{1}(24)where is the electron chemical potential (without the rest energy) in the relativistic units, Here, we have introduced notations and . At strong degeneracy, , , and . In Paper II we also described an alternative expansion in powers of τ, which allows one to avoid numerical cancellations of close terms at small (we switch to this alternative expansion at ).
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^{4} at τ ≲ 10^{2}, 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 whitedwarf 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 electronelectron 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 ionion 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 WignerKirkwood 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 bodycentered cubic (bcc) Coulomb lattice. In the harmonic approximation (e.g., Kittel 1963), the free energy of the lattice is (29)where U_{0} = N_{i}C_{0}(Ze)^{2}/a_{i} is the classical staticlattice energy, C_{0} ≈ −0.9 is the Madelung constant, (30)accounts for zeropoint quantum vibrations, u_{1} = ⟨ ω_{kα} ⟩ _{ph}/ω_{p} ≈ 0.5 is the reduced first moment of phonon frequencies, (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 classicalgas free energy, therefore F_{lat} replaces in Eq. (14).
Beyond the harmoniclattice approximation, the total reduced free energy f_{lat} ≡ F_{lat}/N_{i}k_{B}T can be written as (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_{0} and u_{1} were calculated by Baiko (2000) (see Appendix B.2). For f_{th} = F_{th}/N_{i}k_{B}T, 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_{1}η + η^{2}/24, where the term with u_{1} cancels that in Eq. (32) and the last term represents the WignerKirkwood quantum correction (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 finitetemperature 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 hypernettedchain (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 DebyeHückel limit for the weakly coupled (Γ ≪ 1) electronion plasmas and the ThomasFermi 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. Moleculardynamics and pathintegral 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 smallwavenumber asymptote of the electron dielectric function (Jancovici 1962; Galam & Hansen 1976). The firstorder 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_{i}k_{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), (33)where x_{j} are the number fractions of ions with charge numbers Z_{j} and . In Eq. (33), f_{ex} ≡ F_{ex}/N_{i}k_{B}T 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 chargeneutralizing 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 finitetemperature electron background has been examined by Hansen et al. (1977) in the firstorder 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 interionic 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 strongcoupling 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) (34)where , ⟨ Γ ⟩ = Γ_{e} ⟨ Z^{5/3} ⟩ , δ is defined either as (35)for rigid electron background model, or as (36)for polarizable background, and parameters a, b, α, and β depend on the plasma composition as follows: 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 linearmixing 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_{2}/Z_{1} at a fixed number fraction x_{2} = 1 − x_{1} for a binary ion mixture with Z_{2}/Z_{1} > 2 and x_{2} < 0.5. A much simpler fit, which does not exhibit unphysical behavior, was suggested by DeWitt & Slattery (2003). It can be written as However, the latter fit is valid only for relatively small charge ratios R_{Z} ≲ 3/2. We replace it by the expression (39)where (40)and The approximation in Eq. (40) reproduces reasonably well the results of both Ogata et al. (1993) and DeWitt & Slattery (2003) for random twocomponent ionic bcc lattices. For a multicomponent ion crystal, Medin & Cumming (2010) proposed the extrapolation from the twocomponent plasma case (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) (42)The last term arises from the energy of the magnetic moments of the ions in the magnetic field, (43)where M_{spin} is the ion spin multiplicity, and g_{i} is the gfactor (g_{i} = 5.5857 for protons). For ions with spin onehalf (M_{spin} = 2), the expression in the square brackets in Eq. (43) simplifies to [2cosh(g_{i} ζ_{i}/4)] . For zerospin ions, such as ^{4}He, ^{56}Fe, and other eveneven nuclei in the ground state, F_{spin} = 0.
The ion pressure obeys the nonmagnetic idealgas relation , but expressions for the internal energy and heat capacity are different: Here, the terms u_{spin} and c_{spin} arise from F_{spin}, They simplify at M_{spin} = 2: (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 Bprojections of the spin and the orbital moment and with a fixed Landau number n in a volume V equals (Landau & Lifshitz 1977). Thus one can express the electron number density n_{e} and the grand thermodynamic potential as 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 where and The free energy is given by Eqs. (17), (51), and (52).
The calculation of n_{e}, , and their derivatives at given χ_{e} and τ can be performed using Eqs. (51, (52) and the same analytical approximations to the FermiDirac 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 secondorder 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 FermiDirac 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 (53)where is the zerofield Fermi momentum at the given density. Equation (53) can be written as p_{F} = m_{e}cx_{B}, where (54)is the relativity parameter modified by the field and (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 . Therefore, the field is strongly quantizing at T ≪ T_{cycl} and ρ < ρ_{B}, where T_{cycl} = ħω_{c}/k_{B} ≈ 1.3434 × 10^{8}B_{12} K and (55)The condition n_{e} < n_{B} can be written as . Then Eq. (53) shows that in a strongly quantizing field , except for densities n_{e} close to the threshold n_{B}. Thus T_{F} is reduced, compared to its nonmagnetic value , by factor . In the nonrelativistic limit, , and the parameter θ = T/T_{F} becomes (56)where θ_{0} 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 zerofield 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 neutronstar envelopes at two magnetic field strengths, B = 10^{12} G and 10^{15} G, assuming fullyionized iron (this assumption may be crude in the lower left part of the diagram). In the stronglyquantizing magneticfield 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).
Fig. 1 Characteristic densitytemperature domains at B = 10^{12} G (blue online) and 10^{15} G (red online) for fullyionized iron. Solid lines indicate the Fermi temperature as function of density, the dotted line shows the plasma temperature, the dotdashed line shows the melting temperature as function of density, short and long dashes delimit the domains of strongly quantizing magnetic field and of magnetic condensation, respectively, and the heavy dots mark the critical point for the condensation (Sect. 7.3). 

Open with DEXTER 
4.2.3. Strongly degenerate electrons
If the electrons are strongly degenerate, then one can apply the Sommerfeld expansion (Sect. 3.1) and obtain , where is the zerotemperature value and ΔF is a thermal correction. According to Eqs. (24) and (52), the zerotemperature pressure P_{0} is (57)where is the relativistic unit of pressure, n_{max} is the maximum integer n for which , and . According to Eqs. (24) and (51), the Fermi energy ϵ_{F} is determined by the condition (58)In order to obtain the chemical potential μ_{e} = ϵ_{F} + Δϵ with fractional accuracy ~ , we retain two terms in Eq. (24), insert it into Eq. (51), approximate in the vicinity of by (59)where and , and drop the higherorder terms containing . Then (60)The thermal correction to the pressure equals (61)and the thermal correction to the free energy and internal energy (62)The leading contribution to the heat capacity is . As in the nonmagnetic case, 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 In the nonrelativistic (x_{B} ≪ 1) and ultrarelativistic (x_{B} ≫ 1) limits, we have and , respectively. Compared with the nonmagnetic case (Papers I and II), the dependence 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 The last equation differs from the nonmagnetic case (Paper II) in that x_{B} replaces x_{r} and π^{2}/3 replaces π^{2}.
4.2.5. Nonrelativistic limit
In the nonrelativistic limit (p_{F} ≪ m_{e}c and T ≪ T_{r}), Eqs. (51) and (52) simplify to (69)In the nondegenerate regime (T ≫ T_{F}), one has I_{ν}(χ) ≈ e^{χ} Γ(ν + 1), where Γ(ν + 1) is the gammafunction. Then Eq. (69) yields and (70)which provides the free energy .
As follows from Eq. (70), the reduced internal energy and heat capacity 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 onedimensional. 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 is recovered. In the strongly quantizing, nondegenerate regime (ρ < ρ_{B} and T_{F} ≪ T ≪ T_{cycl}), the last term of Eq. (70) vanishes, which yields (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 , 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 (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 zcomponent of the Lorentz force is absent, and we get the standard equation of hydrostatic equilibrium: . In the second case, the gradients of the kinetic pressure 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: (73)Then the equilibrium condition takes the form 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 fieldfree 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 (74)ξ = [1−(4/ζ)tanh(ζ/4)] ^{1/2}, and θ_{0} and θ_{B} are given by Eqs. (2) and (56), respectively, at fixed n_{e} and T.
5.2. WignerKirkwood 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 WignerKirkwood term. Its expression in an arbitrary magnetic field was obtained by Alastuey & Jancovici (1980): (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}, .
5.3. Electronion correlations
Using the linear response theory in the ThomasFermi 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 zerofield result shows that the strongly quantizing magnetic field () increases the polarization energy at high densities (r_{s} ≪ 1) by a factor of (Potekhin et al. 1999).
Recently, Sharma & Reddy (2011) calculated the screening of the ionion potential due to electrons in a large magnetic field B at T = 0, using the oneloop 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 longrange 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^{12} G, this requires an unusually low temperature of the neutronstar 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^{15} 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 bodycentered cubic (bcc), facecentered cubic (fcc), and hexagonal closelypacked (hcp) OCP lattices. They compared the energies of zeropoint 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 magneticfield 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^{2} 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^{3}. 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_{i}k_{B}T 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.
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_{i}k_{B}T and s_{th} = S_{th}/N_{i}k_{B} are the thermal contributions to the reduced internal energy and entropy. We approximately represent u_{th} by the function (76)where ψ = 12.5 (β/η)^{3/2} + 119 (β/η)^{2} and ζ_{i} = βη, and we represent s_{th} by the function (77)In these equations, and 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 and 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 Then one can calculate thermodynamic functions by differentiating the function f_{th}(η,β). In this way we obtain, for example, Note that the relation between the phonon contributions to pressure and internal energy, 2P_{th}V = 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^{3} and 10^{4}.
Fig. 2 Thermal phonon contribution to the reduced internal energy u_{th} = U_{th}/N_{i}k_{B}T as a function of log (T/T_{p}) = −log η at β = ħω_{ci}/k_{B}T_{p} = 0, 1, 10, 100, and 10^{3} (numbers near the lines). The analytical approximation in Eq. (76) (dotted lines) and in Eq. (78) (shortdashed lines) are compared with the numerical results of Baiko (2009) (solid lines for β = 1, 10, and 100). 

Open with DEXTER 
Fig. 3 Thermal phonon contribution to the reduced entropy s_{th} = S_{th}/N_{i}k_{B}T as a function of log (T/T_{p}) at β = ħω_{ci}/k_{B}T_{p} = 0, 0.1, 1, 10, 100, and 10^{4} (numbers near the lines). The analytical approximations in Eq. (77) (dotted lines) and in Eq. (79) (dashed lines) are compared with the numerical results of Baiko (2009) (solid lines for ). The inset shows the same approximations in the linear scale for s_{th}; here, the solid line corresponds to the numerical data at β = 100. 

Open with DEXTER 
Fig. 4 Thermal phonon contribution to the reduced heat capacity C_{V,lat}/N_{i}k_{B}T as a function of log (T/T_{p}) at β = ħω_{ci}/k_{B}T_{p} = 0, 1, 10, 100, and 10^{3} (numbers near the lines). The analytical approximation in Eq. (80) (shortdashed lines) is compared with the numerical results of Baiko (2009) (solid lines for β = 0, 1, 10, and 100). The dotted lines correspond to the first term on the r.h.s. of Eq. (80). 

Open with DEXTER 
6.2. Zeropoint 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 zeropoint 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 (82)and in Eq. (32) we have , where we have defined (83)The reduced frequency moment still depends on β, because the character of ion vibrations is affected by the magnetic field (they become essentially onedimensional if B is extremely large), but the latter dependence is relatively weak. Having extracted from the available numerical results for u_{1} (Baiko 2009; Baiko & Yakovlev 2012), we can represent it by the simple interpolation (84)where is the zerofield value and is the limit of at β → ∞. Only one of the three phonon branches contributes to in the latter limit, therefore . For the bcc crystal, , whereas varies between 0.18 and 0.19 depending on the orientation of the lattice in the magnetic field. In Fig. 5 we show and the logarithm (base 10) of u_{1} versus β.
Fig. 5 Reduced first moment of phonon frequencies (solid line) and log u_{1} (dashed line) as functions of β = ħω_{ci}/T_{p}. The dotted line shows the asymptote log (β/3). 

Open with DEXTER 
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_{1} between the values of u_{1} for two orientations, where the field lines connect an ion with its nearest neighbor in the first case and with a nextorder nearest neighbor in the second case, is approximately (85)with saturation level (Δu_{1})_{max} = 0.0064 for the bcc lattice.
6.3. Comparison with the BaikoYakovlev 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 BaikoYakovlev 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 fullyionized iron (Z = 26, A = 56) at T = 10^{7} K and B = 10^{12} G (for illustration, the density range is extended to ρ ≲ 10^{5} neglecting the bound states that can be important in this ρ – T domain). We plot the normalized pressure p = P/n_{i}k_{B}T, entropy S/N_{i}k_{B}, heat capacity c_{V} = C_{V}/N_{i}k_{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 () and at B = 10^{12} 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 idealgas 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^{12} G. The latter difference is because at ρ < ρ_{B} the electron gas is effectively onedimensional due to the strong magnetic quantization.
With increasing density, the reduced pressure p first decreases below its idealgas 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^{4} g cm^{3}, the thermodynamic functions approach their zerofield values. The gradually decreasing oscillations correspond to consecutive filling of the electron Landau levels. The magnetic field B = 10^{12} G does not affect the ion contributions in Fig. 6, because it is nonquantizing for the iron nuclei at T = 10^{7} K (ζ_{i} = 0.00342).
The liquidsolid phase transition occurs in Fig. 6 at ρ ≈ 8.25 × 10^{4} g cm^{3}, where we adopt the classical OCP melting condition Γ = 175.2 (Paper I). With further increase in density (ρ ≳ 10^{6}) 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.
Fig. 6 Reduced thermodynamic functions P/n_{i}k_{B}T, S/N_{i}k_{B}, C_{V}/N_{i}k_{B}, χ_{ρ}, and χ_{T} for a fullyionized nonmagnetic (dashed lines) and magnetized (B = 10^{12} K, solid lines) iron plasma at T = 10^{7} K. The vertical dotted lines mark the densities at which (1) ; (2) T_{F} = T; (3) ρ = ρ_{B}; (4) Γ = Γ_{m}; and (5) T_{p} = T. 

Open with DEXTER 
Fig. 7 Characteristics of the melting transition of nonideal carbon and iron plasmas at different field strengths B (marked near the curves). Lower panel: the value Γ_{m} of the Coulomb coupling parameter Γ at the melting point as function of mass density ρ. Upper panel: normalized latent heat per ion at the melting transition. The dotdashed and dashed segments of the curves correspond to the domains of nonperturbative quantum effects (T < 0.5 T_{p}) and electron response (Z^{2} Ry > 0.1 ϵ_{F}), respectively. The dotted horizontal lines mark the OCP values. The filled and open circles mark the positions of the real and virtual condensed surfaces (see text in Sect. 7.3). 

Open with DEXTER 
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_{i}k_{B}T. We plot the data for fully ionized ^{12}C and ^{56}Fe at B = 0 and 10^{13} G for carbon, B = 0, 10^{14} G, and 10^{15} 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 dotdashed 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 hydrogenlike ion exceeds 10% of the Fermi energy, Z^{2} Ry > 0.1 ϵ_{F}, which corresponds to a_{e} ≳ 0.6a_{0}/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 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^{7}–10^{8} g cm^{3} corresponding to the “sensitivity strip” in the neutronstar cooling theory (Yakovlev & Pethick 2004), the stabilization proves to be significant at the magnetar field strengths B = 10^{14}–10^{15} G. The results for ^{56}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^{12}–10^{13} G, the effect is noticeable at lower densities. However, these conclusions remain preliminary in the absence of an evaluation of the magneticfield 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 (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 ionsphere model (Salpeter 1961), which is close to the uniform model of Fushiki et al. (1989). For comparison, the results of the zerotemperature ThomasFermi model for ^{56}Fe at G (Rögnvaldsson et al. 1993) can be approximated (within 4%) by ρ_{s,ξ} with , whereas the finitetemperature ThomasFermi 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 K and critical density at 1 ≲ B_{12} ≲ 10^{3}. According to another study (Lai & Salpeter 1997; Lai 2001), T_{crit} for hydrogen is several times smaller.
Medin & Lai (2006) performed densityfunctional 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 densityfunctional calculations of other authors prompts that Q_{s} may vary within a factor of two at B_{12} ≳ 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 selfconsistently, 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^{2} (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 ^{12}C at and ^{56}Fe at by Eq. (86) with 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 K at 1 ≲ B_{12} ≲ 1000. For comparison, the results of Lai & Salpeter (1997) for H at 10 ≲ B_{12} ≲ 500 suggest 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 fullyionized plasma model. In this model 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})^{5}] (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^{6} 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 neutronstar envelopes with strong magnetic fields. Figure 8 illustrates the structure of a typical magnetar envelope with the groundstate 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^{15} 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.
Fig. 8 Structure of a magnetar envelope having the groundstate nuclear composition, the effective temperature 10^{6.5} K, and magnetic field B = 10^{15} G, inclined at 45° to the surface. Top panel: thermal profiles calculated using the present EOS (solid red line) and the EOSs where the Coulomb nonideality is either neglected (dotted green line) or treated without account of magnetic quantization (dashed blue line, which is superposed on the solid red line). The melting temperature is drawn by the oblique dotdashed line. Thin vertical dotdashed lines mark the points of phase transitions from solid to liquid and back to solid state. Asterisks mark the ends of the convective segment, which is indicated by the arrow. Middle panel: ion charge (Rüster et al. 2006). Bottom panel: reduced pressure, heat capacity, and logarithmic derivatives of pressure. 

Open with DEXTER 
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 heattransport 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 turningup “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 ^{66}Ni and ^{86}Kr at ρ = 1.5 × 10^{9} g cm^{3} (z = 73.8 m). These phase transitions do not cause any substantial breaks in χ_{ρ}, χ_{T}, or C_{V}/N_{i}k_{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_{i}k_{B}T at the interfaces ^{66}Ni/^{86}Kr (ρ = 1.5 × 10^{9} g cm^{3}) and ^{78}Ni/^{124}Mo (ρ = 1.32 × 10^{11} g cm^{3}), which are caused by the decreases in n_{i} with the large jumps in A (of course, the nonnormalized pressure P is continuous). The specific heat per ion C_{V}/N_{i}k_{B} is almost continuous at these interfaces, which means that heat capacity of unit volume abruptly decreases. The drop in χ_{T} at the ^{78}Ni/^{128}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 fullyionized electronion plasmas in magnetic fields and described several improvements to the previously published approximations, taking nonideality attributable to ionion, electronelectron, and electronion 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 110200253a), and the Russian Leading Scientific Schools program (grant NSh3769.2010.2).
References
 Alastuey, A., & Jancovici, B. 1980, Phys. A, 102, 327 [CrossRef] [Google Scholar]
 Baiko, D. A. 2000, Ph.D. Thesis (St. Petersburg: Ioffe Phys.Tech. Inst.) [in Russian] [Google Scholar]
 Baiko, D. A. 2002, Phys. Rev. E, 66, 056405 [NASA ADS] [CrossRef] [Google Scholar]
 Baiko, D. A. 2009, Phys. Rev. E, 80, 046405 [NASA ADS] [CrossRef] [Google Scholar]
 Baiko, D. A., Kaminker, A. D., Potekhin, A. Y., & Yakovlev, D. G. 1998, Phys. Rev. Lett., 81, 5556 [NASA ADS] [CrossRef] [Google Scholar]
 Baiko, D. A., Potekhin, A. Y., & Yakovlev, D. G. 2001, Phys. Rev. E, 64, 057402 [NASA ADS] [CrossRef] [Google Scholar]
 Banerjee, B., Constantinescu, D. H., & Rehak, P. 1974, Phys. Rev. D, 10, 2384 [NASA ADS] [CrossRef] [Google Scholar]
 Berestetskiĭ, V. B., Lifshitz, E. M., & Pitaevskiĭ, L. P., 1982, Quantum Electrodynamics (Oxford: ButterworthHeinemann) [Google Scholar]
 Blandford, R. D., & Hernquist, L. 1982, J. Phys. C: Solid State Phys., 15, 6233 [NASA ADS] [CrossRef] [Google Scholar]
 Blinnikov, S. I., DuninaBarkovskaya, N. V., & Nadyozhin, D. K. 1996, ApJS, 106 171; erratum: 1998, ApJS, 118, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Brami, B., Hansen, J. P., & Joly, F. 1979, Physica A, 95, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Broderick, A., Prakash, M., & Lattimer, J. M. 2000, ApJ, 537, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Caillol, J. M. 1999, J. Chem. Phys., 111, 6538 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V., & Chiu, H.Y. 1968, Phys. Rev., 173, 1220 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G., & Ashcroft, N. W. 1990, Phys. Rev. A, 42, 2284 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G., & Baraffe, I. 2000, ARA&A, 38, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G., & Potekhin, A. Y. 1998, Phys. Rev. E, 58, 4941 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G., Douchin, F., & Potekhin, A. Y. 2002, J. Phys., Condensed Matter, 14, 9133 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1957, An Introduction to the Study of Stellar Structure (New York: Dover, 1957), Chap. X [Google Scholar]
 Cornu, F. 1998, Phys. Rev. E, 58, 5293 [NASA ADS] [CrossRef] [Google Scholar]
 Danz, R. W., & Glasser, M. L. 1971, Phys. Rev. B, 4, 94 [NASA ADS] [CrossRef] [Google Scholar]
 DeWitt, H. E., & Slattery, W. 2003, Contrib. Plasma Phys., 43, 279 [NASA ADS] [CrossRef] [Google Scholar]
 DeWitt, H. E., Slattery, W., & Chabrier, G. 1996, Phys. A, 228, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Farouki, R. T., & Hamaguchi, S. 1993, Phys. Rev. E, 47, 4330 [NASA ADS] [CrossRef] [Google Scholar]
 Fortov, V. E. 2009, Phys. Usp., 52, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Fushiki, I., Gudmundsson, E. H., & Pethick, C. J. 1989, ApJ, 342, 958 [NASA ADS] [CrossRef] [Google Scholar]
 Galam, S., & Hansen, J. P. 1976, Phys. Rev. A, 14, 816 [NASA ADS] [CrossRef] [Google Scholar]
 Girifalco, L. A. 1973, Statistical Physics of Materials (New York: WileyInterscience, 1973), App. V [Google Scholar]
 Griffith, D. J. 1999, Introduction to Electrodynamics, 3rd ed. (London: PrenticeHall), Chap. 6. [Google Scholar]
 Haensel, P., Potekhin, A. Y., & Yakovlev, D. G. 2007, Neutron Stars 1: Equation of State and Structure (New York: Springer) [Google Scholar]
 Hamaguchi, S., Farouki, R. T., & Dubin, D. H. E. 1997, Phys. Rev. E, 56, 4671 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, B. M. S. 2004, Phys. Rep., 399, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, J. P., & Vieillefosse, P. 1975, Phys. Lett. A, 53, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, J. P., Torrie, G. M., & Vieillefosse, P. 1977, Phys. Rev. A, 16, 2153 [NASA ADS] [CrossRef] [Google Scholar]
 Hughto, J., Horowitz, C. J., Schneider, A. S., et al. 2012, Phys. Rev. E, accepted [arXiv:1211.0891] [Google Scholar]
 Ichimaru, S., Iyetomi, H., & Tanaka, S. 1987, Phys. Rep., 149, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Jancovici, B. 1962, Nuovo Cimento, 25, 428 [Google Scholar]
 Johnson, M. H., & Lippmann, B. A. 1949, Phys. Rev., 76, 828 [NASA ADS] [CrossRef] [Google Scholar]
 Kaplan, J. L., & Glasser, M. L. 1972, Phys. Rev. Lett., 28, 1077 [NASA ADS] [CrossRef] [Google Scholar]
 Kittel, C. 1963, Quantum Theory of Solids (New York: Wiley) [Google Scholar]
 Lai, D. 2001, Rev. Mod. Phys., 73, 629 [NASA ADS] [CrossRef] [Google Scholar]
 Lai, D., & Salpeter, E. E. 1997, ApJ, 491, 270 [NASA ADS] [CrossRef] [Google Scholar]
 Landau, L. D. 1930, Z. f. Physik, 64, 629 [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1977, Quantum Mechanics. NonRelativistic Theory, 3rd edn. (Oxford: Pergamon) [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1980, Statistical Physics, 3rd edn. (Oxford: ButterworthHeinemann) [Google Scholar]
 Medin, Z., & Cumming, A. 2010, Phys. Rev. E, 81, 036107 [NASA ADS] [CrossRef] [Google Scholar]
 Medin, Z., & Lai, D. 2006, Phys. Rev. A, 74, 062508 [NASA ADS] [CrossRef] [Google Scholar]
 Medin, Z., & Lai, D. 2007, MNRAS, 382, 1833 [NASA ADS] [Google Scholar]
 Militzer, B., & Graham, R. L. 2006, J. Phys. Chem. Sol., 67, 2136 [NASA ADS] [CrossRef] [Google Scholar]
 Morbec, J. M., & Capelle, K. 2008, Phys. Rev. B, 78, 085107 [NASA ADS] [CrossRef] [Google Scholar]
 Nagai, T., & Fukuyama, H. 1982, J. Phys. Soc. Japan, 51, 3431 [NASA ADS] [CrossRef] [Google Scholar]
 Nagai, T., & Fukuyama, H. 1983, J. Phys. Soc. Japan, 52, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Ogata, S., Iyetomi, H., Ichimaru, S., & Van Horn, H. M. 1993, Phys. Rev. E, 48, 1344 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pollock, E. L., & Hansen, J. P. 1973, Phys. Rev. A, 8, 3110 [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., & Chabrier, G. 2000, Phys. Rev. E, 62, 8554 (Paper I) [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., & Chabrier, G. 2004, ApJ, 600, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., & Chabrier, G. 2010, Contrib. Plasma Phys., 50, 82 (Paper II) [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., Chabrier, G., & Shibanov, Yu. A. 1999, Phys. Rev. E, 60, 2193 [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., Chabrier, G., & Yakovlev, D. G. 2007, Ap&SS, 308, 353 (electronic version with corrected misprints [arXiv:astroph/0611014v3] [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., Chabrier, G., & Rogers, F. J. 2009a, Phys. Rev. E, 79, 016411 [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., Chabrier, G., Chugunov, A. I., DeWitt, H. E., & Rogers, F. J. 2009b, Phys. Rev. E, 80, 047401 [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., Suleimanov, V. F., van Adelsberg, M., & Werner, K. 2012, A&A, 546, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rögnvaldsson, Ö. E., Fushiki, I., Gudmundsson, E. H., Pethick, C. J., & Yngvason, J. 1993, ApJ, 416, 276 [NASA ADS] [CrossRef] [Google Scholar]
 Rosenfeld, Y. 1996, Phys. Rev. E, 54, 2827 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. A. 1971, Phys. Rev. Lett., 27, 1306 [NASA ADS] [CrossRef] [Google Scholar]
 Rüster, S. B., Hempel, M., & SchaffnerBielich, J. 2006, Phys. Rev. C, 73, 035804 [NASA ADS] [CrossRef] [Google Scholar]
 Salpeter, E. E. 1961 ApJ, 134, 669 [Google Scholar]
 Schwarzschild, M. 1958, Structure and Evolution of the Stars (Princeton: Princeton Univ. Press) [Google Scholar]
 Schwinger, J. 1988, Particles, Sources, and Fields (Redwood: AddisonWesley) [Google Scholar]
 Sharma, R., & Reddy, S. 2011, Phys. Rev. C, 83, 025803 [NASA ADS] [CrossRef] [Google Scholar]
 Steinberg, M., Ortner, J., & Ebeling, W. 1998, Phys. Rev. E, 58, 3806 [NASA ADS] [CrossRef] [Google Scholar]
 Steinberg, M., Ebeling, W., & Ortner, J. 2000, Phys. Rev. E, 61, 2290 [NASA ADS] [CrossRef] [Google Scholar]
 Suh, I.S., & Mathews, G. J. 2001, ApJ, 546, 1126 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, S., Mitake, S., & Ichimaru, S. 1896, Phys. Rev. A, 32 [Google Scholar]
 Thorolfsson, A., Rögnvaldsson, Ö. E., Yngvason, J., & Gudmundsson, E. H. 1998, ApJ, 502, 847 [NASA ADS] [CrossRef] [Google Scholar]
 Timmes, F. X., & Arnett, D. 1999, ApJS, 125, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501 [NASA ADS] [CrossRef] [Google Scholar]
 Usov, N. A., Grebenshchikov, Yu. B., & Ulinich, F. R. 1980, Sov. Phys.–JETP, 51, 148 [NASA ADS] [Google Scholar]
 van Leeuwen, J. H. 1921, J. Phys. Radium, Ser. VI, 2, 361 [Google Scholar]
 Wickramasinghe, D. T., & Ferrario, L. 2000, PASP, 112, 873 [NASA ADS] [CrossRef] [Google Scholar]
 Wigner, E. P. 1932, Phys. Rev. 40, 749 [Google Scholar]
 Yakovlev, D. G., & Kaminker, A. D. 1994, in The Equation of State in Astrophysics, eds. G. Chabrier, & E. Schatzman (Cambridge: Cambridge University Press), Proc. IAU Colloq., 147, 214 [Google Scholar]
 Yakovlev, D. G., & Pethick, C. J. 2004, ARA&A, 42, 169 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Yakovlev, D. G., & Shalybkov, D. A. 1989, Astrophys. Space Phys. Rev. 7, 311 [Google Scholar]
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 (A.1)where A = b − gd, B = a − g, C = 2 − d^{2}/f, , and a, b, d, f, and g are the following functions of θ = T/T_{F} (at B = 0): 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_{i}k_{B}T of the classic OCP, we have the following analytical formula (Paper I): (B.1)where A_{1} = −0.907347, A_{2} = 0.62849, B_{1} = 0.0045, B_{2} = 170, B_{3} = −8.4 × 10^{5}, B_{4} = 0.0037, and . The derivative (B.2)reproduces Monte Carlo calculations of the reduced internal energy U_{ii}/N_{i}k_{B}T 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 WignerKirkwood quantum corrections (Wigner 1932; Landau & Lifshitz 1980). The lowestorder correction to the reduced free energy is (B.3)The nextorder correction ∝ ħ^{4} was obtained by Hansen & Vieillefosse (1975). It can be written as (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 lowestorder correction (B.3), i.e., . 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} = −0.895 929 255 68 and u_{1} = 0.511 3875, and for f_{th} the following fitting formula can be used: (B.5)where α_{1} = 0.932446, α_{2} = 0.334547, α_{3} = 0.265764, The Taylor expansion of Eq. (B.5) at small η is consistent with the Wigner correction (B.3). However, the next Taylor term ~η^{3} is absent in the Wigner expansion, and therefore Eq. (B.5) does not reproduce higherorder 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^{5}. 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): (B.6)where a_{1} = 10.9, a_{2} = 247, and a_{3} = 1.765 × 10^{5}. A continuation to arbitrary η, which is consistent with available analytical and numerical results for quantum crystals, reads (Paper II) (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) (C.1)where ensures exact transition to the DebyeHückel limit at Γ → 0, c_{TF} = (18/175) (12/π)^{2/3}Z^{7/3}(1 − Z^{ − 1/3} + 0.2 Z^{−1/2}( fits the numerical data at large Γ and reproduces the ThomasFermi limit (Salpeter 1961) at Z → ∞, the parameters a = 1.11 Z^{0.475}, b = 0.2 + 0.078 (lnZ)^{2}, and ν = 1.16 + 0.08lnZ provide a loworder approximation to F_{ie} for intermediate r_{s} and Γ. The functions improve the fit at relatively large r_{s}. Finally, the function is the relativistic correction, as is 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) (C.2)where the parameter a_{TF} = 0.00352 is related to c_{TF} in Eq. (C.1), and parameters s and b_{1} – b_{4} depend on Z: Here, the numerical parameters are given for the bcc crystal; their values for the fcc lattice are slightly different (Paper I).
All Figures
Fig. 1 Characteristic densitytemperature domains at B = 10^{12} G (blue online) and 10^{15} G (red online) for fullyionized iron. Solid lines indicate the Fermi temperature as function of density, the dotted line shows the plasma temperature, the dotdashed line shows the melting temperature as function of density, short and long dashes delimit the domains of strongly quantizing magnetic field and of magnetic condensation, respectively, and the heavy dots mark the critical point for the condensation (Sect. 7.3). 

Open with DEXTER  
In the text 
Fig. 2 Thermal phonon contribution to the reduced internal energy u_{th} = U_{th}/N_{i}k_{B}T as a function of log (T/T_{p}) = −log η at β = ħω_{ci}/k_{B}T_{p} = 0, 1, 10, 100, and 10^{3} (numbers near the lines). The analytical approximation in Eq. (76) (dotted lines) and in Eq. (78) (shortdashed lines) are compared with the numerical results of Baiko (2009) (solid lines for β = 1, 10, and 100). 

Open with DEXTER  
In the text 
Fig. 3 Thermal phonon contribution to the reduced entropy s_{th} = S_{th}/N_{i}k_{B}T as a function of log (T/T_{p}) at β = ħω_{ci}/k_{B}T_{p} = 0, 0.1, 1, 10, 100, and 10^{4} (numbers near the lines). The analytical approximations in Eq. (77) (dotted lines) and in Eq. (79) (dashed lines) are compared with the numerical results of Baiko (2009) (solid lines for ). The inset shows the same approximations in the linear scale for s_{th}; here, the solid line corresponds to the numerical data at β = 100. 

Open with DEXTER  
In the text 
Fig. 4 Thermal phonon contribution to the reduced heat capacity C_{V,lat}/N_{i}k_{B}T as a function of log (T/T_{p}) at β = ħω_{ci}/k_{B}T_{p} = 0, 1, 10, 100, and 10^{3} (numbers near the lines). The analytical approximation in Eq. (80) (shortdashed lines) is compared with the numerical results of Baiko (2009) (solid lines for β = 0, 1, 10, and 100). The dotted lines correspond to the first term on the r.h.s. of Eq. (80). 

Open with DEXTER  
In the text 
Fig. 5 Reduced first moment of phonon frequencies (solid line) and log u_{1} (dashed line) as functions of β = ħω_{ci}/T_{p}. The dotted line shows the asymptote log (β/3). 

Open with DEXTER  
In the text 
Fig. 6 Reduced thermodynamic functions P/n_{i}k_{B}T, S/N_{i}k_{B}, C_{V}/N_{i}k_{B}, χ_{ρ}, and χ_{T} for a fullyionized nonmagnetic (dashed lines) and magnetized (B = 10^{12} K, solid lines) iron plasma at T = 10^{7} K. The vertical dotted lines mark the densities at which (1) ; (2) T_{F} = T; (3) ρ = ρ_{B}; (4) Γ = Γ_{m}; and (5) T_{p} = T. 

Open with DEXTER  
In the text 
Fig. 7 Characteristics of the melting transition of nonideal carbon and iron plasmas at different field strengths B (marked near the curves). Lower panel: the value Γ_{m} of the Coulomb coupling parameter Γ at the melting point as function of mass density ρ. Upper panel: normalized latent heat per ion at the melting transition. The dotdashed and dashed segments of the curves correspond to the domains of nonperturbative quantum effects (T < 0.5 T_{p}) and electron response (Z^{2} Ry > 0.1 ϵ_{F}), respectively. The dotted horizontal lines mark the OCP values. The filled and open circles mark the positions of the real and virtual condensed surfaces (see text in Sect. 7.3). 

Open with DEXTER  
In the text 
Fig. 8 Structure of a magnetar envelope having the groundstate nuclear composition, the effective temperature 10^{6.5} K, and magnetic field B = 10^{15} G, inclined at 45° to the surface. Top panel: thermal profiles calculated using the present EOS (solid red line) and the EOSs where the Coulomb nonideality is either neglected (dotted green line) or treated without account of magnetic quantization (dashed blue line, which is superposed on the solid red line). The melting temperature is drawn by the oblique dotdashed line. Thin vertical dotdashed lines mark the points of phase transitions from solid to liquid and back to solid state. Asterisks mark the ends of the convective segment, which is indicated by the arrow. Middle panel: ion charge (Rüster et al. 2006). Bottom panel: reduced pressure, heat capacity, and logarithmic derivatives of pressure. 

Open with DEXTER  
In the text 