Fermi operator expansion method for nuclei and inhomogeneous matter with a nuclear energy density functional

概要

Fermi operator expansion method for nuclei and inhomogeneous matter with nuclear
energy density functional
Takashi Nakatsukasa1, 2, 3
1

arXiv:2211.09448v1 [nucl-th] 17 Nov 2022

2

Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577, Japan
Faculty of Pure and Applied Sciences, University of Tsukuba, Tsukuba 305-8571, Japan
3
RIKEN Nishina Center, Wako 351-0198, Japan
(Dated: November 18, 2022)

Background: The nuclear energy density functional method at finite temperature is a useful tool
for studies of nuclear structure at high excitation, and also for researches of nuclear matter involved
in explosive stellar phenomena and neutron stars. However, its unrestricted calculation requires large
computational costs for the three-dimensional coordinate-space solvers, especially for the Hamiltonian matrix diagonalization and (or) the Gram-Schmidt orthonormalization of the single-particle
wave functions.
Purpose: We test numerical performance of a numerical method, that requires neither the diagonalization nor the Gram-Schmidt orthonormalization, for finite nuclei and inhomogeneous nuclear
matter. We examine its advantageous features in future applications.
Methods: The Fermi operator expansion method, which approximates the Fermi-Dirac distribution in terms of the Chebyshev polynomials, is used to construct the one-body density matrix for
the energy density functional calculations at finite temperature. The modified Broyden’s mixing
method is adopted for the self-consistent iteration process.
Results: The method is applied to isolated finite N = Z nuclei and to non-uniform symmetric
nuclear matter at finite temperature, which turns out be very effective with the three-dimensional
coordinate-space representation, especially at high temperature. The liquid-gas transition is clearly
observed in the calculations.
Conclusions: The Fermi operator expansion method is a useful tool for studies of various nuclear
phases at finite temperature with the energy density functional calculations. The method is suitable for massively parallel computing with distributed memory. Furthermore, when the space size
is large, the calculation may benefit from its order-N scaling property.

I.

INTRODUCTION

It is of significant importance to calculate nuclear matter in a variety of phases with different temperature,
utilized in simulation studies of supernovae and neutron
stars. The nuclear energy density functional method at
finite temperature [1, 2] is a desirable choice for studying the inhomogeneous neutron-star matter in outer and
inner crusts. Especially, near the boundary between the
inner crust and the core, various exotic phases, “nuclear
pasta”, are expected to appear.
In order to properly treat thermally dripped nucleons and to study the transition from inhomogeneous to
uniform nuclear matter, the coordinate-space representation is preferable. Furthermore, to find a new exotic structure at finite temperature, it is desired to perform the calculation without assuming any spatial symmetry of the configuration, using the three-dimensional
(3D) coordinate-space representation. Since the 3D
coordinate-space solution is computationally demanding,
most of the finite-temperature mean-field calculations for
nuclei either adopt the harmonic-oscillator-basis (shellmodel-basis) representation [3–6], or are restricted to the
spherical systems1 [8, 9].

1

An exception can be found in a paper by Newton and Stone [7] in

To reduce the computational cost, the finitetemperature Thomas-Fermi approximation has been often adopted [10–13]. The molecular dynamics simulations can be performed with even smaller computational
time, thus, they have been extensively utilized with larger
simulation volumes [14–16]. A major drawback of these
semiclassical approximation is a lack of shell effects. In
the inner crust, the shell effects play a role not only in
protons but also in the band effect for unbound neutrons
scattered by the periodic potentials [17, 18]. An alternative quantum approach to the inner crust is to use the
Wigner-Seitz approximation, pioneered by Negele and
Vautherin [19]. The structure is optimized in a WignerSeitz sphere of radius Rws . using different boundary conditions depending on the parity of the single-particle orbitals. For the inner crust, this trick for the boundary
condition produces a roughly constant neutron density
at the spherical boundary r = Rws . However, some spurious density fluctuations still remain near the edge of
the boundary. Furthermore, the numerical results suffer from ambiguity caused by the choice of the boundary
conditions [20]. It should be noted that a combination of
the Thomas-Fermi and Wigner-Seitz approximations is

which they solve the finite-temperature Hartree-Fock equations
with the BCS treatment on the pairing correlation in the 3D
coordinate space. To reduce the computational time, the states
with the occupation smaller than 10−6 are neglected.

2
frequently used for non-uniform matter at finite temperature based on the microscopic results for the uniform
matter [13, 21, 22]. Another disadvantage of these approximations is that one loses information on the transport properties which may be crucial for understanding
dynamics of neutrons in the inner crust of neutron stars
[23].
In this paper, we perform a feasibility study for the
fully quantum energy density functional (mean-field) calculations without the Wigner-Seitz approximation for
non-uniform nuclear matter at finite temperature. A conventional solution of the finite-temperature mean-field
theory can be summarized as follows: (1) Construct the
mean-field Hamiltonian H which depends on one-body
densities. (2) Diagonalize the Hamiltonian to obtain the
eigenvalues and the eigenvectors, H|ii = i |ii. (3) Calculate the densities, then, go back to (1) to reach the
self-consistency. In the step (3), the Fermi-Dirac distribution
P function f (x) is used to calculate the densities,
ρ = i f (i )|iihi|. The truncation with respect to the
eigenvector |ii may be possible at low temperature, while,
at high temperature, we need to compute all the eigenvalues and eigenvectors. Since this diagonalization is needed
every iteration, it requires a large amount of numerical
resources.
Recently, the shifted Krylov method for the HartreeFock-Bogoliubov (HFB) theory has been proposed [24].
Then, it was extended to the finite-temperature HFB
theory [25]. The method uses the shifted Krylov subspace method for solution of linear algebraic equation,
(z − H)G(z) = 0, where H is the HFB Hamiltonian and
G(z) is the Green’s function. The densities are obtained
from the Green’s function G(z) integrated over complex
energy z. Thus, the diagonalization procedure is unnecessary. This feature is favorable for large systems since
the matrix diagonalization needs the operation of O(N 3 )
where N is the dimension of the matrix. It is shown to
be numerically feasible and efficient in the parallel computation [24, 25]. However, its performance depends on
the required number of iterations of the shifted Krylov
algorithm whose convergence is not guaranteed.
Purposes of the present paper are to study an alternative method for the finite-temperature mean-field calculation, and to examine its performance for nuclear systems.
The methodology is known as the Fermi operator expansion (FOE) method in the condensed matter physics
[26, 27]. It is also known as one of the order-N (O(N ))
method [28], thus, the number of computational operations linearly scales with respect to either the particle
number or the dimension of the one-particle space. In the
O(N ) methods, the “nearsightedness” of many electron
systems play a crucial role [29]. Since the nearsightedness is due to destructive interference effect in quantum
mechanical many-particle systems, we expect that it is
applicable to nuclear systems as well. However, since the
size of a nucleus is roughly ten femtometer at most, the
nearsightedness principle has been assumed not so beneficial in practice. The situation may be different for hot

nuclei and macroscopic neutron-star matters. It is worth
examining the O(N ) methods for calculations of nuclei at
finite temperature and inhomogeneous nuclear matter.
The paper is organized as follows: The finitetemperature mean-field theory is recapitulated in
Sec. II A. In Secs. II B and II C, the Fermi operator expansion method is summarized. In Sec. II D, we propose
an efficient method of computing the entropy without
calculating single-particle energies. The nearsightedness
and the O(N ) method are briefly reviewed in Sec. II E.
Details of the numerical calculations, examination of the
validity of the Chebyshev polynomial expansion, and numerical results for finite nuclei and non-uniform matter
are shown in Sec. III. Concluding remarks are given in
Sec. IV.
II.

THEORY AND NUMERICAL METHODS
A.

Mean-field theory at finite temperature

We recapitulate here the Hartree-Fock (HF) theory at
finite temperature [1]. The partition function and the
statistical density matrix at
β −1 = kB T
h the temperature
i
ˆ0
−β H
ˆ = e−β Hˆ 0 /Z,
are in the form, Z = tr e
, and D
ˆ0 ≡ H
ˆ − µN
ˆ with the one-body HF
respectively, where H
ˆ the particle number operator N
ˆ , and
Hamiltonian H,
the chemical potential µ. The one-body density matrix
is given as
h
i X
ˆ c† cˆi =
ρij = tr Dˆ
hi|αifβµ (α )hα|ji
(1)
j
α

