Thermal evolution of old white dwarfs

This work is devoted to a description of thermodynamic properties of Coulomb crystals which are expected to form in white dwarf interiors. Effects of magnetic field, isotopic impurities, polarization of the electron background and crystal lattice type on the thermal evolution of white dwarfs are discussed. It is shown that the electron polarization could play a noticeable role in the cooling process. While other parameters in concern do not make a significant impact.


Introduction
White dwarfs are compact stars with masses comparable to that of the Sun (but not more than 1.38M ⊙ ) and with radii 100 times less than the Solar radius. Thus, in their inner layers the density (ρ) may reach 10 9 g/cm 3 . It is thought that matter in white dwarf interiors at densities 10 4 g/cm 3 consists of strongly degenerate electron gas and various fully ionized atomic nuclei. These ions are mostly represented by carbon 12 C and oxygen 16 O isotopes but depending on the formation history, inner layers may contain minor quantities of other elements such as neon and iron. Impurities of light elements like hydrogen and helium may occur in outer layers [1]. At early stages of the evolution this mixture of ions is in a liquid state (so-called Coulomb liquid). As star cools the ion plasma crystallizes forming a Coulomb crystal. A partially crystalized white dwarf has been recently observed [2].
The Coulomb crystal could be described as a system of point charges immersed into a uniform neutralizing background of electrons [3,4]. It is important to study thermodynamic properties of Coulomb crystals: the heat capacity of white dwarfs determines the cooling rate and allows one to derive the age of these stars from observations of their surface temperatures [1,3]. Typically, thermodynamic properties of a one component crystal have been considered analytically in the harmonic lattice model or by other numerical methods (e.g., [5,6,7,8,9,10]).
For the first time the problem of Coulomb crystals in white dwarf envelopes was taken into account in 1960s (e.g., [11,12]) and in following works was widely developed (e.g., [1,13,14]). Most authors were interested in the chemical composition of the interior of these stars and processes in its atmosphere. Apparently precise calculations of the heat capacity of Coulomb crystals have never been used to study the thermal evolution of white dwarfs. In this work we try to address this issue and aim to investigate how the accuracy of approximations, magnetic field, isotopic impurities, polarization of the electron background and the type of the crystal lattice affect the heat capacity of the star envelope and, consequently, the cooling rate of these stars.

