Relativistic Degenerate Electron Equation Of State

Preprint LDC-2009-001 4/16/2009 ROUGH DRAFT

ON THE DEGENERATE ELECTRON EQUATION OF STATE Lawrence D. Cloutman [email protected]

Abstract We present a set of informal notes on the equation of state appropriate to deep stellar interiors of ordinary stars. We assume the ions are treated as ideal gases, the radiation is treated in the gray one-temperature (1T) approximation, and the electrons are degenerate. We consider both the thermal and caloric equations of state for arbitrary degrees of degeneracy. Non-relativistic and relativistic formalisms are considered for a small selection of zero-dimensional problems where there is no flow and no chemistry. These cases can be helpful in validating equation of state software.

c

2009 by Lawrence D. Cloutman. All rights reserved. 1

1

Equations of State

We consider thermal and caloric equations of state appropriate to the interiors of ordinary stars. These have been incorporated into an updated version of the COYOTE computer program [1], which was used to produce the numerical results presented here. A short table of physical constants is given in the appendix. All units are CGS and the temperature T is in K. The thermal equation of state is the sum of partial pressures for each species, with the ionic species treated as ideal gases: P =

X Rρα T

Mα

α6=e

+ Pe +

aT 4 , 3

(1)

where P is the total pressure, R is the universal gas constant, ρα is the density of species α, a is the radiation energy density constant, and Mα is the molecular weight of species α. The electron pressure Pe can be calculated from the usual degenerate gas equation of state or as just another ideal gas, depending upon the thermodynamic conditions. The radiation pressure term assumes local thermodynamic equilibrium between the gas and radiation field. The caloric equation of state provides the relationship among temperature, density, internal energy, and composition and is given by ρI =

X

ρα Iα (T ) + ρe Ie + aT 4 ,

(2)

α6=e

where Iα is the ionic species specific thermal internal energy, and Ie is the electron internal energy. In the present application, we assume Iα = Cvα T =

RT (γα − 1)Mα

(α 6= e).

(3)

As in the case of the electron pressure, the electron internal energy can be calculated from the ideal gas case (3) in the nondegenerate limit. The degenerate non-relativistic electron equations of state are [2, 3] ne (η, T ) = Pe (η, T ) =

4π (2me kB T )3/2 F1/2 (η) = 5.44940154 × 1015 T 3/2 F1/2 (η), 3 h 8πkB T (2me kB T )3/2 F3/2 (η) = 0.501583989 T 5/2 F3/2 (η), 3 3h Z ∞ un Fn (η) = du, exp(u − η) + 1 0

(4) (5) (6)

and ρe Ie = 1.5Pe , 2

(7)

where h is Planck’s constant, me is the electron mass (not to be confused with the electron molecular weight Me ), kB is the Boltzmann constant, and the integration variable is u = p2 /2me kB T . The degeneracy parameter is the chemical potential in units of kB T , η = µ/kB T . Given ne and T , first F1/2 (η), then η are computed. Then F3/2 (η) and finally Pe may be computed. This electron equation of state is correct unless relativistic effects become important, which we shall discuss in the next section. The degenerate electron equation of state requires numerical evaluation of the FermiDirac integrals and their inverses. A simple but useful approximation that avoids the calculation of η is [4] i5/3

h

F3/2 (η) = F1/2 (η)

3 1 + 0.1938F1/2 (η) h

i ,

2 1 + 0.12398F1/2 (η)

(8)

which is accurate to within 0.02 percent for η < 30. Accurate tables of the Fermi-Dirac integrals are given by Cloutman [5], and Antia [6] provides accurate rational approximations to the integrals and their inverses. For problems in the ideal gas limit, use of the ideal gas approximation for electrons is computationally less complex and more efficient than the degenerate electron equation of state. For η < −3 (or F1/2 (η) < 4.3366 × 10−2 ), the difference between the two is less than one percent. Using this accuracy criterion, the ideal gas limit is valid for ne < 2.36 × 1014 T 3/2

cm−3 ,

(9)

ρe < 2.15 × 10−13 T 3/2

g/cm3 .

(10)

or

Less restrictive criteria are found commonly in the literature. We define the mean molecular weight per free electron µe by X YZ Z 1 = . µe Z AZ

(11)

The value of µe is usually between 1 and 2. Cox and Giuli [7] (p. 844) show that for η = 0, which is where degeneracy definitely becomes significant, ρ = 6.12 × 10−9 T 3/2 µe

g/cm3 .

(12)

Slightly different choices of accuracy criterion give slightly different numerical factors (for example, see Clayton [3], p. 88).

3

2

Relativistic Degeneracy

At the high densities and temperatures where relativistic effects must be included, equations (4) through (7) must be generalized. The general case is given by equations (161) and (163) of chapter X of Chandrasekhar’s classic book [2]: 1 8π Z ∞ p2 dp , h3 0 exp (E/kB T − η) + 1

ne =

8π Z ∞ Ep2 dp ρe Ie = 3 , h 0 exp (E/kB T − η) + 1 and Pe =

p3 8π Z ∞ ∂E dp, 3 3h 0 exp (E/kB T − η) + 1 ∂p

(13)

(14)

where p is the electron momentum and E is its kinetic energy. Equations (4) through (7) are derived from these general equations using the Newtonian approximation E = p2 /2me . Following Cox and Giuli [7], we note the following assumptions are made in the above equations: 1. Thermodynamic equilibrium 2. Isotropic momentum distribution 3. Fermi-Dirac distribution function dNe (p) = 8πp2 dp/{h3 [exp(E/kB T − η) + 1]} 4. Weakly interacting electrons 5. No e− -e+ pair creation To generalize to the relativistic case, we note that p=

me v [1−( v/c )2 ]1/2

=

vEt c2

(15)

for a single electron with total energy Et (rest energy me c2 plus kinetic E). This is equivalent to Et2 = p2 c2 + m2e c4 =

m e c2 1 − (v/c)2

(16)

Chadrasekhar’s equations (48) and (65) in chapter X give the kinetic energy of a single free electron as 

E = p2 c2 + m2e c4 1

1/2

− me c2

(17)

Note that Chandrasekhar uses a different convention than me for the chemical potential: His α = −η.

4

Note that this equation may be rearranged algebraically to give Chandrasekhar’s equation (160), E2 + 2Eme = p2 . (18) 2 c At low energies, this reduces to the classical Newtonian approximation. The derivative ∂E/∂p is the electron speed, p p2 ∂E = 1+ 2 2 v = ∂p me me c

!−1/2

.

(19)