where fβµ is the Fermi-Dirac function fβµ (x) = {1 +
eβ(x−µ) }−1 . Here, the subscripts i and j denote the
indices for an arbitrary single-particle basis, while α
for the single-particle states to diagonalize H (and H 0 ),
ˆ
H|αi
= α |αi. cˆi (ˆ
c†i ) is an annihilation (creation) operator for a particle at the state |ii. Since the HF Hamilˆ
tonian H[ρ]
is a functional of the one-body density, the
states |αi and energies α depend on ρij . Thus, Eq. (1)
should be iteratively calculated until the self-consistency
is achieved. It is straightforward to extend the theory to
the HFB theory at finite temperature [1, 25].
It should be worth mentioning that the finitetemperature HF theory can be derived by the principle
of maximum entropy, with an assumption that the parˆ
tition function is given in a form, Z = tr[e−β K ] with
ˆ
a one-body operator K. Constraining the energy and
the particle number with Lagrange multipliers (associated with β and µ), it is equivalent to the minimization
of the thermodynamic potential [1].
J = E − T S − µN
h
i
h
i
ˆ ln D
ˆ − µtr D
ˆN
ˆ
= E[ρ] + kB T tr D
h
i
ˆK
ˆ − kB T ln Z,
= E 0 [ρ] − tr D

(2)
(3)
(4)

3
where E[ρ] is the energy density functional and E 0 ≡
E − µN . Taking the variation with respect to the oneˆ it leads to
body operator δ K,
h
i
h
i
δE 0
ˆ K
ˆ − tr Kδ
ˆ D
ˆ − δ ln Z (5)
· δρ − tr Dδ
δJ =
δρ
β
h
i
h
i
0
ˆ δD
ˆ − tr Kδ
ˆ D
ˆ ,
= tr H
(6)
ˆ 0 ≡ P H 0 cˆ† cˆj with H 0 = δE 0 /δρji . Therewhere H
ij i
ij
ij
ˆ =H
ˆ 0.
fore, the principle of maximum entropy gives K
B.

Fermi operator expansion method

According to Eq. (1), the one-body density can be
calculated by diagonalizing H 0 to obtain the eigenstates
and the eigenenergies, |αi and α . However, since we
need to perform the diagonalization every iteration until
the self-consistency is achieved, it is prohibitively difficult
for large systems. In order to reduce the computational
cost, we must avoid the matrix diagonalization which numerically costs O(N 3 ). In this paper, we explore one of
such approaches, the Fermi operator expansion (FOE)
method.
The idea of the FOE method can be easily understood
by rewriting Eq. (1) as ρij = hi|ˆ
ρ|ji, where
X
ˆ
ρˆ =
|αifβµ (α )hα| = fβµ (H).
(7)
α

Thus, the one-body density is nothing but the FermiDirac function whose argument is replaced by the Hamiltonian. In addition, the FOE is based on the polynomial
approximation of the Fermi-Dirac distribution function.
fβµ (x) ≈

M
X

an Tn (x),

(8)

n=0

where Tn (x) is a polynomial function of the n-th degree,
and the summation is truncated at the maximum degree
M . The polynomial approximation should be better at
large T , because the Fermi-Dirac function is smoother at
higher temperature. In contrast, at the zero temperature
limit, the function becomes the Heaviside step function
for which the approximation is not so precise. Nevertheless, in case that there is a gap ∆E at the Fermi surface,
such as the shell gap and the pairing gap, the results of
the finite temperature calculation with kB T  ∆E is
practically identical to the one at zero temperature.
Inserting Eqs. (7) and (8) into ρij = hi|ˆ
ρ|ji, we have
ˆ
ρij = hi|fβµ (H)|ji
=

M
X

an hi|jn i,

(9)

n=0

ˆ
where |jn i ≡ Tn (H)|ji.
If the polynomial function Tn (x)
is simply given by Tn (x) = xn , the state |jn i can be
calculated starting from |j0 i ≡ |ji as
ˆ n−1 i,
|jn i = H|j

n = 0, · · · , M.

(10)

ˆ M times, the
Thus, multiplying the basis state |ji by H
one-body density ρij can be constructed. This is the
basic idea of the FOE method.
In practice, the simple choice of Tn (x) = xn often leads
to a numerical instability, because the functions xn are
diverging function at |x| > 1 for large n. In order to avoid
this numerical problem, a careful choice of the polynomial
functions is required for Tn (x).
C.

Chebyshev polynomials

In the present work, we adopt the Chebyshev polynomials for Tn (x). The Chebyshev polynomials of the first
kind are given by Tn (x) = cos nt with x = cos t, thus,
both x and Tn (x) are bound between −1 and √
1. They
are orthogonal with respect to the weight of 1/ 1 − x2 .
Z 1
dx
Tn (x)Tm (x) √
= Nn δnm ,
(11)
1 − x2
−1
with the normalization constants N0 = π and Nn = π/2
(n 6= 0).
First, we should change the energy scale by transformˆ
c
ˆ into H
ˆ ≡ H−e
where ec ≡ (emax + emin )/2 and
ing H
er
ˆ satisfy
er ≡ (emax − emin )/2. When the eigenvalues of H
ˆ
emin ≤ eα ≤ emax in the adopted model space, those of H
are in the interval [−1, 1]. Instead of expanding fβµ (x)
as Eq. (8), we expand a scaled Fermi-Dirac function f˜(x)
as
M

a0 X
an Tn (x),
+
f˜(x) ≡ fβµ (er x + ec ) ≈
2
n=1
where the coefficients an are given by
Z
2 1
dx
an =
Tn (x)f˜(x) √
.
π −1
1 − x2

(12)

(13)

It is worth noting that f˜(x) and an depend on both β
and µ, for which we omit these subscripts for simplicity.
Instead of Eq. (10), the recursive relations of the
Chebyshev polynomials,
Tn+1 (x) = 2xTn (x) − Tn−1 (x),

n ≥ 1,

(14)

ˆ
lead to recursive formula for |jn i ≡ Tn (H)|ji,
ˆ n i − |jn−1 i,
|jn+1 i = 2H|j

n = 1, · · · , M − 1.

(15)

ˆ
Starting with |j0 i = |ji and |j1 i = H|ji,
all the states
|jn i up to n = M are obtained, then, the one-body density is calculated as
ˆ
ρij = hi|f˜(H)|ji
=

M
X
a0
hi|ji +
an hi|jn i.
2
n=1

(16)

Let us summarize the numerical procedure to reach
the self-consistent solution of the HF problem at a given
temperature T .

4
0. The maximum and minimum energies, emax and
emin , are determined according to a problem of interest. See Sec. III A for those values. The initial
density distribution ρ(r) and the initial chemical
potential µ are given by hand. For a given T , calculate the coefficients an (n = 0, · · · , M ) according
ˆ
to Eq. (13). Set up the initial Hamiltonian H.
1. Calculate |jn i (n = 0, · · · , M ) according to Eq.
(15).
2. Construct the one-body density according to Eq.
(16). Adjust the chemical potential µ if necessary.
ˆ
3. Construct the HF Hamiltonian H[ρ]
using the calculated density ρij .
4. Check the self-consistency between the density and
the Hamiltonian. If it is self-consistent, end the
iteration. Otherwise, go to Step 1 and iterate the
procedure.
In the present formulation, the function fβµ (x) depends on T and µ. Thus, when we change the chemical potential µ, we have to recalculate the coefficients
an in Eq. (12). You may think that it is better to expand the function (1 + eβx )−1 instead of {1 + eβ(x−µ) }−1
ˆ − µ instead of H.
ˆ However, in this case,
and to use H
we need reevaluate the states |jn i using the recursion
relation (15), because |jn i depend on T and µ. Since
the calculation of |jn i requires the major portion of the
computation, we adopt the expansion of fβµ (x). If we
adjust the chemical potential in step 2 of every iteration
to fix the particle number (average density), only the coefficients an in Eq. (13) need to be recalculated. The
additional computation is negligibly small compared to
the calculation of |jn i.

