Self-consistent energy density functional approaches to the crust of neutron stars

In order to study structure of the crust in neutron stars, we develop a finite-temperature Skyrme-Hartree-Fock method in the full three-dimensional coordinate space using the Fermion operator expansion method. It provides us with a possible order-N approach to non-uniform neutron-star matters.


Introduction
The neutron star is a compact star whose radius is about 10 km, which can be regarded as a giant nucleus with macroscopic numbers of hadrons (baryons and mesons). The core part, the central region of the neutron star, is supposed to be uniform nuclear (hadron) matter. However, moving toward the surface, it changes into non-uniform nuclear matter composed of neutron-rich nuclei and neutron matter, called "inner crust". In the vicinity of the boundary between the core and the inner crust, exotic phases, called "nuclear pasta" are predicted. Since a new exotic structure may appear in certain conditions, the calculation without any symmetry restriction is desired.
Structure of these non-uniform exotic nuclear matter have been studied with the Thomas-Fermi (local density) approximation [1] and/or with the Wigner-Seitz approximation [2]. Although these two approximations are complimentary to each other, neither is able to take proper account of transport properties in the non-uniform nuclear matter. Thus, it is highly desirable to go beyond these approximations.
Recently, the three-dimensional (3D) coordinate-space solver of the Hartree-Fock-Bogoliubov (HFB) method has been developed [3]. In order to abstain from diagonalizing the large HFB Hamiltonian matrix, which is computationally demanding, the contour integration in the complex energy plane together with iterative solutions of the shifted linear algebraic equations is adopted to construct densities. The method has been extended to the finite-temperature HFB method, by including the Matsubara frequencies as the imaginary shifts [4]. The method was successfully applied to non-uniform nuclear matter with superfluid neutrons at finite temperature.
In this paper, we present another coordinate-space solver using the Fermion operator expansion (FOE) method [6]. The method was originally developed in the computational condensed matter physics, and known as one of order-N (O(N)) methods [5]. The conventional approaches to the density functional theory scales as N 3 , where N is the number of particles. This N 3 -scaling sets a severe limitation to the system size. In contrast, the O(N) method scales as N 1 in ideal cases. The foundation of the O(N) method is provided by "principle of nearsightedness" [5], namely a consequence of quantum mechanical decoherence among many particles. Calculation of the non-uniform neutron star matter may require us to treat more than thousands of nucleons (baryons) in a large space. In this paper, we demonstrate a test calculation of the FOE method for non-uniform nuclear matter at finite temperature.

Fermion operator expansion method
Let us briefly recapitulate the FOE method.Ĥ is the one-body mean-field Hamiltonian, satisfyingĤ |n = E n |n . The one-body density at a temperature T = 1/β is given bŷ Next, the Fermi-Dirac function is approximated by a series of polynomial functions {T j (x)}; f (x) = j a j T j (x). Then, the density, represented by f (Ĥ), can be written aŝ where T j (x) is a polynomial of the j-th degree. Here, the summation is truncated at the maximum degree M. In order to obtain a reasonable description of the Fermi-Dirac function, roughly speaking, M should increase at lower T . The density can be constructed by the M times operation of the Hamiltonian. We adopt the Chebyshev polynomials following ref. [6]. They can be determined by the following recursion relations: with T 0 (x) = 1 and T 1 (x) = x. Therefore, the operation of the polynomial of the Hamiltonian on a state |ψ , |ψ ( j) ≡ T j (Ĥ) |ψ , is calculated as with |ψ (0) = |ψ and |ψ (1) =Ĥ |ψ . The density is given bŷ We need to perform the calculation for, at least, N different states, {|ψ k , k = 1, · · · , N}. If the Hamiltonian is sparse, the total computational cost is estimated as O(N 2 ). However, when the "nearsightedness" exists, it can be reduced to O(N) [5].

Application of the FOE method
In this paper, we apply the FOE method to the BKN energy density functional [7] which is a simplified Skyrme-like functional without the spin-orbit coupling and assumes the identical densities of protons and neutrons. For simplicity, we neglect the pairing correlation, thus, it is the finite-temperature Hartree-Fock calculation, instead of finite-temperature HFB. The calculation is performed using the 3D coordinate-space representation discretized in the mesh [8]. The number of mesh points is 25 3 with the cubic mesh of (1 fm) 3 . The baryon density and the proton ratio are fixed at ρ B = 0.03 fm −3 and at Y p = 0.5 (symmetric nuclear matter), respectively.
We first set the temperature T = 5 MeV, then, perform the self-consistent iterative calculation starting from the body-centered configuration (bcc) as the initial state. The converged state also has a bcc-like structure: In the "sea" of symmetric nuclear matter whose density is about 8.3 × 10 −4 fm −3 , a "nucleus" (excess density) is embedded at the center and 1/8 of a nucleus is at each corner of the cubic box ( Fig. 1(b)). In total, there are two nuclei in the adopted space of volume V = 25 3 fm 3 . The number of nucleons in the space is ρ B × V = 468.75, so that, at every iteration, the chemical potential µ is adjusted to produce this particle number. Increasing the temperature ( Fig. 1(a)), the nucleons are gradually leaking from "nucleus" to the "sea", approaching to the uniform matter. It should be noted that, since the symmetric matter is assumed, the result is not realistic for neutron star inner crust. However, it is still meaningful for testing performance of the FOE method applied to nuclear matter.
The FOE method is suitable for parallel computing, because calculation of eq. (5) can be performed independently for different |ψ . The FOE method is extremely efficient at high temperature, since the Fermi-Dirac function becomes smooth and is well approximated by a small number of polynomials. The computational cost is roughly proportional to M in eq. (5). M is order of thousands at T ∼ 100 keV, while it is order ten at T ∼ 10 MeV. This work is supported by JSPS KAKENHI Grant No.18H01209 and No.19H05142. This research used computational resources provided by Joint Center for Advanced High Performance Computing (JCAHPC) through the HPCI System Research Project (Project ID: hp200069) and through Multidisciplinary Cooperative Research Program in Center for Computational Sciences, University of Tsukuba.