This reduces to p/me and c in the classical and extreme relativistic regimes, respectively. Divine [8] recasts the degenerate electron equation of state into a form valid for all speeds, from the slow classical limit to the extreme relativistic limit. To do this, the FermiDirac integrals are generalized to the following three cases: IN 1/2 (η, β) =

∞

Z 0

u1/2 du 5 β 3 u3 2 2 1 + βu + 2β u + exp(u − η) + 1 2 2

IU 3/2 (η, β) =

Z 0

∞

∞

Z

IP 3/2 (η, β) =

0

u3/2 du 1 1 + βu exp(u − η) + 1 2 

!1/2

(20)

3/2

u3/2 du 5 β 3 u3 1 + βu + 2β 2 u2 + exp(u − η) + 1 2 2

(21) !1/2

.

(22)

The integration variable is the scaled energy of an individual electron u = E/kB T , and β = kB T /me c2 =

T . 5.93 × 109

(23)

The degenerate electron equations of state are 4π (2me kB T )3/2 IN 1/2 (η, β), 3 h

(24)

8πkB T (2me kB T )3/2 IP 3/2 (η, β), 3 3h

(25)

ne (η, T ) = Pe (η, T ) = and

4πkB T (2me kB T )3/2 IU 3/2 (η, β), (26) h3 In the nonrelativistic limit (β → 0), we recover equations (4) through (7). See also Blinnikov, ρe Ie = Ue (η, T ) =

et al. [9] for expansions appropriate to all levels of degeneracy and relativistic effects. We note that most of the literature uses a form of the generalized Fermi-Dirac integrals (20) through (22) written in terms of another generalized Fermi-Dirac integral, Fn (η, β) =

Z 0

∞

un (1 + 0.5βu)1/2 du (n > −1). exp(u − η) + 1 5

(27)

Then IN 1/2 (η, β) = F1/2 (η, β) + βF3/2 (η, β)

(28)

β F5/2 (η, β) 2 IU 3/2 (η, β) = F3/2 (η, β) + βF5/2 (η, β).

IP 3/2 (η, β) = F3/2 (η, β) +

(29) (30)

While this form is mathematically elegant, it has no clear computational advantage over Divine’s form. When are relativistic effects important? Part of the conventional wisdom has it that relativistic effects are important for densities greater than 106 − 107 g/cm3 . Clayton [3] derives this limit assuming that relativistic effects are important once the Fermi energy is twice the rest mass energy of an electron. This limit is far too lax if one needs pressure and internal energy to accuracies of a few percent or less. That is, this criterion is a sufficient condition for relativistic effects to be significant. However, it is clear from the definitions of the generalized Fermi-Dirac integrals that relativistic effects are important if β is large, regardless of the density. There is an additional limit emphasized by Mitalas [10], and a derivation of this limit may be found on page 807 of Cox and Giuli. For relativistic effects to be negligible for η > 0, βη << 1. The product βη is the chemical potential in units of the electron rest mass, which clearly is a physically meaningful relativity parameter. Figure 24.4 of Cox and Giuli (p. 847) is a useful version of the traditional ρ − T plane showing where relativistic effects are important for the electron equation of state. Mitalas [10] suggests a slightly different approximation, Fn (η, β) = Fn (η) +

β Fn+1 (η). 4

(31)

He also notes F5/2 (η) ≈ 5ηF3/2 (η)/7 for η >> 1, in which case F3/2 (η, β) ≈ F3/2 (η)(1 + 15 β η / 28). Given the ease with which the integrals in equations (20) through (22) may be evaulated and tabulated using the methods of Cloutman [5], we shall not investigate the accuracy of such approximations at the present time. The integrals obey the following recursion relation (for example, Cox and Giuli, p. 810): (n + 1)Fn (η, β) =

dFn+1 (η, β) dFn (η, β) −β dη dβ

For large β, Fn (η, β) ≈ (β/2)1/2 Fn+1/2 (η, β) for all η.

6

(k > −1)

(32)

3

Partial and Pressure Ionization, Coulomb Corrections, Pair Production

There are several effects that we have not yet included that affect the equation of state under certain conditions. We shall discuss some of these in this section. Clayton [3] (pp. 139-155) provides an elementary discussion of these issues. The first issue is partial and pressure ionization. In this study, we assumed total ionization. We argue on energetic grounds that this should be a good approximation for the conditions of deep stellar interiors where the thermal energy per particle exceeds the ionization potential. In principle, partial ionization can be included in my program by using either the equilibrium chemistry package to solve the Saha equation or the kinetic chemistry package using ionization and recombination rates. There is presently no mechanism in the code to account for pressure ionization (lowering of the continuum). Eggleton, Faulkner and Flannery [11] provide a first pass at computing variable ionization, including pressure ionization. Updates to this model are presented by Pols, et al. [12]. Weaver, Zimmerman, and Woosley [13] briefly describe an average atom model for partial and pressure ionization of a multicomponent mixture. Another physical effect is the change to pressure and internal energy due to the Coulomb interactions between charged particles in the plasma: “Coulomb interactions between charged particles provide the major non-ideal correction to the pressure at the densities and temperatures encountered in stars of around a solar mass or less, while they also crucially influence the properties of matter at high densities and low temperatures.” Pols, et al. [12], p. 964. Clayton provides an introduction to this effect, and Iben, Fujimoto, and MacDonald [14] provide a more detailed and complete model. Pols, et al. also discuss Coulomb corrections. For very high temperatures, pair creation of electron-positron pairs must be included in the equation of state. This occurs when the thermal energy of fluid particles approaches or exceeds the rest mass energy of an electron, or kB T > me c2 [15]. Timmes and Arnett [16] also provide expressions for arbitrary levels of degeneracy and relativistic effects, including electron-positron pair creation. Timmes and Swesty [17] describe a numerical method for implementing tables of such equations of state while insuring thermodynamic consistency. Researchers at LLNL created the more detailed OPAL equation of state tables [18, 19, 20]. Similar work has been reported by Mihalas, Hummer, and D¨appen (MHD) [21, 22, 23].

7

4

Non-Relativistic Numerical Examples