D.

Calculation of entropy

It is of significant importance to calculate the entropy
of systems at finite temperature. The calculation of the
free energy requires the evaluation of the entropy as well.
For the product wave functions, the entropy S is given
by
X
S = −kB
[fα ln fα + (1 − fα ) ln (1 − fα )] ,
(17)

with the polynomial expansion as
M0

b0 X
g˜(x) ≈
+
bn Tn (x),
2
n=1

analogous to Eq. (12). The coefficients bn are determined
in the same manner as Eq. (13). Then, the entropy can
be calculated as
h
i
ˆ
S = kB tr g˜(H)
(20)


0
M
X b0 X
 +
≈
bn hj|jn i .
(21)
2
n=1
j
ˆ
Since the states |jn i = Tn (H)|ji
are calculated in Eq.
(15) in order to construct the density, almost no extra
cost is needed for evaluation of the entropy provided that
M 0 ≤ M . In fact, we find that the condition M 0 ≤ M is
well satisfied in practice (See Sec. III B). At small temperature, g˜(x) has a sharp peak at x = (µ − ec )/er , which
demands large value of M 0 . However, in this case, M
must be also large, because f˜(x) also produces a sharp
transition from 1 to 0. At T = 0, f˜(x) becomes a discontinuous Heaviside function, f˜(x) = θ(x), while g˜(x) is a
constant function, g˜(x) = 0.
In the FOE method, we end up the vectors, |jn i, which
contains the information of the Chebyshev polynomials
ˆ
of the Hamiltonian, |jn i = Tn (H)|ji.
Therefore, quantities that are continuous functions of the single-particle
energies, including the density and the entropy, can be
evaluated in principle from |jn i. The number of the vectors |jn i is N × M , where N is the dimension of the
single-particle space (system size). Since M is inversely
proportional to the temperature T as Eq. (22) below,
the FOE is more efficient at higher temperature.

E.

Nearsightedness and order-N method

The FOE method is regarded as one of the linear
system-size scaling method, namely, the order-N (O(N ))
method. According to Ref. [30], the degrees of polynomials necessary for an accuracy of 10−D (D > 1) is estimated as
M=

α

where fα = fβµ (α ). In order to calculate this, normally we need all the eigenvalues of the Hamiltonian ,
α , which requires an additional computation, namely
the diagonalization of the Hamiltonian. It demands a
large computational cost of O(N 3 ).
In this paper, we propose another manner to approximate a function
n
o n
o
g˜(x) ≡ −f˜(x) ln f˜(x) − 1 − f˜(x) ln 1 − f˜(x) , (18)

(19)

2
(D − 1)er β.
3

(22)

Assuming a gaussian basis functions of range σ centered
at mesh points whose spacing comparable to σ, the maˆ at a large separation have the same
trix elements for H
ˆ M are
width σ, and √
the long-range matrix elements for H
estimated as M σ [30]. Therefore, the range of the density matrix of Eq. (16) is approximately given as
rN ∼

√

r
Mσ ∼

~2
(D − 1)β,
3m

(23)

5
where we use Eq. (22) and er ∼ emax ∼ ~2 σ −2 /(2m) at
a small value of σ. In other words, the density matrices
ρ(r, r0 ) are localized, namely, ρ(r, r0 ) ≈ 0 at |r − r0 | >
rN . It becomes more “nearsighted” (rN → 0) for higher
temperature β → 0.
At the zero temperature limit, rN can stay finite if
there is a gap δe at the Fermi surface, although Eq. (23)
diverges. Taking the chemical potential µ as the mid
value of the gap, the condition that discrepancies between
the Heaviside function and the Fermi-Dirac function are
smaller than 10−D except for the gap interval leads to
β > (2 ln 10)

D
.
δe

(24)

For instance, when we have a shell gap of δe = 2 MeV
at the Fermi surface and require the accuracy of D = 3
(error smaller than 10−3 ), the calculation at T = 100
keV is practically identical to that at T = 0.
The nearsightedness of the density enables us to perform the O(N ) calculation. The calculation of |jn i
(n = 1, · · · , M ) in Eq. (15) can be performed in a truncated space whose dimension does not depend on the
system size. Since the nonlocal (off-diagonal) densities
ρij with Rij > rN vanish, the matrix-vector product
in Eq. (15) can be performed in the restricted active
subspace. Here, Rij mean that spatial distance between
two basis states |ii and |ji. For the coordinate-space
basis, they are trivially Rij = |ri − rj |. For the singlecenter harmonic-oscillator basis, which are common and
efficient in calculation of finite nuclei, it is difficult to find
a pair of states with Rij > rN . Thus, the applicability of
the O(N ) method also relies on the choice of the basis.
Before finishing this section, we emphasize the following advantageous features of the method in numerical
computation. First of all, in order to construct the onebody density ρij , only the matrix-vector product, operation of the Hamiltonian on a state, is necessary in Eq.
(15). Although the self-consistency between the density
ρ and the Hamiltonian H[ρ] requires the iteration, neither the matrix diagonalization nor the linear algebraic
equations are involved to achieve the mean-field solution
at the temperature T . Second, the calculations of ρij in
Eqs. (15) and (16) for different j can be independently
performed. It is suitable for massively parallel computing for large systems. Last, but not the least, the method
may receive benefits from its nearsightedness, and the
computational cost could linearly scale with the system
size (Sec III E).
III.
A.

NUMERICAL RESULTS

Energy density functional and numerical details

We use the BKN energy density functional [31], which
is a functional of the isoscalar kinetic and local densities
and assumes the spin-isospin symmetry without the spinorbit interaction. The pairing correlation is neglected.

Since the BKN functional is not suitable for description
of the neutron-rich matter, we study only the symmetric
nuclear matter and finite nuclei with N = Z. Nevertheless, it serves for the purposes of the present paper,
namely, to examine applicability and performance of the
FOE method for the mean-field (energy-density) calculation for nuclei and nuclear matter at finite temperature.
We adopt the 3D Cartesian grid representation [32] of
the square box with periodic boundary condition. The
3D grid size is set to be h3 = (1.0 fm)3 . The differentiation is evaluated with the nine-point finite difference.
For calculation of isolated finite nuclei, the center-of-mass
correction is taken into account by modifying the nucleon’s mass as m → m ∗ A/(A − 1). For the non-uniform
nuclear matter calculation (A → ∞), we use the bare
nucleon’s mass. The fast Fourier transform is utilized for
calculation of the Coulomb potential, which is well suitable for periodic systems. For the calculation of the isolated finite nucleus, we use the method same as Ref. [33]
following the idea given in Ref. [34].
In order to make the Chebyshev polynomial expansion,
we need to set the maximum and the minimum singleparticle energies. The minimum energy is taken as emin =
−50 MeV, which is safe enough in the cases of N = Z
nuclei and the symmetric nuclear matter. The maximum
energy is set as
~2  π  2
emax = 3 ×
,
(25)
2m h
according to the maximum kinetic energy for the 3D grid
of h3 .
For the self-consistent iteration of the finitetemperature HF calculation, we use the modified Broyden’s method [35, 36]. We use the HF potential v(r) as
the Broyden’s vector to update [36].
B.

Validity check for the polynomial expansion

Let us first examine the accuracy of the expansion with
the Chebyshev polynomials. According to Eq. (22), the
maximum degrees of the polynomials M are adopted as
M = 1.5 × er β, which corresponds to the accuracy of
10−5.5 . In Fig. 1, we show the approximated Fermi-Dirac
distribution function fβµ () with the chemical potential
µ = −10 MeV (panel (a)), and the deviation from the
exact values (panel (b)),


M
a0 X
 − ec
df () =
+
an Tn
− fβµ ().
(26)
2
er
n=1
We assume that the maximum and minimum singleparticle energies are max = 500 MeV and min = −100
MeV, respectively, which leads to M = 900β (β in units
of MeV−1 ). The largest deviation appears around  = µ
and its value is order of 10−6 , which is consistent with
the estimation of Eq. (22). We can clearly see the importance of the temperature-dependent maximum degrees
M.

6
1
0.8

T=10 MeV
T= 1 MeV
T=0.1 MeV