The heat capacity of Coulomb crystals
The potential energy of a one component Coulomb crystal without polarization of the electron background is given by where Z is the charge number of ions, while N , n e ≡ Zn and n are the total number of ions, the electron and the ion number densities, respectively; r i is the radius vector of the i -th ion. We expand the potential energy in powers of ion displacements u i from their equilibrium positions X i (so that r i = X i + u i ) and retain terms up to the second order (harmonic approximation) The zero-order term in this series is the electrostatic energy, U 0 . The second order term determines the phonon spectrum of the crystal ω ν (k). At any fixed wavevector k it could be calculated by solving the dispersion equation: det{D αβ pp ′ (k) − ω 2 ν (k)δ αβ δ pp ′ } = 0, where indices α, β denote Cartesian components, indices p, p ′ enumerate ions in the elementary cell, ν enumerates oscillation modes at given k (ν = 1, . . . , 3N cell ), N cell is the number of ions in the primitive cell. The expression for the dynamic matrix D αβ pp ′ (k) can be found in [15]. Any phonon thermodynamic function can be calculated based on the information about the phonon spectrum. For instance, the heat capacity is C = N ∑ kν w 2 (4 sinh 2 (w/2)) −1 , where w ≡ ω ν (k)/T , summation is over all ν and k in the first Brillouin zone. It is well known that C ≈ 3N at high temperatures (T > T p , where T p ≡ ω p is the ion plasma temperature, ω p ≡ √ 4πnZ 2 e 2 /M is the ion plasma frequency, M is the mass of ions). At low temperatures Among all lattices the body-centered cubic (bcc) lattice is known as the lattice with the lowest Helmholtz free energy and it is thought that this lattice forms in white dwarf cores, but the difference between the energies of different lattices is very small [16]. If the electron polarization is taken into account the face-centered cubic (fcc) lattice can become more energetically favorable than the bcc lattice [17]. In the present paper we study bcc, fcc and hexagonal-close packed (hcp) lattices. The last one is the lattice with the highest packing efficiency (along with fcc lattice).
In Fig. 1A the ratios of heat capacities of the bcc, fcc, and hcp lattices are plotted as functions of t ≡ T /T p . The difference between lattices at low and intermediate temperatures is noticeable. In particular, in the quantum limit C hcp /C bcc ≈ 1.190 and C hcp /C fcc ≈ 1.089. It is also clear that the heat capacity of the hcp lattice approaches the Debye limit at much lower temperatures than those of the bcc and fcc lattices. At intermediate temperatures all ratios have minima. For the hcp to fcc ratio the minimum is equal to 0.9 at t ≈ 0.018 while C hcp /C bcc ≈ 0.75 at t ≈ 0.025 [15]. Notice that the heat capacity of more complex lattices can differ considerably from the heat capacity of the bcc lattice but it happens only at T ≪ T p . In addition we compared the previous fit of the heat capacity of the Coulomb crystal C fit from the Monte Carlo simulations [1] (it isn't connected with any particular lattice type) with our precise calculations. The ratio C fit /C bcc is plotted in Fig. 1A by dashed line. It lies between 0.86 and 1.12 and at low temperatures equal 0.977.
The phonon spectrum of the magnetized Coulomb crystal can be found from the dispersion equation spectrum and thermodynamics of such systems were widely investigated in [18] and [19]. The heat capacity of the bcc lattice in the external magnetic field is plotted in Fig. 1B for   Previously in this paper we thought that the electron background is uniform but it is not quite correct. The electron density changes in response to the presence of ions. The simplest way to take it into account is to use the linear response formalism and the static longitudinal dielectric function of electrons ε(k) ≡ 1 + κ 2 TF /k 2 , where κ 2 TF ≡ 4πe 2 ∂n e /∂µ e , µ e is the electron chemical potential (e.g., [5,17] and references therein). In Fig. 2A we plot for the bcc lattice the ratio (C(t, κ TF a) − C(t, 0))/C(t, 0) as function of t for different κ TF a, where a = (4πn/3) −1/3 is the ion-sphere radius and C(t, 0) is the heat capacity of the bcc lattice with the uniform background. This ratio is always greater than 0. At high temperatures it decreases as t −2 . At low temperatures it is independent of t and C(t, κ TF a)/C(t, 0) − 1 ∝ (κ TF a) 2 at κ TF a ≪ 1.
The linear response formalism assumes that κ TF a < 1 which is appropriate for the white dwarfs interior. For example, in fully ionized 12 C crystal at density ρ = 10 6 g cm −3 κ TF a ≈ 0.425, while κ TF a = 1 at ρ ≈ 3000 g cm −3 . At κ TF a = 1, C(t, κ TF a)/C(t, 0) does not exceed 1.48, at κ TF a = 0.5 it is below 1.105. Another factor influencing thermodynamic properties of white dwarfs' cores is the presence of impurities and dislocations in the crystal structure. If all impurities are isotopic and substitutional, the heat capacity of the crystal can be calculated by the perturbation theory of disordered crystal spectra [20]. The ratio of an impurity contribution to the crystal heat capacity (∆c is the heat capacity change per one impurity ion) over the perfect crystal heat capacity (c is the perfect crystal heat capacity per one ion) decays as t −2 in the classic regime as seen in Fig.  2B [21]. At low temperatures, this ratio tends to a constant 1.5ϵ, where ϵ ≡ (M ′ − M )/M ≥ −1, M ′ is the impurity mass and M is the mass of "basic" ions. If ϵ ≫ 1 the maximum impurity effect takes place at temperatures 10 −2 T p < T < 10 −1 T p . For impurities with M ′ = 11M maximum ∆c/c ≈ 40 reaches at t ≈ 0.0124T p . Hence, if the concentration of such impurities is 100 times less than the concentration of "basic" ions the total heat capacity will be amplified by 40%. On the other hand formation of such massive impurities in white dwarf envelopes is unlikely. Contribution of more realistic impurities is less noticeable. When the mass of impurities is not considerably different from "basic" ions the heat capacity changes don't exceed a few per cent. For instance, at low temperatures even 10% of impurities with M ′ = 0.9M diminish the total crystal heat capacity only by 1.5%.

Cooling of white dwarfs
Now we can estimate how isotopic impurities, polarization of the electron background and the type of the crystal lattice affect the thermal evolution of old white dwarfs. Consider a simple model of an isothermal white dwarf with completely ionized crystallized carbon core and helium surface. All effects associated with ion mixtures and crystallization process are discarded (we assume that envelopes are fully crystalized). Crystallized white dwarfs have virtually no own sources of thermonuclear energy and shine exclusively due to cooling. Consequently, their luminosity and thermal evolution are determined by thermodynamic properties of envelopes: −L s = C(T )dT /dt, where C(T ) is the total heat capacity of the star. At high temperatures heat capacities of all lattices tend to 3N and distinctions between them can be neglected. At low and intermediate temperatures the heat capacity of crystals can significantly differ from the heat capacity of the bcc lattice which is taken as a basis. For instance, C hcp /C bcc and C fcc /C bcc lay between 0.75 and 1.2. The first ratio was fitted in [15], while C bcc in [16]. Using these equations we have calculated the thermal evolution of white dwarfs (masses 1.2M ⊙ and 0.6M ⊙ ) if their cores are modeled by crystals with the bcc lattice (thick solid line in Fig. 3) or by crystals with the hcp lattice (dashed line in Fig. 3). The massive white dwarfs (∼ 1.2M ⊙ ) crystalize at age ∼ 10 7 years. At an age ∼ 10 9 years crystals become quantum and the cooling rate of stars increases (it could be see in the left part of Fig.  3). Less massive stars crystallize at age > 10 7 years. The increase of the heat capacity leads to a decrease of the cooling speed and vice versa. Since T > 0.0068T p C bcc > C hcp and at T < 0.0068T p C bcc < C hcp the total difference between the bcc and hcp lattices is relatively small and does not exceed a few per cent (∼ 3 % at an age of 10 10 years for the massive white dwarf). Also we have used the old fit of the heat capacity to calculate the thermal evolution of stars [1]. Results are plotted by thin solid lines in Fig. 3. They lie even closer to new calculations for the bcc lattice than curves for the hcp lattice.
As shown previously the heat capacity can increase by about fifty present at low temperatures and densities due to the electron polarization. For simplicity we will not take into account that κ TF a depends on ρ. The maximum possible variation of the heat capacity of the bcc lattice can be described like a Heaviside function with the discontinuity at t ∼ 0.5 (it is equal 1.5 at t < 0.5 and 1 at t > 0.5). For this situation the cooling curves are plotted by dot-and-dash lines. The opposite case is described as a Heaviside function equal 1/1.5 at t < 0.5 and 1 at t > 0.5 and illustrated by dotted lines. It is not connected with any physical parameters but fully covers any possible isotopic impurities. All other cooling lines (due to the uncertainty in the lattice type or other parameters in consideration) will lie in the shaded (green) region. Hence the surface temperature of a 10 10 years old massive white dwarf can vary less than 1.2 times. It is a very optimistic range. Realistic isotopic impurities can change the surface temperature by less than 1 percent. Uncertainties in the heat capacity values will affect the thermal evolution of less massive stars (0.6M ⊙ ) at an age of ∼ 10 11 years which is greater than the age of the Universe.
For a more realistic analysis, one should take into account that white dwarfs may consist of a number of different chemical elements. Also we have't discussed the process of crystallization of stars and how it depends on the lattice type. We plan to address these problems in the future work.

Summary
We have studied heat capacity of different Coulomb crystals. The effects of the magnetic field, isotopic impurities, polarization of the electron background and the type of the crystal lattice are discussed. At high temperatures, the difference between lattices is not significant. At low temperatures, the heat capacity of some crystals may differ from the heat capacity of the bcc lattice with uniform background by up to 50%. The magnetic field does not play any role in the thermal evolution of white dwarfs while isotopic impurities and the type of the crystal lattice change the cooling rate by a few per cent only. The polarization of the electron background play more noticeable role in the cooling process of massive white dwarfs at age of 10 9.5 years or more and should be taken into account. It was also shown that the difference between the previously used heat capacity model and more precise calculations does not noticeably affect the thermal evolution of white dwarfs.