The COYOTE computer program was used to run several zero-dimensional test cases in order to check out the software that evaluates the equations of state described in the previous section. We shall consider four sets of physics options and two densities: one nondegenerate, one degenerate. The first example is conditions for helium burning via the triple alpha process. Here the temperature is above 108 K, which corresponds to an average thermal energy per particle of 104 ev. 2 Here is a short list of the ionization energies of ions with a single electron: 3 1. H I (Z=1): 13.59844 ev 2. He II (Z=2): 54.41778 3. O VIII (Z=8): 871.4101 4. S XVI (Z=16): 3494.1892 5. Cr XXIV (Z=24): 7894.81 6. Cu XXIX (Z=29): 11567.617 We should expect that if the thermal energy (temperature) is larger than the ionization energy, then the electron will be stripped from most atoms. In the cases I will be considering first, helium is burned to carbon and oxygen. My temperatures will always exceed 104 ev, but the ionization energy of O VIII is less than 103 ev. Therefore, the approximation of total ionization should be adequate, and that approximation will be assumed in the following examples.

4.1

Ideal Gas Case Without Radiation

We consider a mixture of fully ionized 4 He and its electrons at a temperature of 5.0 × 108 K and densities of 103 and 106 g/cm3 . Here is selected output from the program. isp

label

wt

gamma

charge

C_v