0.6

0.4

0.2

0.2

-40

-20

0

20

0

40

(b) Error

2x10-6

-20

0

20

40

-20

0
ε [ MeV ]

20

40

(b) Error
5x10-6
dg(ε)

1x10-6
0
-1x10-6

0
-5x10-6

-2x10-6
-3x10-6

-40

1x10-5

3x10-6

df(ε)

0.6

0.4

0

T=10 MeV
T= 1 MeV
T=0.1 MeV

(a) Entropy function
0.8

g(ε)

fβµ(ε)

1

(a) Fermi-Dirac function

-40

-20

0
ε [ MeV ]

20

-1x10-5

40

FIG. 1.
(a) Fermi-Dirac distribution function fβµ () as
a function of the single-particle energy, calculated with the
Chebyshev polynomial expansion. The degrees of polynomials are M = 90, 900, 9,000 for T = 10, 1, and 0.1 MeV,
respectively. (b) Error in the polynomial approximation, Eq.
(26).

-40

FIG. 2. Same as Fig. 1, but for the entropy function g() of
Eq. (27) instead of fβµ ().

2x10-6

M=90
M=100
M=120

1.5x10-6
1x10-6

We perform the same analysis on the function,
df(ε)

g() = −fβµ () ln fβµ () − {1 − fβµ ()} ln {1 − fβµ ()} ,
(27)
which is used for calculation of the entropy, Eq. (17),
and show the result in Fig. 2. The maximum degrees
M 0 is taken as M 0 = M = 900β. The Chebyshev expansion for g() is well approximated with the deviation
is smaller than 10−5 for any temperature. They are well
controlled as far as the maximum degrees M are adjusted
in proportion to β.
If we increase M at the fixed temperature, the accuracy is significantly improved as shown in Fig. 3. We
find an agreement with Eq. (22); M = 100 (M = 120)
corresponds to D = 6 (D = 7) in Eq. (22).
In order to check the accuracy in the final results, we
perform the FOE calculation using M = k × er β with
different values of k. We use the Woods-Saxon potential,
Vws (r) = Vws /(1 + exp(r − Rws )/aws ), with Vws = −50
MeV, Rws = 3 fm, and aws = 0.5 fm. The adopted space
size is (13 fm)3 and the chemical potential is fixed at µ =
−15 MeV. The calculated quantities with T = 1 MeV and
10 MeV are shown in Table I. The nucleon number is
calculated as the integration of density over the adopted
space.R The “Woods-Saxon energy” is defined as Ews =
(1/2) Vws (r)ρ(r)dr. The difference between k = 1.5
and k = 2.0 is negligible, less than 100 eV in energy. In
this paper, we use the temperature-dependent maximum

5x10-7
0
-5x10-7
-1x10-6
-1.5x10-6
-2x10-6
-100

-50

0
ε [ MeV ]

50

100

FIG. 3. Error in the Fermi-Dirac function at T = 10 MeV
with different values of the maximum degrees M . The magenta line is calculated with M = 90, the same as that in
Fig. 1(b). Those with M = 100 and M = 120 are shown by
orange and light blue lines, respectively.

degrees, M = 1.5 × er β, which provides a reasonable
accuracy.
With the same chemical potential, the nucleon number A increases approximately threefold from T = 1 MeV
to 10 MeV. On the other hand, difference in the WoodsSaxon energy is only about 30%. This can be understood
from the nucleon density profile shown in Fig. 4. A significant portion of nucleons is dripped from the WoodsSaxon potential at T = 10 MeV. The particles far out of
the potential range Rws = 3 fm do not contribute to Ews .

7
0.18

[35]. In Ref. [36], its performance has been examined for
finite nuclei with the Skyrme-HFB calculations at zero
temperature using the matrix diagonalization. We perform similar study for the FOE calculation at finite temperature. In Fig. 5, we show the convergence behaviors of
the modified Broyden’s method, compared to the linear
mixing method. Here, we show the difference in diagonal
density ρ(r) between the current (n) and the previous
(n − 1) iteration steps.
Z




(28)
|∆Fn | ≡ A−1 dr ρ(n) (r) − ρ(n−1) (r) .

T=1 MeV
T=10 MeV

0.16
0.14

ρ(r) [ fm-3 ]

0.12
0.1
0.08
0.06
0.04
0.02
0

0

1

2

3

4
r [ fm ]

5

6

7

8

FIG. 4. Calculated density profiles at T = 1 MeV and 10
MeV as a function of the radial coordinate. Symbols indicate calculated values at the square mesh points and lines are
obtained by the spline interpolation. See text for details.

C.

Isolated nuclei at finite temperature
1.

Modified Broyden’s method

In the HF calculations, the self-consistency between
the densities and the potentials (Hamiltonian) is reˆ (n) at the nquired. For a given HF Hamiltonian H
in
th iteration, the one-body density ρij is obtained using the FOE calculation of Eq. (16), which produces
(n)
ˆ out
. The self-consistency is
a new HF Hamiltonian H
achieved when we reach the fixed point for the Hamiltonian, Hin = Hout , which is equivalent to the density
fixed point. Since naive iteration with the total replace(n)
ˆ (n+1) = H
ˆ out
ment of the Hamiltonian as H
does not conin
verges in many cases, the linear mixing is often adopted
(n)
ˆ (n+1) = (1−α)H
ˆ (n) +αH
ˆ out
with a mixing parameter
as H
in
in
α. Although the divergence can be avoided if we choose
the parameter α small (0 < α  1), the convergence can
be very slow.
In this paper, we use the modified Broyden’s method

Since the chemical potential is adjusted every iteration
to reproduce either the average nucleon density ρav , or
the nucleon number, the baryon (nucleon) number A is
fixed during the iteration.
In the linear mixing, the result depends on the magnitude of the mixing parameter α. Figure 5 shows the case
of 16 O at T = 1 MeV. In this case, the calculation with
α = 1 does not converge, while that with α = 0.5 gives
the fastest convergence among α = 1, 0.8, 0.5, and 0.2.
The optimum value of α varies and is difficult to predict.
For instance, in the case of T = 10 MeV, the calculation
with α = 1 converges faster than that with α = 0.5. In
order to guarantee the convergence, we need to choose a
small value of α, typically α < 0.2, which leads to a slow
convergence of the iterative procedure.
The modified Broyden’s method provides a faster and
a stable convergence. The modified Broyden’s algorithm
[35] also contains two parameters we need to choose,
namely, the mixing parameter α and the maximum number of stored vectors m. It turns out that the result does
not strongly depend on the choice of α and m. For the
mixing parameter α, we can safely choose α ≈ 1. Larger
values of m give slightly better convergence, but only a
few iteration number difference between m = 10 and 100.
In the present paper, we adopt α = 0.8 and m = 100.
Although storing m Broyden’s vectors may require large
memory resources when the system size is large, the computational time for the Broyden’s procedure is negligible.

2.
TABLE I. Calculated values of nucleon number A, the total
kinetic energy Ekin , the Woods-Saxon energy Ews , and the
entropy S/kB , using different values of k. The temperature is
T = 1 MeV for the upper three rows, while T = 10 MeV for
the rest. See text for details.
k

A

1.0
1.5
2.0

9.8125
9.8125
9.8125

1.0
1.5
2.0

27.3152
27.3152
27.3152

Ekin [ MeV ]
Ews [ MeV ]
T = 1 MeV
159.5155
−183.8843
159.5157
−183.8843
159.5157
−183.8843
T = 10 MeV
532.9096
−248.7950
532.9055
−248.7954
532.9055
−248.7954

S/kB
8.3121
8.3121
8.3121
78.6638
78.6634
78.6634

Isolated doubly magic nuclei at finite temperature

First, let us show results of 16 O at finite temperature.
We use the space size of (13 fm)3 with the 3D cubic grid
of (1 fm)3 . In Fig. 6, the total energy E and the free energy F at every iteration are plotted. At T = 1 MeV, as
the iteration number increases, both E and F decrease
to the final values, E = −132.6 MeV and F = −132.9
MeV. The calculated entropy is very small, about 0.3kB .
This is due to the doubly closed-shell nature of 16 O. At
T = 10 MeV, in contrast, the total energy E increases
to reach the converged value, E = 210.0 MeV. Note that
the line in Fig. 6 is shifted downwards by 330 MeV to be
presented in the same panel as T = 1 MeV. Nevertheless,
the free energy F decreases, because the entropy gradu-

