EDP Sciences
Free Access
Issue
A&A
Volume 550, February 2013
Article Number A43
Number of page(s) 16
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201220082
Published online 23 January 2013

© 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 107−109 G (e.g., Wickramasinghe & Ferrario 2000, and references therein) and neutron stars with typical B ~ 108−1014 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 ne and ni 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 ne = Zni) 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 ne and temperature T. Instead of ne it is convenient to introduce the dimensionless density parameter rs = ae/a0, where a0 is the Bohr radius and . The parameter rs can be quickly evaluated from the relations where n24 ≡ ne/1024   cm-3 and ρ0 = 2.6752   (A/Z) g cm-3. The analogous density parameter for the ions is RS = aimi(Ze)22 = 1822.89AZ7/3   rs, where mi is the ion mass and is the ion sphere radius.

At stellar densities it is convenient to use, instead of rs, the nonmagnetic relativity parameter (1)where pF = ħ   (3π2ne)1/3 is the electron Fermi momentum in the absence of a magnetic field, and ρ6 ≡ ρ/106 g cm-3. The Fermi energy (without the rest energy) for the electron gas is and the Fermi temperature TF ≡ ϵF/kB = Tr   (γr − 1), where Tr ≡ mec2/kB = 5.93 × 109   K, , and kB is the Boltzmann constant. A useful measure of electron degeneracy is θ = T/TF. In the nonrelativistic limit (xr ≪ 1),  K, and (2)where (3)In the opposite ultrarelativistic limit (xr ≫ 1), θ ≈ (263   Γe)-1. The strength of the Coulomb interaction of nonrelativistic ions is characterized by the Coulomb coupling parameter (4)where T6 ≡ T/106 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 ≳ ai or at T ≪ Tp, where Tp ≡ ħωp/kB is the ion plasma temperature and ωp = (4πe2   niZ2/mi(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. Magnetic-field parameters

In the nonrelativistic theory (Landau & Lifshitz 1977), the energy of an electron in magnetic field B equals , where pz is the momentum component along B, ωc = eB/mec is the electron cyclotron frequency, characterizes a Landau level, σ =  ±1 determines the spin projection on the field, and nL 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 pz are inter-related as The levels are double-degenerate with respect to σ. Their splitting due to the anomalous magnetic moment of the electron is  ≈ (mec2αf/2π)   b at b ≪ 1 and  ~(mec2α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 kBT: (9)where B0 = 2.3505 × 109   G, (10)where αf = e2c is the fine-structure constant, and (11)where B12 ≡ B/1012 G. The magnetic length gives a characteristic transverse scale of the electron wave function.

For the ions, the cyclotron energy is ħωci = Z   (me/mi)   ħω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 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 Fii as the difference between F and  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 ai (Salpeter 1961). The electron exchange-correlation 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 ion-electron (electron polarization) contribution Fie, 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 CV = (∂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 Ni = niV nonrelativistic classical ions is (15)where Mspin 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 Ne = neV is the number of electrons and μe is the electron chemical potential without the rest energy mec2. The pressure and the number density are functions of μe and T: where χe ≡ μe/kBT, τ ≡ T/Tr, and (20)is the Fermi-Dirac 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 second-order 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 104, 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) yields1(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 104 at τ ≲ 102, 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 Fee 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 ~ Tp. 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 ≪ Tp (see Chabrier et al. 2002, for references and discussion).

3.2.2. Coulomb crystal

At T < Tm, where Tm 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 (29)where U0 = NiC0(Ze)2/ai is the classical static-lattice energy, C0 ≈ −0.9 is the Madelung constant, (30)accounts for zero-point quantum vibrations, u1 =  ⟨ ω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 classical-gas free energy, therefore Flat replaces in Eq. (14).

Beyond the harmonic-lattice approximation, the total reduced free energy flat ≡ Flat/NikBT can be written as (32)Here, the first three terms correspond to the three terms in Eq. (29), and fah is the anharmonic correction. The most accurate values of the constants C0 and u1 were calculated by Baiko (2000) (see Appendix B.2). For fth = Fth/NikBT, we use the highly precise fit of Baiko et al. (2001) (Appendix B.2). In the classical limit η ≪ 1, it reduces to fth ≃ 3lnη − 2.49389−1.5u1η + η2/24, where the term with u1 cancels that in Eq. (32) and the last term represents the Wigner-Kirkwood quantum correction (Eq. (B.3)), which is the same in the liquid and solid phases (Pollock & Hansen 1973). In the opposite limit T ≪ Tp, we have fth ≃ −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 CV = 3NikB independent of the polarization).

In Paper I we calculated Fie using the semiclassical perturbation theory of Galam & Hansen (1976) with a model structure factor, and fit the results by an analytical function of xr 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.13.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 xj are the number fractions of ions with charge numbers Zj and . In Eq. (33), fex ≡ Fex/NikBT is the reduced nonideal part of the free energy, Fex is the excess free energy, which is equal to Fii in the case of the rigid charge-neutralizing electron background and to Fii + Fie + Fee 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 − fLM, in Coulomb liquids for arbitrary coupling parameters Γj reads (Potekhin et al. 2009b) (34)where ,  ⟨ Γ ⟩  = Γe ⟨ Z5/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, Δfsol, 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 RZ = Z2/Z1 at a fixed number fraction x2 = 1 − x1 for a binary ion mixture with Z2/Z1 > 2 and x2 < 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 RZ ≲ 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 two-component ionic bcc lattices. For a multicomponent ion crystal, Medin & Cumming (2010) proposed the extrapolation from the two-component plasma case (41)where the indices are arranged so that Zj < Zj + 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 Mspin is the ion spin multiplicity, and gi is the g-factor (gi = 5.5857 for protons). For ions with spin one-half (Mspin = 2), the expression in the square brackets in Eq. (43) simplifies to  [2cosh(gi   ζi/4)] . For zero-spin ions, such as 4He, 56Fe, and other even-even nuclei in the ground state, Fspin = 0.

The ion pressure obeys the nonmagnetic ideal-gas relation , but expressions for the internal energy and heat capacity are different: Here, the terms uspin and cspin arise from Fspin, They simplify at Mspin = 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 Δpz 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 (Landau & Lifshitz 1977). Thus one can express the electron number density ne and the grand thermodynamic potential as where ϵn(pz) 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 ne, , 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 ne 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 CV 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 (53)where is the zero-field Fermi momentum at the given density. Equation (53) can be written as pF = mecxB, where (54)is the relativity parameter modified by the field and (Sect. 2.1). With increasing ne at constant B and zero temperature, the electrons start to populate the first excited Landau level when ne reaches . Therefore, the field is strongly quantizing at T ≪ Tcycl and ρ < ρB, where Tcycl = ħωc/kB ≈ 1.3434 × 108B12 K and (55)The condition ne < nB can be written as . Then Eq. (53) shows that in a strongly quantizing field , except for densities ne close to the threshold nB. Thus TF is reduced, compared to its nonmagnetic value , by factor . In the nonrelativistic limit, , and the parameter θ = T/TF becomes (56)where θ0 is the nonmagnetic value given by Eq. (2).

The opposite case of a nonquantizing magnetic field occurs at T ≫ TB, where TB 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 TB = Tcycl, if ρ ≤ ρB and TB = Tcycl/γr, if ρ > ρB (a more sophisticated but qualitatively similar definition of TB 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 = 1012 G and 1015 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 Tcycl, the dependence TF(ρ) is steeper than at B = 0, in agreement with Eq. (56). The line Tm(ρ) separates Coulomb liquid from Coulomb crystal. Near the lower right corner of the figure, where T ≪ Tp, 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).

thumbnail Fig. 1

Characteristic density-temperature domains at B = 1012 G (blue online) and 1015 G (red online) for fully-ionized iron. Solid lines indicate the Fermi temperature as function of density, the dotted line shows the plasma temperature, the dot-dashed 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 zero-temperature value and ΔF is a thermal correction. According to Eqs. (24) and (52), the zero-temperature pressure P0 is (57)where is the relativistic unit of pressure, nmax 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 higher-order 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 < TF), then In the nonrelativistic (xB ≪ 1) and ultrarelativistic (xB ≫ 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 ne ≈ nB. 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 xB replaces xr and π2/3 replaces π2.

4.2.5. Nonrelativistic limit

In the nonrelativistic limit (pF ≪ mec and T ≪ Tr), Eqs. (51) and (52) simplify to (69)In the nondegenerate regime (T ≫ TF), one has Iν(χ) ≈ eχ   Γ(ν + 1), where Γ(ν + 1) is the gamma-function. 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 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  is recovered. In the strongly quantizing, nondegenerate regime (ρ < ρB and TF ≪ T ≪ Tcycl), 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 Pkin, 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 jm = 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 jm × 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: . 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 field-free expression for Fee, which matches available exact limiting expressions, including the cases of nonquantizing, strongly quantizing degenerate, and strongly quantizing nondegenerate plasmas. The modification consists in replacing Fee(θ,Γe) by Fee(θ ∗ e), where (74)ξ =  [1−(4/ζ)tanh(ζ/4)] 1/2, and θ0 and θB are given by Eqs. (2) and (56), respectively, at fixed ne 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): (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. 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 () increases the polarization energy at high densities (rs ≪ 1) by a factor of (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 πħ/pF in a cylinder of a radius  ~πħ/pF 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 ≪ TF, i.e., deep within the triangular domains formed by the lines TF and ρB in Fig. 1. At the typical pulsar magnetic fields B ~ 1012 G, this requires an unusually low temperature of the neutron-star crust. On the other hand, the conditions T ≪ TF and ρ < ρB can be easily fulfilled in the outer crust of magnetars at B ~ 1015 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 Fie 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 RS 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α ∝ k2 near the center of the Brillouin zone, which leads to the unusual dependence of the heat capacity of the lattice CV,lat ∝ T3/2 at T → 0 instead of the Debye law CV,lat ∝ T3. 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 fth to the reduced free energy of a Coulomb crystal F/NikBT 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 > kBT, 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 fth in Eq. (32) can be rewritten as fth = uth − sth, where uth = Uth/NikBT and sth = Sth/NikB are the thermal contributions to the reduced internal energy and entropy. We approximately represent uth by the function (76)where ψ = 12.5   (β/η)3/2 + 119   (β/η)2 and ζi = βη, and we represent sth by the function (77)In these equations, and are the values of uth and sth at β = 0. Equations (76) and (77) exactly reproduce the known asymptotic limits: uth = 3 in the classical nonmagnetic limit (η ≪ β ≪ 1), uth = 2 in the classical magnetic limit (η ≪ β-1 ≪ 1), uth = 1 in the case where β ≫ η ≫ 1, and uth = 0.6sth ∝ η−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 fth(η,β). In this way we obtain, for example, Note that the relation between the phonon contributions to pressure and internal energy, 2PthV = Uth, 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. 24. Their reasonable behavior beyond the range of available numerical data is demonstrated by plotting them also at larger β = 103 and 104.

thumbnail Fig. 2

Thermal phonon contribution to the reduced internal energy uth = Uth/NikBT as a function of log (T/Tp) = −log η at β = ħωci/kBTp = 0, 1, 10, 100, and 103 (numbers near the lines). The analytical approximation in Eq. (76) (dotted lines) and in Eq. (78) (short-dashed lines) are compared with the numerical results of Baiko (2009) (solid lines for β = 1, 10, and 100).

Open with DEXTER

thumbnail Fig. 3

Thermal phonon contribution to the reduced entropy sth = Sth/NikBT as a function of log    (T/Tp) at β = ħωci/kBTp = 0, 0.1, 1, 10, 100, and 104 (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 sth; here, the solid line corresponds to the numerical data at β = 100.

Open with DEXTER

thumbnail Fig. 4

Thermal phonon contribution to the reduced heat capacity CV,lat/NikBT as a function of log (T/Tp) at β = ħωci/kBTp = 0, 1, 10, 100, and 103 (numbers near the lines). The analytical approximation in Eq. (80) (short-dashed 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. 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 Uq is given by Eq. (30).

In the case where the crystal is placed in a magnetic field, Uq includes contributions due to both magnetic and lattice confinements of the ion motion. However, since the magnetic contribution Niħω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 one-dimensional if B is extremely large), but the latter dependence is relatively weak. Having extracted from the available numerical results for u1 (Baiko 2009; Baiko & Yakovlev 2012), we can represent it by the simple interpolation (84)where is the zero-field 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 u1 versus β.

thumbnail Fig. 5

Reduced first moment of phonon frequencies (solid line) and log u1 (dashed line) as functions of β = ħωci/Tp. 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 Δu1 between the values of u1 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 (85)with saturation level (Δu1)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 ≪ Tp and T ≫ Tp and have similar accuracies within several percent points in the intermediate range 0.1Tp/β ≲ T ≲ 10Tp. 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 = 107 K and B = 1012 G (for illustration, the density range is extended to ρ ≲ 105 neglecting the bound states that can be important in this ρ – T domain). We plot the normalized pressure p = P/nikBT, entropy S/NikB, heat capacity cV = CV/NikB, 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 = 1012 G (TF = T), population of excited Landau levels (ρ = ρB), melting point with formation of a classical Coulomb crystal (Tm = T), and quantum effects in the crystal (Tp = T).

At low densities, the ideal-gas values are approached: p = 1 + Z, χρ = χT = 1, cV = (3 + 3Z)/2 at B = 0, and cV = (3 + Z)/2 at B = 1012 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 × 104 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 = 1012 G does not affect the ion contributions in Fig. 6, because it is nonquantizing for the iron nuclei at T = 107 K (ζi = 0.00342).

The liquid-solid phase transition occurs in Fig. 6 at ρ ≈ 8.25 × 104 g cm-3, where we adopt the classical OCP melting condition Γ = 175.2 (Paper I). With further increase in density (ρ ≳ 106) 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 cV = 3 characteristic of the classical simple crystal. At still higher density the ion motions become quantized (Tp ≫ T) which leads to the further decrease in the heat capacity and the entropy.

thumbnail Fig. 6

Reduced thermodynamic functions P/nikBT, S/NikB, CV/NikB, χρ, and χT for a fully-ionized nonmagnetic (dashed lines) and magnetized (B = 1012 K, solid lines) iron plasma at T = 107 K. The vertical dotted lines mark the densities at which (1) ; (2) TF = T; (3) ρ = ρB; (4) Γ = Γm; and (5) Tp = T.

Open with DEXTER

thumbnail 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 dot-dashed and dashed segments of the curves correspond to the domains of nonperturbative quantum effects (T < 0.5   Tp) and electron response (Z2   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 Qm = Uliq − Usol), divided by NikBT. We plot the data for fully ionized 12C and 56Fe at B = 0 and 1013 G for carbon, B = 0, 1014 G, and 1015 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   Tp. 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, Z2   Ry > 0.1   ϵF, which corresponds to ae ≳ 0.6a0/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 Qm 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 ρ ~ 107–108 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 = 1014–1015 G. The results for 56Fe 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 = 1012–1013 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 Tcrit, 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 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 56Fe at  G (Rögnvaldsson et al. 1993) can be approximated (within 4%) by ρs with , 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  K and critical density at 1 ≲ B12 ≲ 103. According to another study (Lai & Salpeter 1997; Lai 2001), Tcrit for hydrogen is several times smaller.

Medin & Lai (2006) performed density-functional calculations of the cohesive energy Qs 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 Qs may vary within a factor of two at B12 ≳ 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 Tcrit 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 = mi/4aR2 (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 12C at  and 56Fe at by Eq. (86) with and ξ = 0.55 ± 0.11, respectively.

Medin & Lai (2007) found that the critical temperature is Tcrit ≈ 0.08Qs/kB. Their numerical results for He, C, and Fe can be roughly (within a factor of 1.5) described as  K at 1 ≲ B12 ≲ 1000. For comparison, the results of Lai & Salpeter (1997) for H at 10 ≲ B12 ≲ 500 suggest  K. The discrepancies between different estimates of ρs and Tcrit 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  K and ρcrit ≈ ρs,0.47. With decreasing temperature T below Tcrit, the surface density increases and tends to the limit ρs,1 given using Eq. (86) as ρs ≈ ρs,1/ [1 + 1.1   (T/Tcrit)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 Tm > Tcrit. In this case, the open circles mark the density that the condensed surface would have at much smaller temperatures T ~ 106   K ≪ Tcrit. 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 106.5 K and magnetic field B = 1015 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.

thumbnail Fig. 8

Structure of a magnetar envelope having the ground-state nuclear composition, the effective temperature 106.5 K, and magnetic field B = 1015 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 dot-dashed line. Thin vertical dot-dashed 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 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 66Ni and 86Kr at ρ = 1.5 × 109 g cm-3 (z = 73.8 m). These phase transitions do not cause any substantial breaks in χρ, χT, or CV/NikB, 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/nikBT at the interfaces 66Ni/86Kr (ρ = 1.5 × 109 g cm-3) and 78Ni/124Mo (ρ = 1.32 × 1011 g cm-3), which are caused by the decreases in ni with the large jumps in A (of course, the non-normalized pressure P is continuous). The specific heat per ion CV/NikB is almost continuous at these interfaces, which means that heat capacity of unit volume abruptly decreases. The drop in χT at the 78Ni/128Mo interface is due to the same decrease in ni, 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.


1

The multiplier was accidentally omitted in Paper II.

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 ≪ TF was improved by Ichimaru et al. (1987A.1). The result reads (A.1)where A = b − gd, B = a − g, C = 2 − d2/f, , and a, b, d, f, and g are the following functions of θ = T/TF (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 fii ≡ Fii/NikBT of the classic OCP, we have the following analytical formula (Paper I): (B.1)where A1 = −0.907347, A2 = 0.62849, B1 = 0.0045, B2 = 170, B3 = −8.4 × 10-5, B4 = 0.0037, and . The derivative (B.2)reproduces Monte Carlo calculations of the reduced internal energy Uii/NikBT 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 ≫ Tp. One can extend the applicability range of the analytical EOS to T ~ Tp by the Wigner-Kirkwood quantum corrections (Wigner 1932; Landau & Lifshitz 1980). The lowest-order correction to the reduced free energy is (B.3)The next-order 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 lowest-order 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 C0 = −0.895   929   255   68 and u1 = 0.511   3875, and for fth 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 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 105. 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 a1 = 10.9, a2 = 247, and a3 = 1.765 × 105. 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 Debye-Hückel limit at Γ → 0, cTF = (18/175)   (12/π)2/3Z7/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   Z0.475, b = 0.2 + 0.078   (lnZ)2, and ν = 1.16 + 0.08lnZ provide a low-order approximation to Fie for intermediate rs and Γ. The functions improve the fit at relatively large rs. 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 aTF = 0.00352 is related to cTF in Eq. (C.1), and parameters s and b1 – b4 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

thumbnail Fig. 1

Characteristic density-temperature domains at B = 1012 G (blue online) and 1015 G (red online) for fully-ionized iron. Solid lines indicate the Fermi temperature as function of density, the dotted line shows the plasma temperature, the dot-dashed 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
thumbnail Fig. 2

Thermal phonon contribution to the reduced internal energy uth = Uth/NikBT as a function of log (T/Tp) = −log η at β = ħωci/kBTp = 0, 1, 10, 100, and 103 (numbers near the lines). The analytical approximation in Eq. (76) (dotted lines) and in Eq. (78) (short-dashed lines) are compared with the numerical results of Baiko (2009) (solid lines for β = 1, 10, and 100).

Open with DEXTER
In the text
thumbnail Fig. 3

Thermal phonon contribution to the reduced entropy sth = Sth/NikBT as a function of log    (T/Tp) at β = ħωci/kBTp = 0, 0.1, 1, 10, 100, and 104 (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 sth; here, the solid line corresponds to the numerical data at β = 100.

Open with DEXTER
In the text
thumbnail Fig. 4

Thermal phonon contribution to the reduced heat capacity CV,lat/NikBT as a function of log (T/Tp) at β = ħωci/kBTp = 0, 1, 10, 100, and 103 (numbers near the lines). The analytical approximation in Eq. (80) (short-dashed 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
thumbnail Fig. 5

Reduced first moment of phonon frequencies (solid line) and log u1 (dashed line) as functions of β = ħωci/Tp. The dotted line shows the asymptote log (β/3).

Open with DEXTER
In the text
thumbnail Fig. 6

Reduced thermodynamic functions P/nikBT, S/NikB, CV/NikB, χρ, and χT for a fully-ionized nonmagnetic (dashed lines) and magnetized (B = 1012 K, solid lines) iron plasma at T = 107 K. The vertical dotted lines mark the densities at which (1) ; (2) TF = T; (3) ρ = ρB; (4) Γ = Γm; and (5) Tp = T.

Open with DEXTER
In the text
thumbnail 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 dot-dashed and dashed segments of the curves correspond to the domains of nonperturbative quantum effects (T < 0.5   Tp) and electron response (Z2   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
thumbnail Fig. 8

Structure of a magnetar envelope having the ground-state nuclear composition, the effective temperature 106.5 K, and magnetic field B = 1015 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 dot-dashed line. Thin vertical dot-dashed 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.