[#1]

1

e-

N

5.485798959D-04

1.666666667D+00

-1

2.273463736d+11

2

4He III

N

4.001502840D+00

1.666666667D+00

2

3.116770247d+07

3 4

12C VII 16O IX

N N

1.199670852D+01 1.599052636D+01

1.666666667D+00 1.666666667D+00

6 8

1.039598901d+07 7.799471207d+06

idealg = 1, eosform = 2.0, rad1T = 0.0 2

An average thermal energy per particle of 1 ev corresponds to a temperature of 11,605.9 K. A note on nomenclature: The roman numeral following a chemical symbol is the number of electrons stripped from the ion, plus 1. Therefore, H I is a neutral hydrogen atom and He II is singly-ionized helium. 3

8

mixture mean molecular weight = 1.3342D+00 mixture C_v = 9.3483D+07

[#2]

mixture sound speed = 2.2788D+08 cm/s i j P rho I

T

3

3

3.11677D+19 1.00027D+03 4.67415D+16 5.00000D+08

3

3

3.11677D+22 1.00027D+06 4.67415D+16 5.00000D+08

Mass Fractions i j 1

2

3

4

3

3 2.74112D-04 9.99726D-01 0.00000D+00 0.00000D+00

3

3 2.74112D-04 9.99726D-01 0.00000D+00 0.00000D+00 Hand calculation produces these results for the ρ = 1000 case:

• ρe = 0.27418601, ρHe = 9.99995926 × 102 g/cm3 • ne = ρe /me = 3.009927328 × 1026 cm−3 • Cv T = 4.67415 × 1016 erg/g using Cv at [#2] • Pe = 2.077839836 × 1019 , PHe = 1.038919184 × 1019 , P = Pe + PHe = 3.11676 × 1019 dynes/cm2 • ρe Ie = 3.116759753 × 1019 , ρHe IHe = 1.558378775 × 1019 , ρI = ρe Ie + ρHe IHe = 4.675138528 × 1019 ergs/cm3 • Ie = 1.136731868 × 1020 , IHe = 1.558385124 × 1016 , I = (ρI)/ρ = 4.673876581 × 1016 ergs/g Note that I have presented 10 digits in some of the above numbers. This runs counter to standard practice of reporting results only to the number of significant figures determined by the least accurate input number. In the present example, we are limited to 4 or 5 significant digits by the accuracy of the physical constants. The computer doesn’t care how many digits of R or Mα are significant, it carries all our double precision arithmetic out to about 16 digits. I did the hand calculations with a 10-digit scientific calculator using the full significance, just as the computer did, in an attempt to duplicate the computer output. Therefore, differences between my hand calculations and the computational results should be limited by rounding in the hand calculator or in the printout if it is has fewer than 10 digits. Note that the pressure calculations look correct, but the internal energy is off in the fourth digit. This small discrepancy appears to be due to the number of significant figures carried in the physical constants. The hand calculation used Cvα for each species computed from equation (3). These are listed in the printout above in the column labelled [#1]. The 9

enthalpy tables generated by gasnuc and used in the cfd program to compute the mixture Cv at [#2] uses a specific heat at constant presssure CP = 4.968 kcal/mol-K = 20.786 kJ/mol-K. The ρ = 106 g/cm3 case also appears to be working correctly. The pressure scales linearly with density, and I is independent of density.

4.2

Degenerate Gas Case Without Radiation

This is the same problem as in the previous subsection except that idealg has been changed from 1 to 2 (electrons treated as degenerate rather than as an ideal gas). i

j

3 3

3 3

p

rho

I

T

3.11881D+19 1.00027D+03 4.67703D+16 5.00000D+08 4.98584D+22 1.00027D+06 7.47681D+16 5.00000D+08

Mass Fractions i

j

1

2

3

4

3 3

3 2.74112D-04 9.99726D-01 0.00000D+00 0.00000D+00 3 2.74112D-04 9.99726D-01 0.00000D+00 0.00000D+00 Hand calculation produces these results for the ρ = 1000 case:

• ρe = 0.27418601, ρHe = 9.99995926 × 102 g/cm3 • ne = ρe /me = 3.009927328 × 1026 cm−3 • F1/2 (η) = 4.940287164 × 10−3 , η ≈ −5.2, F3/2 (η) = 7.417716111 × 10−3 • Pe = 2.079882898×1019 , PHe = 1.038919184×1019 , P = Pe +PHe = 3.118802082×1019 dynes/cm2 • ρe Ie = 3.119824347 × 1019 , ρHe IHe = 1.558378775 × 1019 , ρI = ρe Ie + ρHe IHe = 4.678203122 × 1019 ergs/cm3 • I = (ρI)/ρ = 4.676940348 × 1016 ergs/g The pressure is correctly computed by the code. We are at a low enough density that degeneracy effects are not important. The discrepancy in the internal energy is the same small value as in the previous case and is undoubtedly due to the same cause. For ρ = 106 , electron degeneracy is important, which can be seen by comparing P and I for the two densities. Here are the hand calculations: • ρe = 2.7418601 × 102 , ρHe = 9.99995926 × 105 g/cm3 10

• ne = ρe /me = 3.009927328 × 1029 cm−3 • F1/2 (η) = 4.940287164, η ≈ 3.55, F3/2 (η) = 14.07626571 • Pe = 3.946900081×1022 , PHe = 1.038919184×1022 , P = Pe +PHe = 4.985819265×1022 dynes/cm2 • ρe Ie = 5.920350122 × 1022 , ρHe IHe = 1.558378775 × 1022 , ρI = ρe Ie + ρHe IHe = 7.478728897 × 1022 ergs/cm3 • I = (ρI)/ρ = 7.476710185 × 1016 ergs/g

4.3

Degenerate Gas Case With Radiation

This is the same problem as in the previous subsection, except radiation pressure has been included (idealg = 2, rad1T = 0.0 changed to 1.0). i

j

p

3

3

1.88751D+20 1.00027D+03 5.19328D+17 5.00000D+08

3

3

5.00160D+22 1.00027D+06 7.52406D+16 5.00000D+08

Mass Fractions i j 1

rho

2

I

3

T

4

3

3 2.74112D-04 9.99726D-01 0.00000D+00 0.00000D+00

3

3 2.74112D-04 9.99726D-01 0.00000D+00 0.00000D+00 The electron and ion contributions are the same as for the case in the previous sub-

section. However, comparison with those results shows that the radiation dominates in the lower density case and is not very important in the higher density case. For ρ = 103 g/cm3 , • Pr = 1.575625 × 1020 , P = 1.887505208 × 1020 dynes/cm3 • ρIr = 4.726875 × 1020 , ρI = 5.194695312 × 1020 ergs/cm3 • Ir = 4.725599 × 1017 , I = 5.193293123 × 1017 ergs/g For ρ = 106 g/cm3 , • P = 5.001575515 × 1022 dynes/cm3 • ρI = 7.525997647 × 1022 ergs/cm3 • I = 7.523966176 × 1016 ergs/g 11

We get good agreement between the code and the hand calculation, so we conclude that the equations of state have been implemented correctly. Note one peculiarity: The internal energy density I decreases with density while ρI increases. This feature is due to the radiation pressure term, for which I does indeed decrease with increasing density.

4.4

Ideal Gas Case Without Radiation Revisited

The ideal gas equation of state has been extensively tested and used for combustion research over many years, but this is the first serious test of the astrophysical version. The results in the previous three subsections used the input file gasdata generated by the physics library program gasnuc. The results in this subsection used an input file generated with the gasmak library, which is intended for use under less extreme conditions of density and temperature. However, it should produce the same results as in subsection 4.1. To test this idea, we reran the case in subsection 4.1 with the code “switches” all set to typical combustion physics values, and produced a new gasdata from gasmak. I had to include the wt and gamma arrays in the input file citape and use eosform(k)=1.0. These results are very close to those in subsection 4.1. The difference is that gasnuc assumes all nuclei are totally ionized, whereas gasmak does not. gasdata input isp label

wt

htform

gamma

charge

1

e-

5.485800000D-04

0.000000000D+00

1.666667000D+00 -1

2

He++

4.000763040D+00

1.821902486D+03

1.666667000D+00

2

3 4

C O+

1.201120000D+01 1.599880000D+00

1.699772945D+02 3.730241396D+02

1.666667000D+00 1.666667000D+00

0 1

idealg = 1

rad1T = 0.0D+00

rgas = 8.31451D+07

species input k=1 wt=5.485800D-04 gamma=1.66667D+00 eosform=1.0D+00 ze=-1.0 k=2 wt=4.000763D+00 gamma=1.66667D+00 eosform=1.0D+00 ze= 2.0 k=3 wt=1.201120D+01 gamma=1.66667D+00 eosform=1.0D+00 ze= 0.0 k=4 wt=1.599880D+00 gamma=1.66667D+00 eosform=1.0D+00 ze= 1.0 nreg = nreg =

1 mmw, P, C_v = 1.3340D+00 3.1173D+19 9.3495D+07 1 gamma, sound spd = 1.6667D+00 2.2791D+08

i

j

P

rho

I

3 3

3 3.11735D+19 1.00027D+03 4.67474D+16 5.00000D+08 3 3.11735D+22 1.00027D+06 4.67474D+16 5.00000D+08 12

T

Mass Fractions i

j

1

2

3

3 3

3 2.74163D-04 9.99726D-01 0.00000D+00 0.00000D+00 3 2.74163D-04 9.99726D-01 0.00000D+00 0.00000D+00

13

4

5

Relativistic Numerical Examples

Tables of relativistic Fermi-Dirac integrals were calculated that were used to develop software for the relativistic equation of state. The numerical method is the same as in [5] and is summarized in appendix B. The tables cover the ranges −30 ≤ η ≤ 70 and 0.0 ≤ β ≤ 1.2. The software is based on bilinear interpolation, which is easy to implement and has the 3/2

advantage that it is reversible. That is, if one interpolates a table to find η for a given IN 3/2

and β, then interpolation for that η and β will produce the original value of IN to within rounding errors. The question of thermodynamic consistency, however, needs examination. It may be necessary to change to another algorithm in the future. Non-Relativistic Fermi-Dirac Integrals, beta = 0.000000D+00 eta F_1/2(eta) F_3/2(eta) F_5/2(eta) F_3/2/F_1/2 -2.0D+01 1.8266498364D-09 2.7399747556D-09 6.8499368901D-09 1.500000 -1.0D+01 4.0233994367D-05 6.0351475898D-05 1.5087929519D-04 1.500012 -5.0D+00 5.9571769052D-03 8.9463822604D-03 2.2379248359D-02 1.501782 0.0D+00 6.7809389515D-01 1.1528038371D+00 3.0825860828D+00 1.700065 5.0D+00 7.8379760573D+00 2.7802446216D+01 1.2748954491D+02 3.547146 1.0D+01 2.1344471492D+01 1.3427015996D+02 1.0346842542D+03 6.290629 2.0D+01 5.9812795370D+01 7.2656828397D+02 1.0590639177D+04 12.14737 3.0D+01 1.0969481834D+02 1.9853113777D+03 4.2929257585D+04 18.09850 7.0D+01 3.9053966669D+02 1.6419179065D+04 8.2233568906D+05 42.04228 The above table contains selected values of the non-relativistic Fermi-Dirac integrals from non-degenerate to highly degenerate conditions. The last column is the ratio F3/2 (η)/F1/2 (η), which has a value of 1.5 in the non-degenerate regime. Degeneracy effects begin to be important in the range −5 ≤ η ≤ 0. Relativistic Fermi-Dirac Integrals, beta = eta I_N^1/2 I_P^3/2 F_1/2(eta,beta) F_3/2(eta,beta) -2.0D+01 1.8300763027D-09 2.7451144550D-09 F_n 1.8273346162D-09 2.7416864914D-09 -1.0D+01 4.0309466704D-05 6.0464684858D-05 F_n 4.0249077525D-05 6.0389179238D-05 -5.0D+00 5.9683647766D-03 8.9631740393D-03 F_n 5.9594128020D-03 8.9519746268D-03 0.0D+00 6.7953557400D-01 1.1551168247D+00 F_n 6.7838199987D-01 1.1535741345D+00 5.0D+00 7.8727569873D+00 2.7898127777D+01 F_n 7.8449226902D+00 2.7834297156D+01 1.0D+01 2.1512535321D+01 1.3504700934D+02 14

1.000000D-03 (T = 5.93d+06 K) I_U^3/2 I_P / I_N F_5/2(eta,beta) I_U / I_P 2.7485424186D-09 1.50000 6.8559272181D-09 1.00125 6.0540190479D-05 1.50001 1.5101124068D-04 1.00125 8.9743734518D-03 1.50178 2.2398825018D-02 1.00125 1.1566595149D+00 1.69986 3.0853804114D+00 1.00134 2.7961958398D+01 3.54363 1.2766124123D+02 1.00229 1.3556546596D+02 6.06484

F_n 2.0D+01 F_n 3.0D+01 F_n 7.0D+01 F_n

2.1378006768D+01 6.0723318472D+01 5.9994107800D+01 1.1218582469D+02 1.1018981251D+02 4.1124248826D+02 3.9461910790D+02

1.3452855273D+02 7.3452712403D+02 7.2921067173D+02 2.0176032258D+03 1.9960121877D+03 1.7040118103D+04 1.6623380365D+04

1.0369132308D+03 7.3984357633D+02 1.0632904599D+04 2.0391942639D+03 4.3182076140D+04 1.7456855841D+04 8.3347547647D+05

1.00384 12.0963 1.00724 17.9845 1.01070 41.4357 1.02446

The above table contains selected values of the relativistic Fermi-Dirac integrals from non-degenerate to highly degenerate conditions for β = 0.001, which corresponds to a temperature of only 5.93 × 106 K. There are two kinds of lines in this table. Lines with a numerical value of η in the first column contain values of IN 1/2 , IP 3/2 , IU 3/2 , and the ratio IP 3/2 /IN 1/2 . This ratio should be 1.5 in the non-degenerate, non-relativistic (ideal gas) limit, and that is the case. Lines that begin with Fn contain the generalized Fermi-Dirac integrals F1/2 (η, β), F3/2 (η, β), F5/2 (η, β), and the ratio IU 3/2 /IP 3/2 . This ratio should be 1.0 in the non-relativistic ideal gas limit. These values should be non-relativistic to a high degree of accuracy. Indeed, even the most extreme table entry F5/2 (70, 0.001) = 8.33 × 105 is very close to the non-relativistic value of F5/2 (70) = 8.22 × 105 . As pointed out by Cox and Giuli [7] and Mitalas [10], βη is the relevant relativity parameter if η > 0, and this parameter has a value of 0.07 in this case. Relativistic corrections of a few percent are therefore the most that we should expect unless the degeneracy becomes even stronger. Relativistic Fermi-Dirac eta I_N^1/2 -2.0D+01 1.8610488070D-09 -1.0D+01 4.0991676049D-05 -5.0D+00 6.0694944118D-03 0.0D+00 6.9257111528D-01 5.0D+00 8.1882795811D+00 1.0D+01 2.3045277346D+01 2.0D+01 6.9122745126D+01 3.0D+01 1.3542777639D+02 4.0D+01 2.2204918983D+02 7.0D+01 6.1283442222D+02

Integrals, beta = 1.000000D-02 I_P^3/2 I_U^3/2 2.7915732115D-09 2.8261209177D-09 6.1488002972D-05 6.2248963800D-05 9.1149584335D-03 9.2278286230D-03 1.1760276835D+00 1.1915796198D+00 2.8765028924D+01 2.9411003478D+01 1.4211333702D+02 1.4739706693D+02 8.0756464659D+02 8.6259256375D+02 2.3166022714D+03 2.5435570133D+03 4.9718830712D+03 5.5992962195D+03 2.2989230017D+04 2.7626633947D+04

I_P/I_N 1.5000 1.5000 1.5018 1.6981 3.5130 6.1667 11.6831 17.1058 22.3909 37.5130

I_U/I_P 1.0124 1.0124 1.0124 1.0132 1.0225 1.0372 1.0681 1.0980 1.1262 1.2017

The above table contains selected values of the relativistic Fermi-Dirac integrals from non-degenerate to highly degenerate conditions for β = 0.01 (T = 5.9 × 107 K). The ratios IP 3/2 (η, β)/IN 1/2 (η, β) = 1.5 and IU 3/2 (η, β)/IP 3/2 (η, β) = 1.0124 for η = −30 already show 15

a small correction to the internal energy even at this low value of β. Both ratios show a deviation from the non-relativistic equation of state for η ≥ 0. The deviations in the pressure and internal energy are in opposite directions: Relativistic effects lower the pressure from the non-relativistic value, but raise the internal energy. The ratio IU 3/2 (η, β)/IP 3/2 (η, β) is unity in the non-relativistic limit. If one needs the equation of state to be highly accurate, say to one or two percent, then the relativistic form of the Fermi-Dirac integrals must be used for all advanced stages of thermonuclear burning in stars. That is the case for our He-burning deflagrations, which have β values ranging from 0.02 to 0.2 and η ranging from around 0 on up to (more typically) very large values. Relativistic Fermi-Dirac eta I_N^1/2 -2.0D+01 2.1099706215D-09 -1.0D+01 4.6474500680D-05 -5.0D+00 6.8822748209D-03 0.0D+00 7.9756571221D-01 5.0D+00 1.0789258068D+01 1.0D+01 3.6127876763D+01 2.0D+01 1.4571631007D+02 3.0D+01 3.5969711021D+02 4.0D+01 7.0982871156D+02 7.0D+01 2.8965620597D+03

Integrals, beta = 8.000000D-02 I_P^3/2 I_U^3/2 3.1649559333D-09 3.4573522515D-09 6.9712272321D-05 7.6152704776D-05 1.0334840467D-02 1.1290138286D-02 1.3442792541D+00 1.4761587432D+00 3.5848901929D+01 4.1468602335D+01 2.0141514315D+02 2.4938290334D+02 1.4544288748D+03 1.9962686266D+03 5.0947576926D+03 7.4855216620D+03 1.2926075918D+04 1.9916917487D+04 8.6778679108D+04 1.4526126426D+05

I_P/I_N 1.5000 1.5000 1.5017 1.6855 3.3226 5.5751 9.9812 14.1640 18.2101 29.9592

I_U/I_P 1.0924 1.0924 1.0924 1.0981 1.1568 1.2382 1.3725 1.4693 1.5408 1.6739

Relativistic Fermi-Dirac eta I_N^1/2 -2.0D+01 3.8647559200D-09 -1.0D+01 8.5125935008D-05 -5.0D+00 1.2612456657D-02 0.0D+00 1.5444667734D+00 5.0D+00 3.0920462644D+01 1.0D+01 1.4791258199D+02 2.0D+01 8.9445052736D+02 3.0D+01 2.7410012980D+03 4.0D+01 6.1875555471D+03 7.0D+01 3.1127224576D+04

Integrals, beta = 5.000000D-01 I_P^3/2 I_U^3/2 5.7971338816D-09 8.1250650143D-09 1.2768963684D-04 1.7896555141D-04 1.8934790886D-02 2.6540859491D-02 2.5360735014D+00 3.5928741185D+00 8.9053253616D+01 1.3762655541D+02 6.8542850125D+02 1.1427387425D+03 7.4406200481D+03 1.3328810905D+04 3.3019001872D+04 6.1061718833D+04 9.7670675641D+04 1.8384141223D+05 8.4106570839D+05 1.6224610601D+06

I_P/I_N 1.5000 1.5000 1.5013 1.6420 2.8801 4.6340 8.3186 12.0463 15.7850 27.0203

I_U/I_P 1.4016 1.4016 1.4017 1.4167 1.5454 1.6672 1.7914 1.8493 1.8823 1.9291

The values of β = 0.08 and 0.5 were chosen to be in the middle of the range of T encountered in He burning and near the upper limit of T found all but the most advanced stages of nuclear burning. The ratios clearly show the importance of relativistic effects at these temperatures regardless of the density (which determines the degeneracy parameter η for a given T ). Note also that the ratio IU 3/2 (η, β)/IP 3/2 (η, β) is approaching the extreme 16

relativistic limit of 2.0 as η and β increase. Note also that relativistic effects become important at lower values of η and β for the caloric equation of state than for the thermal equation of state. The question of accuracy of the numerical integrations used to make these tables is legitimate. The code is the same one used to make extensive tables for β = 0 that are accurate to at least 10 digits [5]. The numerical method is based on Simpson’s rule with Aitken extrapolation to obtain values for finite integration intervals, followed by two extrapolation techniques to obtain the value on the infinite interval. Internal checks indicate that the tables in this appendix were computed to an accuracy of 12-15 digits. This is consistent with the 16 digit accuracy of double precision artithmetic on my computer. A cross-check was made by using Gauss-Laguerre quadrature with 8, 16, 32, and 64 nodes to compute some of these integrals, an approach proposed by Kippenhahn and Thomas [24]. We obtained agreement to about 10 digits in the non-degenerate limit. However, the Gauss-Laguerre quadrature failed miserably at high degeneracy, with the values of the integrals not converging smoothly as the number of integration nodes was increased. While the source of this behavior was not investigated, there are three obvious suspects: Round-off errors, coding errors, and (most likely) the fact that the error term is proportional to the n-th order derivative of the integrand, where n is the number of nodes. If the integrand is such that derivatives grow without bound as the order increases, then the truncation error is also unbounded. My code has been checked in the limits of large and small η by comparing to analytic series expansions of the integrals and by comparing to earlier published values. Now let us consider the relativistic effects by examining the governing equations. Consider equation (4) rewritten as F1/2 (η) = 1.105 × 108

ρ −3/2 T µe

(33)

as a measure of the importance of degeneracy effects. If F1/2 (η) < 6.0 × 10−3 , or η < −5, then degeneracy is negligible (approximately 0.1 percent effect on the pressure). Relativistic effects are important when β, defined by equation (23), is much less than unity. For β = 0.01, we have T = 5.93×107 K. We saw from the table above that even at this low (by the standards of He burning) temperature, relativistic effects are present at the one percent level. Even for η ≈ −5, the effect is about 2 percent, and it gets larger the more degenerate the electrons. The density for these conditions is ρ/µe = 24.6 g/cm3 . The molecular weight per electron is a weak function of composition for H-depleted mixtures, and its value is approximately 2.0.

17

6

Conclusions

The results of this study suggest the following conclusions: 1. It is possible, with relatively elementary numerical techniques, to provide Fermi-Dirac integrals, both relativistic and non-relativistic, with an accuracy of at least one part in 1010 . 2. It is easy to generate extensive tables of the required generalized Fermi-Dirac integrals. Only a few selected values are presented here. 3. Relativistic electron degeneracy is important in all advanced phases of stellar evolution, including helium burning. 4. The degenerate electron equation of state for arbitrary degrees of degeneracy and relativity parameter is a straightforward generalization of the non-relativistic equation of state and is little more difficult to implement into a fluid dynamics code. 5. Some simple Pade approximants or similar numerical device would be convenient for computing the integrals as functions of η and β, as well as for their inverses.

18

References [1] L. D. Cloutman, COYOTE: A Computer Program for 2D Reactive Flow Simulations, Lawrence Livermore National Laboratory Report UCRL-ID-103611, 1990. [2] S. Chandrasekhar, An Introduction to the Study of Stellar Structure (University of Chicago Press, Chicago, 1939). [3] D. D. Clayton, Principles of Stellar Evolution and Nucleosynthesis (McGraw-Hill, New York, 1968). [4] R. B. Larson and P. R. Demarque, “An application of Henyey’s approach to the integration of the equations of stellar structure,” Ap. J. 140, 524 (1964). [5] L. D. Cloutman, “Numerical evaluation of the Fermi-Dirac integrals,” Ap. J. Suppl. 71, 677 (1989). [6] H. M. Antia, “Rational function approximations for Fermi-Dirac integrals,” Ap. J. Suppl. 84, 101 (1993). [7] J. P. Cox and R. T. Giuli, 1968, Principles of Stellar Structure, Vol. 2. Physical Principles (New York: Gordon and Breach); chapt. 24. [8] N. Divine, “Numerical evaluation of the degenerate equation of state,” Ap. J. 142, 1652 (1965). [9] S. I. Blinnikov, N. V. Dunina-Barkovskaya, and D. K. Nadyozhin, “Equation of state of a Fermi gas: Approximations for various degrees of relativism and degeneracy,” Ap. J. Suppl. 106, 171 (1996). [10] R. Mitalas, “Relativity parameter at high degeneracy,” Ap. J. 172, 179 (1972). [11] P. P. Eggleton, J. Faulkner, and B. P. Flannery, “An approximate equation of state for stellar material,” Astron. Astrophys. 23, 325 (1973). [12] O. R. Pols, C. A. Tout, P. P. Eggleton, and Z. Han, “Approximate input physics for stellar modelling,” MNRAS 274, 964 (1995). [13] T. A. Weaver, G. B. Zimmerman, and S. E. Woosley, “Presupernova evolution of massive stars,” Ap. J. 225, 1021 (1978). [14] I. Iben, Jr., M. Y. Fujimoto, and J. MacDonald, “Diffusion and mixing in accreting white dwarfs,” Ap. J. 388, 521 (1992). 19

[15] H.-Y. Chiu, 1968, Stellar Physics (Waltham, MA: Baisdell). [16] F. X. Timmes and D. Arnett, “The accuracy, consistency, and speed of five equations of state for stellar hydrodynamics,” Ap. J. Suppl. 125, 277 (1999). [17] F. X. Timmes and F. D. Swesty, “The accuracy, consistency, and speed of an electronpositron equation of state based on table interpolation of the Helmholtz free energy,” Ap. J. Suppl. 126, 501 (2000). [18] C. A. Iglesias, F. J. Rogers, and B. G. Wilson, “Spin-orbit interaction effects on the Rosseland mean opacity,” Ap. J. 397, 717 (1992). (OPAL) [19] F. J. Rogers and C. A. Iglesias, “Radiative atomic Rosseland mean opacity tables,” Ap. J. Suppl. 79, 507 (1992). (OPAL) [20] F. J. Rogers, 1994, in G. Chabrier and E. Schatzman, eds., Proc. IAU Colloq. 147, The Equation of State in Astrophysics. Cambridge Univ. Press, Cambridge. (OPAL) [21] D. G. Hummer and D. Mihalas, “The equation of state for stellar envelopes. I. An occupation probability formalism for the truncation of internal partition functions,” Ap. J. 331, 794 (1988). (MHD) [22] D. Mihalas, W. D¨appen, and D. G. Hummer, “The equation of state for stellar envelopes. II. Algorithm and selected results,” Ap. J. 331, 815 (1988). (MHD) [23] W. D¨appen, D. Mihalas, D. G. Hummer, and B. W. Mihalas, “The equation of state for stellar envelopes. III. Thermodynamic quantities,” Ap. J. 332, 261 (1988). (MHD) [24] R. Kippenhahn and H.-C. Thomas, “Integralapproximationen f¨ ur die Zustandsgleichung eines entarteten Gases,” Zs. f. Ap. 60, 19 (1964). [25] E. Isaacson and H. B. Keller, 1966, Analysis of Numerical Methods (New York: John Wiley). [26] H. L. Gray and T. A. Atchison, 1968, in Proc. ACM 23rd National Conf. (Princeton: Brandon/Systems Press), p. 73. [27] M. Abramowitz and I. A. Stegun, 1965, Handbook of Mathematical Functions (New York: Dover). [28] A. Sommerfeld, Z. Phys. 47, 1 (1928). [29] J. McDougall and E. C. Stoner, Phil. Trans. 237, 67 (1939).

20

A

Physical Constants • Radiation energy density constant: a = 7.563 × 10−15 erg cm−3 K−4 (astef in the code) • Atomic mass unit: 1.6605402 × 10−24 g (atmass) • Universal gas constant: R = 8.314510 × 107 erg mol−1 K−1 (rgas) • Planck’s constant: h = 6.6260755 × 10−27 erg s (hplanck) • Reduced Planck’s constant: h ¯ = 1.05457266 × 10−27 erg s (hbar) • Avogadro’s number: NA = 6.0221367 × 1023 mol−1 (avogad) • Speed of light: c = 2.99792458 × 1010 cm/s (clight) • Boltzmann’s constant: kB = 1.380658 × 10−16 erg/K (boltzk) • Electron mass: me = 9.1093897 × 10−28 g (emass) • Proton charge: e = 4.8029 × 10−10 esu (echarge) • π = 3.141592654 (pi)

21

B

Excerpts from ApJS Paper [5]

B.1

Numerical Integration Techniques

In general numerical quadrature must be used to evaluate the Fermi-Dirac integrals for conditions of partial degeneracy (η near zero), and attention must be paid to several complications. For n < 0, the integrand has a singularity at the origin. For nonintegral n, certain derivatives of the integrand have a singularity at the origin. Finally, the interval of integration is infinite. This section describes the techniques and strategies that avoid numerical difficulties associated with these features and allow one to obtain excellent accuracy with only a small computational effort. The algorithm we adopt is based on a pair of extrapolation procedures and readily produces accurate results without manual input beyond setting up the program and checking the output for consistency of the extrapolations. The difficulty with singularities in the derivatives is due to the truncation error term of the chosen quadrature rule. Consider the example of Simpson’s rule applied to the evaluation of a definite integral. Simpson’s rule is derived by fitting a parabola to three equidistant points and integrating the parabola. For m (an odd integer) equally spaced points, m ≥ 5, the three-point rule may be applied to a sequence of three-point subintervals to obtain Z

b

a

(m−1)/2 (m−1)/2 X X H f (x) dx = f (a + 2jH) f (a + [2j − 1] H) + 2 f (a) − f (b) + 4 3 j=1 j=1





(b − a) 4 (4) H f (ξ), (34) 180 where H is the integration step size, and ξ is some point in the open interval (a, b) = −

(a, a + [m − 1] H) [25]. The truncation error term is valid only if the first four derivatives of the integrand are bounded in the interval of integration. For n = 1/2 and 3/2, the first and second derivatives, respectively, are singular at the origin, and Simpson’s rule is less accurate than shown in equation (34). This is because, for n = 1/2, the slope of the integrand at the origin is infinite, and the parabola used to fit the integrand at the first three integration nodes must have a finite slope and therefore cannot make a very accurate fit. The n = 3/2 case is less accurate because the parabolic fit has similar trouble with the infinite second and higher derivatives of the finite integrand. This difficulty is not peculiar to Simpson’s rule. For example, the trapezoidal rule Z

b

a





m−2 X H  (b − a) 2 (2) f (x) dx = f (a) + f (b) + 2 f (a + jH) − H f (ξ) 2 12 j=1

(35)

requires the first two derivatives to be finite. Gaussian quadrature of order N requires that the first 2N derivatives be bounded, which is even more restrictive. 22

The difficulty with singular integrands and derivatives for half-integer values of n ≥ −1/2 may be eliminated by the simple change of variables z 2 = x applied to equation (6): Fn (η) = 2

Z

z 2n+1 dz .) 1 + exp(z 2 − η)

∞

0

(36)

This same change of variables is used for the relativistic integrals. Although this transformation can be used for integral values of n, it is not necessary. Evaluation of equation (6) or equation (36) begins with application of Simpson’s rule over the finite interval [0, t] for two values of t. If the values of t are properly chosen, then we can extrapolate to the integral over [0, ∞]. Consider the general case of numerically evaluating improper integrals of the form S(p) =

Z

∞

a

f (p, x) dx = lim St (p), t→∞

where St (p) =

Z

(37)

t

f (p, x) dx,

(38)

a

and where a is a finite constant. The limit in equation (37) is obtained by means of extrapolation procedures. This may be done by means of any of several transforms that were originally devised for doing Laplace transforms numerically. Define the B and G transforms [26] as G[S(p); t, k] = and B[S(p); t, k] =

St+k (p) − Rt (p, k) St (p) , Rt 6= 1, 1 − Rt (p, k)

(39)

Skt (p) − ρt (p, k) St (p) , ρt 6= 1, 1 − ρt (p, k)

(40)

where Rt (p, k) =

f (p, t + k) , k > 0, f (p, t)

(41)

ρt (p, k) =

k f (p, kt) , k > 1. f (p, t)

(42)

and

These transforms have the property that lim G[S(p); t, k] = lim B[S(p); t, k] = S(p).

t→∞

t→∞

(43)

Gray and Atchison [26] show that if limt→∞ Rt (p, k) 6= 0 or 1, G converges to S(p) faster than St+k (p). A similar theorem was given for the B transform. It can be shown that, if St is evaluated exactly, the G transformation is exact for f = exp(−x), and the B transformation is exact for f = x−s , s > 1 for finite t and k. Experience has shown that these 23

transforms do help the convergence in evaluating the Fermi-Dirac integrals. Because the Fermi-Dirac integrand decays exponentially, the G transformation is more accurate than the B transformation. The St (p) are evaluated by numerical integration and contain truncation and roundoff errors. Aitken extrapolation is used to improve the accuracy of the St (p) and is based on the assumption that, for sufficiently large m, the truncation error vanishes smoothly and monotonically as m increases. If we use µ, 2µ, and 4µ subintervals, where µ = m − 1, then we assume I = Iµ + c H p , I = I2µ + c (H/2)p ,

(44)

and I = I4µ + c (H/4)p , where I is the exact value of the integral, Iµ is the approximate value using µ subintervals, and c and p are free parameters. Equations (15) can be solved for I, c, and p given H and the Im : I = I4µ −

(I4µ − I2µ )2 , I4µ − 2I2µ + Iµ

(45)

p = log10 [(I − Iµ )/(I − I2µ )]/ log10 (2),

(46)

c = (I − Iµ )/H p .

(47)

and

For integrands with finite fourth derivatives, p ≈ 4 for Simpson’s rule. Because extrapolations are usually considered risky at best, this algorithm for evaluating Fermi-Dirac integrals has been extensively tested. A sample of the test results are detailed by Cloutman [5]. Included in the tests were the single known analytic expression for a Fermi-Dirac integral [28]. If we make the change of variables w = 1 + exp(x − η) and set n = 0, equation (6) becomes F0 (η) =

Z

∞

1+exp(−η)

  w − 1 ∞ dw = ln = ln[exp(η) + 1]. 1+exp(−η) w(w − 1) w

(48)

We obtained 12-digit accuracy with only 121 nodes in the Simpson quadratures.

B.2

Fermi-Dirac Integral Evaluation in Applications Programs

It would be inefficient to use the techniques of the previous section to evaluate a FermiDirac integral every time an applications program needs a value. Therefore the numerical integration program was used to generate accurate tables to be used with a separate function 24

subprogram to evaluate Fn in applications programs. Tables were published for the nonrelativistic case for n = -1/2, 1/2, 3/2, and 5/2 and for −5 ≤ η ≤ 25 in steps of 0.05 [5]. Outside that range, series expansions provide good accuracy. Cox and Giuli [7] show that for η ≤ 0, a good approximation is Fn (η) = Γ(n + 1) exp(η)

∞ X

(−1)r

r=0

exp(r η) , n > −1, (r + 1)n+1

(49)

where Γ is the gamma function. Using Γ(z + 1) = zΓ(z) and Γ(1/2) = π 1/2 (for example, Abramowitz and Stegun [27]), we specialize this series to F1/2 (η) =

∞ (−1)j+1 exp(j η) π 1/2 X , 2 j=1 j 3/2

(50)

F3/2 (η) =

∞ exp(j η) 3π 1/2 X (−1)j+1 . 4 j=1 j 5/2

(51)

and

Only five terms are used for η ≤ −5, which provide 12 digit accuracy for the worst cases, F1/2 (−5) and F3/2 (−5). For η > 25, another asymptotic form is used. For integrals of the form I(η) =

Z 0

∞

φ0 (u) du exp(u − η) + 1 ∞ X

≈ φ(η) + 2

C2j φ(2j) (η),

(52)

j=1

for large η. This series expansion has errors of order exp(−η). The C2j are given by C2j =

∞ X

(−1)i+1 i−2j = [1 − 2−2j+1 ] ζ(2j),

(53)

i=1

where ζ is the Riemann zeta function. The first five coefficients are C2 = π 2 /12, C4 = 7π 4 /720, C6 = 31π 6 /30240, C8 = 127π 8 /1209600, and C10 = 511π 10 /47900160. After evaluating the necessary derivatives, the final expansions are 







∞ n+1 X Y η n+1  Fn (η) = 1 + 2 C2r  k  η −2r  , n > 0, η >> 1. n+1 r=1 k=n−2r+2

(54)

This expression specializes to F1/2 (η) ≈

π 2 −1/2 7π 4 −5/2 31π 6 −9/2 1397π 8 −13/2 2 3/2 η + η + η + η + η 3 12 960 4608 81920

(55)

F3/2 (η) ≈

2 5/2 π 2 1/2 7π 4 −3/2 31π 6 −7/2 381π 8 −11/2 η + η − η − η − η . 5 4 960 10752 81920

(56)

and

25

For intermediate values of η, we interpolate in the tables. Hermite interpolation (for example, Isaacson and Keller [25]) is highly accurate and fits both function values and first derivatives at the table’s nodes. The Fermi-Dirac integrals obey the relation [29] dFn = n Fn−1 (η), n > 0, dη

(57)

so we need the table for F−1/2 for a routine that evaluates F1/2 and F3/2 . Cubic Hermite interpolation, which fits the function values and derivatives at two nodes, provides seven digits accuracy near η = -5 and 12 digits near η = 20. We finally adopted fifth order Hermite interpolation, which fits function values and derivatives at three nodes. The accuracy is at least ten digits near η = -5 and least 12 near η = 25.

26