8
1

0.16
α=1
α=0.8
α=0.5
α=0.2
Broyden

0.01

0.12

ρ(r) [ fm-3 ]

0.0001
|ΔFn|

T=0.5 MeV
T=1.0 MeV
T=4.0 MeV
T=7.0 MeV
T=7.4 MeV
T=7.46 MeV
T=7.47 MeV
T=8.0 MeV

0.14

1x10-6

0.1
0.01

0.08
0.06

1x10-8

0.005

0.04
0

1x10-10

1x10-12

0.02

4

5

4
r [ fm ]

5

6

7

8

0
0

20

40

60

80
100
120
Itenation number

140

160

180

200

FIG. 5. Comparison of convergence between the linear mixing (black lines) and the modified Broyden’s method (red
symbols) for 16 O calculated at T = 1 MeV. The space size
is (13 fm)3 . See text for details.
-120
-122

E (T=1 MeV)
F (T=1 MeV)
E (T=10 MeV)
F (T=10 MeV)

-124

0

1

2

3

6

7

8

FIG. 7.
Nucleon density distribution of 16 O at different
temperature. The horizontal axis is the distance from the
center of mass. Since those calculated at T = 0.5 MeV and 8
MeV are indistinguishable from those at T = 1 MeV and 7.47
MeV, respectively, they are shown as symbols. The space size
is (13 fm)3 .
Inset: Density distributions in the outer region of r > 4 fm
is shown, but those at T = 0.5 MeV and 8 MeV are omitted.
See text for details.

E or F [ MeV ]

-126
-128
-130
-132
-134
-136
-138
0

2

4

6

8
10
12
Iteration number

14

16

18

20

FIG. 6.
Total energy E (solid lines) and free energy F
(dashed) of 16 O at T = 1 MeV (blue) and T = 10 MeV
(red) as functions of self-consistent iteration number with the
modified Broyden’s method. The lines for T = 10 MeV are
shifted by −330 MeV for E and by +300 MeV for F . See text
for details.

ally increases as the iteration proceeds. The entropy is
calculated as S = 64.6kB .
The calculated nucleon density distributions are presented in Fig. 7. The center of mass of 16 O is located
at the center of the cubic box of (13 fm)3 . The density
values at the calculated grid points are shown by circles
for T = 0.5 and 8 MeV. The spline interpolation is used
to show the smooth lines in Fig. 7. At low temperature,
such as T = 0.5 and 1 MeV, we find a signature of the
shell effect as a dip at the center of the nucleus. This is
due to the full occupation of the 0p orbitals. Higher is
the temperature, more fractional is the occupation, leading to weakening of the shell effect. At T = 7 MeV, the
density hole at the center disappears.
The phase transition to the uniformed nuclear matter
takes place at the critical temperature Tc = 7.47 MeV.
More precisely speaking, it is 7.46 < Tc ≤ 7.47 MeV.

A discontinuous change in the density profile suddenly
occurs at T = Tc . This is a consequence of the selfconsistent evolution of the mean-field potential, which
gives a striking contrast to the density change in the fixed
potential (Fig. 4). The 16 O nucleus is in a liquid phase at
T = 0. Since there are some dripped nucleons at T 6= 0, it
is a coexistence phase of liquid and vapor at 0 < T < Tc .
Then, it is transformed to the gas phase at T > Tc . We
should note here that the critical temperature Tc depends
on the adopted volume V that is (13 fm)3 in the present
calculation. Tc for the isolated nucleus should be given
as the value at V → ∞. See Sec. III D 2 for more details.
The dripped nucleons at different temperature can be
seen in the inset of Fig. 7. In the gas phase, the density
should be ρgas = 16/(13 fm)3 = 7.28 × 10−3 fm−3 . The
uniform density obtained at T > Tc is very close to ρgas ,
however, the calculated density is not perfectly constant.
It has a minimum value at the center and slightly increases as r increases. This strange behavior is an artifact
due to the finiteness of the box size. Following the idea
of Ref. [34], the Coulomb potential for the isolated system is calculated by assuming that there exist no charge
outside of the adopted space ((13 fm)3 in the present
case). Therefore, the charged particles (protons) tend to
move toward the edge of the box, in order to reduce the
Coulomb repulsive energy. We have also confirmed that
the density is perfectly constant if the Coulomb potential
is neglected.
Figure 8 shows the density profiles for 40 Ca. The
model space is taken as (17 fm)3 with the cubic grids
of (1 fm)3 . The shell effect opposite to 16 O is seen
at low temperature, namely a bump at the center of
the nucleus. This is due to the full occupation of 1s
orbital. The shell effect becomes invisible at T = 7

9
50

0.18

ρ(r) [ fm-3 ]

0.14
0.12

MeV
MeV
MeV
MeV
MeV

0.1
0.08
0.06
0.04
0.02
0

0

1

2

3

4

5

6

7

axial

30

20
triaxial

8

0

FIG. 8.
Nucleon density distribution of O at different
temperature. The horizontal axis is the distance from the
center of mass. The space size is (17 fm)3 .
800
600
400
200

E
Ekin

0

-ST+600MeV
F

-200
-400

-344.35

-600

-344.4

-800

-344.45
0

0

0.5
2

4

6

8

10

T [ MeV ]

FIG. 9. Total energy E (red solid) and free energy F (blue
solid) are shown as functions of T , together with kinetic energy Ekin (black dotted) and −ST (purple dash-dotted) for
40
Ca. Note that the line of −ST is shifted upward by 600
MeV.
Inset: E and F in the zero temperature region are magnified.

MeV. The critical temperature of the liquid-gas transition is located in 8.5 < Tc ≤ 8.6 MeV. The discontinuous density change is seen at T = Tc . The density of
the uniform matter in the present calculation should be
ρgas = 40/(17 fm)3 = 8.14 × 10−3 fm−3 .
In Fig. 9, the energy E and the free energy F are shown
as functions of T . At T = Tc ≈ 8.6 MeV, the total energy
E shows a kink because of the abrupt density change.
However, this kink is almost canceled by an opposite kink
behavior in the entropy, and F = E − ST behaves rather
smoothly. E is a monotonic increasing function of T ,
while F is a decreasing function. The zero temperature
limit is easily achieved in this case, because the 40 Ca
nucleus is doubly magic with large shell gaps at the Fermi
surface. In the inset panel of Fig. 9, we find that E =
F holds in very high accuracy at T = 0.1 MeV. The
difference is within 0.1 eV. Even at T = 0.5 MeV, it is
within 50 keV. Note that the upper (lower) limit of the

axial
0

0.5

1

1.5

2

2.5

3

T [ MeV ]

16

E [ MeV ]

|Q2|
|Q22|

10

r [ fm ]

-1000

triaxial

40
Q momdent [ fm2 ]

T=1.0
T=3.0
T=7.0
T=8.5
T=8.6

0.16

FIG. 10. Calculated intrinsic quadrupole moments (Q2 and
|Q22 |) for 24 Mg at finite temperature. At low temperature
(T . 0.5 MeV), the self-consistent iterations starting from
different initial states result in different solutions, labeled by
“axial” and “triaxial”, respectively. See text for details.

total energy at T = 0 is given by E (F ) at T > 0.
It is worth mentioning that the coordinate-space representation is essential to describe the dripped nucleons and the liquid-gas phase transition. Most of the
finite-temperature mean-field calculations in the past
have been performed with the harmonic oscillator basis
[3, 4, 6]. Those studies are focused on the shape change
and the pairing properties at finite temperature, however, it is difficult to describe the uniform matter and
the dripped nucleons. In Refs. [8, 9], adopting the spherical Wigner-Seitz approximation, the finite-temperature
Skyrme Hartree-Fock calculation was performed in the
radial coordinate representation. Our results on properties of the liquid-gas phase transition turns out to be substantially different from Refs. [8, 9]. For instance, they
showed that ignoring the Coulomb potential for 208 Pb
leads to a significant increase in Tc (about 5 MeV) and
a smooth continuous transition from the liquid to the
gas phase. In our calculation, properties of the liquidgas transition is almost invariant, even if we neglect the
Coulomb potential: We find a slight increase by only a
few hundreds of keV with a discontinuous transitions of
the density profile into the uniform matter. Although
the the results in this paper are on N = Z nuclei only, it
would be important to perform the detailed comparison
in future to identify origins of the discrepancies.

3.

Isolated deformed nuclei at finite temperature

We calculate an isolated 24 Mg nucleus at finite temperature, which is known to be deformed in the ground
state. In Fig. 10, we present the calculated quadrupole

10

400
-186

E [ MeV ]

200

axial

-188
-190

0

-192

triaxial
0

0.2

0.4

0.6

0.8

E
Ekin
F

1

-200

-400
0

1

2

3

4

5

6

7

8

T [ MeV ]

FIG. 11. Total energy E (red solid line) and free energy F
(blue solid) are shown as functions of T , together with kinetic
energy Ekin (black dotted) for 24 Mg.
Inset: E and F near the zero temperature are magnified.

moment, which is defined as
v
u 2
Z
uX
Q2 ≡ t
|Q2µ |2 , Q2µ ≡ drr2 Y2µ ρ(r).

(29)

µ=−2

At low temperature, when we start the self-consistent
iteration with a Hamiltonian corresponding to an axially
symmetric deformed density distribution, the calculation
converges to an axially symmetric prolate nucleus (Q2µ =
0 except for µ = 0). However, near the zero temperature,
this does not correspond to the state with the minimum
free energy. If we start with a triaxial shape, it ends
up with a triaxially deformed nucleus, characterized by
Q22 6= 0. The shape transition from triaxial to axial
shapes takes place at temperature T = Ttri with 0.5 <
Ttri < 0.6 MeV. The axial prolate shape persists till the
second shape transition to the spherical shapes, which
takes place at temperature T = Tdef with 2.7 < Tdef <
2.8 MeV. This is shown in Fig. 10 as two lines.
In Fig.11, we show the temperature dependence of the
energy E and the free energy F . A kink of the energy E is
caused by the liquid-gas phase transition. The calculated
critical temperature is 6.5 < Tc < 6.6 MeV. Again, this
kink is almost canceled by an opposite kink behavior in
the entropy, and a kink in the free energy F is much
smaller. The effect of the shape transition at T ≈ 2.7
MeV is invisible in the temperature dependence of E and
F , while that of the axial-triaxial transition at T ≈ 0.5
MeV can be seen in the inset of Fig.11. An extrapolation
to T = 0 using calculations at T > 0.5 MeV may lead to
a wrong answer. The zero temperature limit should be
carefully examined when a structure change is expected
at very low temperature.
Another interesting feature is a property of the kinetic
energy. At T = 0.1 MeV, the kinetic energy for the triaxial solution is smaller than that of the axial one by
about 6.3 MeV, while the difference in the total energy
is about 2.4 MeV. This clearly indicates that the triax-

ial shape in 24 Mg is realized by significant decrease in
the kinetic energy, although it is unfavored by the potential energy. As the deformation decreases with increasing
temperature, the kinetic energy monotonically increases
up to T = Tdef . This also suggests that the deformation reduces the kinetic energy, while the potential energy favors the sphericity. This is consistent with the fact
that the Thomas-Fermi method cannot produce the deformation. When there are more nucleons moving along
z direction than x and y directions, according to the uncertainty principle, the kinetic energy can be reduced by
elongating a potential in the z direction. This effect is
completely lost in the local density approximation.
At T > Tdef , the nucleus is spherical and the kinetic
energy decreases as increasing T . At T = Tc , since the
nucleus suddenly breaks up into a gas phase, the kinetic
energy shows a discontinuous drop. This is because number of dripped nucleons increases as a function of T up
to T = Tc , and their momenta are smaller than those of
nucleons confined inside the nucleus, which is also due to
the uncertain principle.

D.

Non-uniform periodic nuclear matter at finite
temperature

Next, we apply the method to the non-uniform symmetric nuclear matter. The only difference from the calculations in Sec. III C is the treatment of the Coulomb
potential. For periodic non-uniform nuclear matter, we
assume the uniform distribution of electrons, to guarantee the charge neutrality. This results in the vanishing
k = 0 Fourier component of the Coulomb potential. In
the present calculations, the electron energy does not affect the structure of nuclear matter, since we calculate
the nuclear matter at given baryon density ρ and proton
ratio (Yp = 0.5).

1.

A = 32 in a cell of (17 fm)3

First, we calculate the symmetric nuclear matter at
average baryon density ρ = 6.51 × 10−3 fm−3 with a simple cubic initial configuration in which a 32 S nucleus is
located at the center of a cubic box of (17 fm)3 . At low
temperature, we find the 32 S nucleus in a self-consistent
solution. The 32 S nucleus is deformed at low temperature
T < Tdef . The deformation disappears at T = Tdef with
1.6 < Tdef < 1.7 MeV. The dripped nucleons increase
with temperature, then, the liquid-gas phase transition
takes place at T = Tc with 7.7 < Tc < 7.8 MeV. See
Fig. 12 for evolution of the density distribution as a function of temperature.
In addition to the simple cubic configuration, we also
perform calculations with the body-centered-cubic (bcc)
configuration as the initial state. This leads to two 16 O
nuclei in the same cell (17 fm)3 . At low temperature,
the self-consistent calculation converges to the bcc phase.

11
T=2.0
T=5.0
T=7.0
T=7.7
T=7.8

0.16
0.14

MeV
MeV
MeV
MeV
MeV

y [ fm ]

0.18

ρ(r) [ fm-3 ]

0.12
0.1
0.08
0.06
0.04

0

0

1

2

3

4

5

6

7

8

r [ fm ]

FIG. 12. Density distributions at ρ = 6.51 × 10−3 fm−3 in
a cubic configuration at finite temperature. The cell size is
(17 fm)3 , and the horizontal axis represents the distance from
the center of the cell.

Since the (free) energy is larger than that of the cubic
configuration, the bcc state exists as a metastable equilibrium. Panels (a) and (b) in Fig. 13 show the density
distributions at T = 0.1 MeV in the xy plane at z = 8 fm
and at z = 0. In contrast, panels (c) and (d) in Fig. 13
show those at T = 5 MeV, indicating that the bcc state
is no longer stable at higher temperature. During the
self-consistent iterations starting from the bcc state, the
16
O nucleus at the center disappears leading to the cubic
configuration, namely, a single 32 S nucleus in the cell of
(17 fm)3 . The stability of the bcc state seems to be lost
around T = 4 MeV.
Another calculation with the initial configuration of a
40
Ca nucleus located at the center of the cell of (17 fm)3
is performed. This corresponds to the average density
of ρ = 8.14 × 10−3 fm−3 . The variation of the density
distribution as a function of temperature is similar to the
one for the isolated 40 Ca nucleus in Fig. 8. However, the
critical temperature for the liquid-gas transition slightly
increases, 8.6 < Tc < 8.7 MeV.
In the present calculation, the number of particles is
irrelevant to the computational cost. Thus, as far as the
cell and the grid sizes are invariant, the computing time
is roughly the same for any density and particle numbers
in the cell. It should be noted that there is no spurious effect in dripped nucleons, namely a rise up of the
density near the boundary observed in cases of isolated
nuclei (Sec. III C 2) , such as Fig. 7. The Coulomb potential in the periodic systems is influenced by the periodic
presence of other nuclei outside of the cell (17 fm)3 . The
dripped nucleons produce perfectly flat density distribution outside of the nucleus.

2.

A = 32 in a cell of (23 fm)3

Enlarging the cell size into (23 fm)3 keeping the baryon
number A = 32 in the cell, we perform the same calcula-

y [ fm ]

0.02

�
�
�
�
�
-�
-�
-�

�
�
�
�
�
-�
-�
-�
-�

(a)

(b)
�.��
�.�

(c)

-� -� -� -� � � � � �

x [ fm ]

(d)

�

-� -� -� � � � � �

x [ fm ]

FIG. 13. Density distributions in the xy plane at ρ = 6.51 ×
10−3 fm−3 with a cell of (17 fm)3 , calculated with the bcc
initial state. (a) T = 0.1 MeV and z = 8 fm, (b) T = 0.1
MeV and z = 0, (c) T = 5 MeV and z = 8 fm, (d) T = 5
MeV and z = 0.

tions with the cubic and bcc initial configurations. The
average density is ρ = 2.63 × 10−3 fm−3 . Both the cubic and bcc configurations exist at low temperature at
T . 3.7 MeV. The solution with a single 32 S nucleus
in the cell has lower energy than the bcc solution. In
the calculation with T ≥ 3.8 MeV, the bcc metastable
solution seems to disappear, since we end up with the
single 32 S nucleus in the cell even if we start with the bcc
configuration with two 16 O nuclei. The density profiles
obtained with calculations starting from the bcc initial
configuration are shown in Fig. 14 at T = 2 MeV (panels
(a) and (b)) and at T = 5 MeV ((c) and (d)).
Figure 15 presents the free energy per nucleon F/A for
various phases. The cubic configuration of 32 S has the
lowest free energy at T < Tc ≈ 5.1 MeV. At T > Tc ,
the uniform symmetric matter becomes the lowest. This
critical temperature Tc is significantly smaller than Tc
for the average density ρ = 6.51 × 10−3 fm−3 with the
cell size (17 fm)3 . This can be understood as follows:
For the uniform phase at T > 5 MeV, the system is well
approximated by the classical gas. The entropy of the
classical ideal gas has the volume dependence as S ∼
AkB ln(V /A). Thus, the entropy per nucleon S/A for
the cell of (23 fm)3 is larger than that of (17 fm)3 , by
δ(S/A) = kB ln(233 /173 ) ≈ 0.9kB . This leads to a shift
of 0.9kB T in the free energy of the uniform matter in
(17 fm)3 , shown by a dashed line in Fig. 15. Since the
entropy in the localized phases, such as cubic and bcc, is
scarcely affected by the volume change, Tc , given by the
crossing point of the uniform and cubic phases, decreases
as the volume increases.

12

y [ fm ]

11

(b)

(c)

(d)

0

-11
11

y [ fm ]

(a)

0.1

0.01

1e-3

0

1e-4
1e-5

-11
-11

0

11 -11

x [ fm ]

11

0

x [ fm ]

FIG. 14. Density distributions in the xy plane at ρ = 2.63 ×
10−3 fm−3 with a cell of (23 fm)3 , calculated with the bcc
initial state. (a) T = 2 MeV and z = 11 fm, (b) T = 2 MeV
and z = 0, (c) T = 5 MeV and z = 11 fm, (d) T = 5 MeV
and z = 0. Note that the color map is given in logarithmic
scale.

F/A [ MeV ]

0
-5

bcc

-10

F/A

(un

cubic

2

m)

+0

.9 T

un

-15

-20

ifor

ifo

3

4

5

T [ MeV ]

6

rm
7

FIG. 15.
Free energy per particle of symmetric nuclear
matter at ρ = 2.63 × 10−3 fm−3 with a cell of (23 fm)3 , for
cubic, bcc, and uniform phases. The dashed line is given by
shifting the line of the uniform matter by 0.9 ∗ T . See text for
details.

E.

Nearsightedness

Finally, let us check properties of the ”nearsightedness” in calculations of nuclear matter at finite temperature, then, examine whether it benefits calculations
of neutron star matter. The O(N ) calculation can be
achieved if the one-body density matrix, ρ(r, r0 ), is localized in a space considerably smaller than the cell size. In
the present calculation, we can truncate the Hamiltonian
matrix in Eq. (15) into a space only nearby |jn i. See also

arguments in Sec. II E.
Figure 16 presents the off-diagonal behaviors of the
density matrix for a 40 Ca nucleus located at the center of
the cell (23 fm)3 of the simple cubic lattice. For comparison, those for the uniform matter is shown in the bottom
panels (b) and (d). We adopt the center of the cell r = 0
as a reference point and show ρ(r, 0) as a function of the
distance r. The magnitude of the off-diagonal density
exponentially decays. For the uniform matter, the calculated behaviors indicate the decay constant proportional
to the temperature T . This is known in studies of finitetemperature density matrix for electrons in metals [37].
In contrast, for non-uniform matter with 40 Ca in a cell,
the decay is significantly faster than the uniform matter
and is insensitive to the temperature. This is not entirely
attributed to the finite radius of the nucleus 40 Ca. From
T = 1 MeV to T = 5 MeV, the radius of 40 Ca is reduced
by about 0.5 − 1 fm (cf. Fig. 8), because more nucleons are dripped to form low-density matter. However,
the effect of this reduction in the nuclear radius is not
visible in Fig. 16 (a) and (c). The fast decay may be,
at least partially, due to large density inside the nucleus.
At zero temperature, the uniform matter is expected to
show an oscillating pattern of ρ(r, 0) ∼ kF cos(kF r)/r2
[28], where kF is the Fermi momentum. Thus, at larger
density (larger kF ), the off-diagonal density goes to zero
more quickly.
Eventually, the localization of the density matrix is
more prominent in the non-uniform phase than in the
uniform matter. Adopting the cut-off value for the relative magnitude as 10−4 (the dashed line in Fig. 16 (d)),
the cut-off distance for the uniform matter is given by
Rc ≈ 13 fm at T = 5 MeV, and it is considerably larger
than 20 fm at T = 1 MeV. In contrast, for the nonuniform matter, the cut-off distance is Rc ≈ 10 fm at
T = 5 MeV and Rc ≈ 13 fm at T = 1 MeV. When
we calculate the ρ(r, r0 ) with the recursion relation (15),
truncating the space into a local subspace, |r − r0 | < Rc ,
may lead to a sizable reduction in the computational cost
if the cell size is larger than Rc3 .

IV.

CONCLUSION

We examine the applicability and the usefulness of the
Fermi operator expansion (FOE) method in nuclear energy density functional approaches at finite temperature.
The one-body density matrix, which is identical to the
Fermi operator, is expanded in terms of the Chebyshev
polynomials up to the finite order. The maximum degree
of the polynomials is inversely proportional to the temperature. Thus, it becomes extremely efficient for calculations at high temperature. For the self-consistent iteration procedure, we adopt the modified Broyden’s mixing
method. The same idea of the polynomial expansion is
applied to calculations of the entropy, which enables us
to estimate the free energy without diagonalization of
the Hamiltonian matrix. The FOE method is applied to

13
1

1

0.6
0.2

|ρ( r,0)| / ρ (0,0)

T=1 MeV
T=3 MeV
T=5 MeV

0.4

ρ(r,0)/ρ(0,0)

(c) 40Ca

(a) 40Ca

0.8

0

(b) uniform

0.8
0.6
T=

0.4

T=1 MeV
T=3 MeV
T=5 MeV

10 -8

10 -2
T=
1

10 -6

10

0

M

0.2

10 -4

eV

0
0

2

4

6

r [ fm ]

8

10

10 -10

M
eV

(d) uniform
0

5

10

r [ fm ]

15

20

FIG. 16. Normalized density matrix ρ(r, 0)/ρ(0, 0) as a function of r for the simple cubic configuration ((a) and (c)) and
for the uniform matter ((b) and (d)). Panels (c) and (d)
show the absolute values in logarithmic scale, with the dashed
lines indicating the value of 10−4 . The average density is
ρ = 8.22 × 10−3 fm−3 which corresponds to A = 40 (40 Ca) in
a cell of (23 fm)3 .

calculations of isolated nuclei and non-uniform nuclear
matter, using the 3D coordinate-space representation.
We investigate thermal properties of isolated nuclei in
a cell of (13 fm)3 . For 24 Mg, the triaxial shape has the
minimum energy at zero temperature. The triaxial state
exists as a solution at T . 0.6 MeV, beyond that, the
state disappears. The axial deformed solution survives
till T ≈ 2.7 MeV, beyond which the nuclear shape is
spherical. The liquid-gas transition takes place around
Tc ≈ 6.5 MeV. For doubly-magic spherical nuclei, such
as 16 O and 40 Ca, the critical temperature of the liquidgas transition has slightly higher values, Tc = 7 − 9 MeV.
However, the detailed values of the critical temperature
may not have a significant meaning for isolated nuclei, because they depends on the volume of the adopted space.
We need to take the infinite volume limit. Nevertheless,
it is a great advantage of the coordinate-space representation to be capable of describing both the spatially localized nucleus and the extended matter.
For periodic non-uniform nuclear matter, the calculations are performed with different cell sizes, (13 fm)3
and (23 fm)3 , with the same nucleon number A = 32.
We start the self-consistent iteration with different initial
states, such as the simple cubic and the bcc configurations. At low temperature, both the simple cubic and
the bcc states exist as self-consistent solutions. The cubic state is lower in free energy than the bcc state. The
transition to the uniform matter takes place at Tc , the
value of which is smaller for a larger cell. This volume
effect on the critical temperature Tc is due to the volume
dependence of the entropy of the uniform matter. For the
inner crust of neutron stars in the beta equilibrium, the
cell size is supposed to decrease as the density increases
[19]. Since the entropy of the classical gas behaves as
S/A ∼ kB ln(V /A) ∼ −kB ln ρ, Tc may become larger at
larger densities. This is somewhat opposite to our naive

expectation, because the density profile becomes flatter
at higher density. It may be of interest to investigate the
critical temperature Tc at different density regions.
Advantageous features of the FOE method in computational point of view can be summarized as follows: (1)
The matrix diagonalization is not involved in the calculation, including the calculation of the entropy. (2)
The calculation of the density matrix ρij is independent
with respect to the index j. Thus, it is suitable for the
distributed-memory parallel computing. (3) The computational cost could scale linearly with respect to the space
dimension N , when N is large enough. Here, N is the
dimension of the matrix ρij .
The last point (3) above is numerically investigated
by examining the decay of the density matrix ρ(r, r0 ) at
large |r − r0 |. For the uniform matter, the decay length
is shorter at higher temperature, which has been known
for electron systems [28]. In addition, for non-uniform
matter with localized nuclei, the decay is significantly
faster than the uniform matter. The decay pattern of the
non-uniform matter at T = 1 MeV is close to that of the
uniform matter at T ≈ 10 MeV. The short decay length
of the density matrix could lead to the O(N ) calculation
by truncating the space in the matrix operation. The
O(N ) method may be more useful in the non-uniform
matter than in the uniform matter.
The calculations in the present paper use the BKN
energy density functional. It is straightforward to extend this to realistic Skyrme functionals, which is under
progress. The proper treatment of the matter in the periodic potential requires the band calculation. The density should be constructed by averaging the calculated
results over different values of Bloch wave numbers k,
which should be relatively easy to perform. The calculation can be further parallelized with respect to different
k.
The extension of the FOE method to the finitetemperature HFB calculation is formally straightforward
as well. This can be done by replacing the single-particle
Hamiltonian H by the HFB Hamiltonian,


H −µ
∆
Hµ =
,
(30)
−∆∗ −(h − µ)∗
in Eq. (7) to achieve the generalized density matrix R.
X
R=
|αifβ (Eα )hα| = fβ (Hµ ),
(31)
α≷0

where |αi are the quasiparticle eigenstates, Hµ |αi =
Eα |αi, and the summation is taken over both positive
and negative quasiparticle energies. However, there is
a practical issue to be examined in future, namely, the
truncation of the pairing model space. Since most of the
pairing energy functional has been constructed with a
cut-off energy, it is preferable to develop a prescription
to allow the truncation of the pairing model space.
The FOE method may open a new possibility for studies of the non-uniform baryonic matter at finite temperature and neutron-star matter in the crust region.

14
ACKNOWLEDGMENTS

This work is supported in part by JSPS KAKENHI
Grant No. 18H01209. This research in part used com-

[1] J.-P. Blaizot and G. Ripka, Quantum Theory of Finite
Systems (MIT Press, Cambridge, 1986).
[2] N. Schunck, ed., Energy Density Functional Methods for
Atomic Nuclei, 2053-2563 (IOP Publishing, 2019).
[3] A. L. Goodman, Nuclear Physics A 352, 30 (1981).
[4] J. L. Egido and P. Ring, Journal of Physics G: Nuclear
and Particle Physics 19, 1 (1993).
[5] G. Bertsch and J. Mehlhaff, Computer Physics Communications 207, 518 (2016).
[6] W. Zhang and Y. F. Niu, Phys. Rev. C 97, 054302 (2018).
[7] W. G. Newton and J. R. Stone, Phys. Rev. C 79, 055801
(2009).
[8] P. Bonche, S. Levit, and D. Vautherin, Nuclear Physics
A 427, 278 (1984).
[9] P. Bonche, S. Levit, and D. Vautherin, Nuclear Physics
A 436, 265 (1985).
[10] M. Brack, C. Guet, and H.-B. H˚
akansson, Physics Reports 123, 275 (1985).
[11] M. Onsi, H. Przysiezniak, and J. M. Pearson, Phys. Rev.
C 55, 3139 (1997).
[12] M. Okamoto, T. Maruyama, K. Yabana, and T. Tatsumi,
Phys. Rev. C 88, 025801 (2013).
[13] C.-J. Xia, T. Maruyama, N. Yasutake, and T. Tatsumi,
Phys. Rev. D 106, 063020 (2022).
[14] G. Watanabe, K. Sato, K. Yasuoka, and T. Ebisuzaki,
Phys. Rev. C 69, 055805 (2004).
[15] C. J. Horowitz, D. K. Berry, C. M. Briggs, M. E. Caplan,
A. Cumming, and A. S. Schneider, Phys. Rev. Lett. 114,
031102 (2015).
[16] M. E. Caplan and C. J. Horowitz, Rev. Mod. Phys. 89,
041002 (2017).
[17] A. Bulgac and P. Magierski, Nuclear Physics A 683, 695
(2001).
[18] N. Chamel, Nuclear Physics A 747, 109 (2005).
[19] J. Negele and D. Vautherin, Nuclear Physics A 207, 298
(1973).

putational resources provided through the HPCI System
Research Project (Project ID: hp200069), and by Multidisciplinary Cooperative Research Program in Center for
Computational Sciences, University of Tsukuba.

[20] M. Baldo, E. Saperstein, and S. Tolokonnikov, Nuclear
Physics A 775, 235 (2006).
[21] H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi, Nuclear Physics A 637, 435 (1998).
[22] H. Togashi, K. Nakazato, Y. Takehara, S. Yamamuro,
H. Suzuki, and M. Takano, Nuclear Physics A 961, 78
(2017).
[23] N. Chamel and P. Haensel, Living Rev. Relativity 11, 10
(2008).
[24] S. Jin, A. Bulgac, K. Roche, and G. Wlazlowski, Phys.
Rev. C 95, 044302 (2017).
[25] Y. Kashiwaba and T. Nakatsukasa, Phys. Rev. C 101,
045804 (2020).
[26] S. Goedecker and L. Colombo, Phys. Rev. Lett. 73, 122
(1994).
[27] S. Goedecker and M. Teter, Phys. Rev. B 51, 9455 (1995).
[28] S. Wu and C. Jayanthi, Physics Reports 358, 1 (2002).
[29] W. Kohn, International Journal of Quantum Chemistry
56, 229 (1995).
[30] R. Baer and M. Head-Gordon, Phys. Rev. Lett. 79, 3962
(1997).
[31] P. Bonche, S. Koonin, and J. W. Negele, Phys. Rev. C
13, 1226 (1976).
[32] T. Nakatsukasa and K. Yabana, Phys. Rev. C 71, 024301
(2005).
[33] J. Maruhn, P.-G. Reinhard, P. Stevenson, and A. Umar,
Computer Physics Communications 185, 2195 (2014).
[34] J. Eastwood and D. Brownrigg, Journal of Computational Physics 32, 24 (1979).
[35] D. D. Johnson, Phys. Rev. B 38, 12807 (1988).
[36] A. Baran, A. Bulgac, M. M. Forbes, G. Hagen,
W. Nazarewicz, N. Schunck, and M. V. Stoitsov, Phys.
Rev. C 78, 014318 (2008).
[37] S. Goedecker, Phys. Rev. B 58, 3501 (1998).