Three-dimensional mesh calculations for covariant density functional theory

In contrast to the non-relativistic approaches, three-dimensional (3D) mesh calculations for the {\it relativistic} density functional theory have not been realized because of the challenges of variational collapse and fermion doubling. We overcome these difficulties by developing a novel method based on the ideas of Wilson fermion as well as the variational principle for the inverse Hamiltonian. We demonstrate the applicability of this method by applying it to $^{16}$O, $^{24}$Mg, and $^{28}$Si nuclei, providing detailed explanation on the formalism and verification of numerical implementation.


I. INTRODUCTION
The self-consistent mean-field theory with phenomenological interactions has been successfully employed in nuclear physics in the past decades [1]. This theory is intimately related to the density functional theory (DFT) originally developed for many electron systems [2][3][4]. The DFT is based on the Hohenberg-Kohn theorem [2,4] for the existence of a universal energy density functional for many-body systems, which completely contains the many-body correlations in principle. In nuclear systems, a large part of the many-body correlations is taken into account through the parameters of phenomenological interactions, which are determined by fitting to a selected set of experimental data. The resultant Hartree-Fock equation can thus be regarded as the Kohn-Sham equation [3,4] in the DFT, which already contains the correlation effects even though it appears an equation for noninteracting systems. Notice that, starting from an energy density functional, one does not have to rely on a specific nucleon-nucleon interaction, which is an important feature in considering the density dependence of the energy functional. Energy functionals of Skyrme [5] and Gogny [6] types have been developed for the non-relativistic nuclear many-body calculations. The relativistic variant of DFT, referred to as the covariant density functional theory (CDFT), has also been developed [4,[7][8][9][10][11].
The nuclear DFT has several advantages over other theoretical methods. Firstly, since it reduces a problem of an interacting system to a problem of a non-interacting system, the numerical cost increases moderately with the number of particles in the system. Because of this, the DFT is the only method at present which is able to describe atomic nuclei in the whole nuclear chart at a reasonable computational cost. Secondly, formulated in the body-fixed frame, the DFT provides an intuitive view of nuclear deformation as a spontaneous breaking of symmetries. As a matter of fact, a variety of deformation de-grees of freedom are expected to play an important role in nuclear phenomena. For instance, it has been pointed out that the non-axial and reflection-asymmetric deformations play important roles in the fission barriers as well as the fission paths [12,13]. Also, the cluster states are expected to emerge in the excited states of sd-shell nuclei [14][15][16][17][18][19], and the exotic octupole and hexadecapole deformations are expected in several regions in the nuclear chart [20][21][22][23][24][25][26][27][28][29][30][31][32].
In order to describe various nuclear deformations flexibly and efficiently with DFT, a coordinate-space representation, in which the full real space is discretized into three-dimensional (3D) lattice, is most suitable. Notice that with the basis expansion method, in which wave functions are expanded on a finite basis set, the solutions strongly reflect the properties of the basis used, and the method is less powerful as compared to the coordinate-space representation. For the nuclear DFT calculations with the non-relativistic Skyrme functional in the 3D coordinate-space representation, the so-called imaginary time method has been successfully developed [33,34]. In this method, starting with arbitrary singleparticle wave functions, one obtains self-consistent solutions after evolving the wave functions along the imaginary time axis. It has been extensively applied to calculations not only for the ground state [35,36] but also for excited states of atomic nuclei. The latter includes exotic excitation modes in neutron-rich nuclei [37] and exotic structure in high-spin states [38][39][40]. The method has also been applied to deformation of Λ hypernuclei [41] as well as to deformed Hartree-Fock-Bogoliubov calculations with the two-basis method [26,42].
In contrast to the non-relativistic case with Skyrme functional, the CDFT calculations for deformed nuclei have been performed only with the basis expansion method, see, e.g., Refs. [12,13,31,[43][44][45][46][47]. That is, the CDFT calculations in the 3D coordinate space have not been realized yet. There are two main reasons for this, which are associated to the following essential challenges. The first challenge is the existence of the Dirac sea, in which the negative-energy spectrum is not bound from the bottom, in a relativistic Hamiltonian. In most of relativistic treatments, the densities are constructed with the occupied states in the Fermi sea only, explicitly excluding by hand all the states in the Dirac sea. This is referred to as the no-sea approximation, and is equivalent to neglecting the vacuum polarization effect [7]. In the no-sea approximation, the ground-state solution always corresponds to a saddle point on the energy surface, while in the non-relativistic case the ground state corresponds to the absolute minimum. It prevents a direct application of the imaginary time evolution for the coordinate-space calculations to the relativistic systems. That is, an iterative solution inevitably dives into the Dirac sea during the imaginary time evolution, leading to a divergence of the solution [48][49][50][51]. This problem is known as variational collapse. The second challenge in the relativistic calculations is the fermion doubling problem [52][53][54][55][56][57][58]. When one tries to solve a Dirac equation on discretized lattice in the real space, one tends to obtain spurious solutions with rapid oscillations as a function of coordinate, despite that the corresponding expectation value of energy is small. This problem arises as long as the first derivative in a Dirac Hamiltonian is approximated by a finite difference, and has been in fact well known in the field of lattice quantum chromodynamics (QCD) [52,53].
In this paper, we propose a novel and practical method for the CDFT calculations on 3D lattice, overcoming the problems of the variational collapse and the fermion doubling. For the first problem, we employ the inverse Hamiltonian method, which have been developed in the previous publication [51]. While in Ref. [51] we treated only the radial coordinate for a spherical relativistic Hamiltonian, in this paper we extend it to a 3D space. For the second problem, we extend the method of Wilson fermions, which has been widely adopted in the lattice QCD calculations [52,53]. With these strategies, we realize for the first time the relativistic calculations on 3D lattice without assuming any spatial symmetry.
The paper is organized as follows. In Sec. II, we will briefly recall the relativistic point-coupling model for the CDFT calculations. In Secs. III and IV, the two major difficulties, that is, the variational collapse and the fermion doubling, in the coordinate-space CDFT calculations will be explained in detail. We will provide and discuss the strategy to resolve each of these difficulties. In Sec. V, we will numerically apply our method firstly to a spherical nucleus 16 O and check the validity of the calculations. We will then apply the method to deformed nuclei 24 Mg and 28 Si and discuss its applicability. We will summarize the paper in Sec. VI.

II. RELATIVISTIC POINT-COUPLING MODEL
The relativistic variant of DFT, referred to as CDFT [4,[7][8][9][10][11], has been employed as widely as the nonrelativistic DFT in the studies of nuclear structure. A covariant energy density functional is obtained from a Lorentz invariant effective Lagrangian, which describes an effective interaction among nucleons via meson ex-changes. The relativistic treatment of nuclei has several characteristic features mainly due to the Lorentz invariance, which is the most important underlying symmetry of QCD. The CDFT has attracted much attention during the past decades because of its success in many aspects in nuclear physics. For instance, the saturation mechanism of nuclear matter [59,60], the large spin-orbit splittings which yield the shell structure of nuclei [60,61], and the origin of the pseudospin symmetry [62][63][64] are consistently understood with CDFT as the consequences of a delicate balance between the large attractive scalar and repulsive vector fields in nuclear medium. The timeodd components in the mean field, which are relevant to discussions of, e.g., rotating nuclei and nuclear magnetic moment, are entirely fixed by the Lorentz symmetry [65][66][67][68][69], while those in the non-relativistic models may have some ambiguity [70].
Throughout this paper we employ a relativistic energy density functional based on the relativistic point-coupling model. This is a zero-range interaction analogous to the non-relativistic Skyrme interaction. The effective Lagrangian is given by where is the free nucleon part, with ψ and m being the nucleon field and nucleon mass, respectively. The electromagnetic part L em is given by where A µ is the electromagnetic field and F µν = ∂ µ A ν − ∂ ν A µ is the electromagnetic field strength tensor. The nuclear interaction L int characterizes the relativistic point-coupling model. It consists of the four-fermion couplings L 4f , derivative terms L der , and higher-order terms L hot [71]: Here, L 4f is the leading-order of zero-range approximation to the meson-exchange interaction [72], while L der is the next-to-leading-order terms, which simulate the finite range of the meson-exchange interactions, The L hot corresponds to the self-couplings of the scalar and vector mesons, which introduce a density dependence into the N -N contact couplings, (7) Notice that the indices of the coupling constants, S, V , T S, and T V , denote the spin-isospin properties of the vertices, and they correspond to the exchanges of the σ, ω, ρ, and δ mesons, respectively.
Under the assumption of time-reversal symmetry, the energy density functional is obtained from the Lagrangian with the Hartree and the no-sea approximations as where K runs over the four spin-isospin channels, and the corresponding densities read Here ψ i (r) is the wave function for the i-th nucleon. The relativistic Hartree equation, or the relativistic Kohn-Sham equation, is obtained by taking a variation of the energy functional in Eq. (8) with respect to the single-particle wave functions as with where ρ is the proton density. These equations are solved self-consistently to obtain the ground state of atomic nuclei. Finally, one obtains the total binding energy as where the center-of-mass energy E CM is calculated by taking the expectation value of the kinetic energy for the center-of-mass motion of the whole nucleus with respect to the many-body ground-state wave function as [73]

A. Variational collapse
Our aim in this paper is to carry out three-dimensional mesh calculations with the energy density functional given by Eq. (8). To this end, we have to deal with the two challenges mentioned in Sec. I. In this section, we consider the first problem, i.e., the variational collapse, and discuss the practical solutions for that.
The variational principle is a simple but powerful guiding principle to find approximate solutions in nonrelativistic quantum mechanical problems. According to the variational principle, one obtains a better solution by minimizing the energy as much as possible. The imaginary time method, which has been successfully employed in 3D coordinate-space calculations for the nonrelativistic systems, is entirely based on such a variational principle [33]. That is, the variational principle guarantees that the evolution in imaginary time, where {|ψ k } denotes a set of single-particle wave functions, decreases the total Hartree-Fock energy as a function of τ , and eventually leading to a self-consistent solution.
In contrast, in the relativistic systems, the imaginary time evolution inevitably breaks down due to the presence of the Dirac sea states below the Fermi sea. If one naively applies the imaginary time method, the iterative solution dives into the continuum in the Dirac sea, which has been numerically confirmed in Ref. [50]. This occurs since the imaginary time evolution seeks for the lowest single-particle states. This is not what one wants, since only the lowest states in the Fermi sea are required in usual mean-field calculations with the no-sea approximation, whereas the Dirac sea states have to be explicitly excluded.
Such a breakdown of variational calculations, that is, the energy minimization, for relativistic systems is called variational collapse. The variational collapse problem has long been known and discussed in the field of relativistic quantum chemistry [74][75][76][77][78], mainly in the case of basis-expansion calculations. It has recently been discussed also in the context of nuclear physics [48][49][50][51], in connection to a realization of CDFT calculations on 3D coordinate-space representation.

B. Inverse Hamiltonian method
In order to avoid the variational collapse in the relativistic calculations, we follow the idea of Hill and Krauthauser [77], which is based on the variational principle for operator 1/(h − W ), where W is a real number not equal to any of the eigenvalues of h. As shown in Fig. 1(a), the ordinary variational principle is not applicable to a relativistic Hamiltonian because it has a negative-energy spectrum down to negative infinity as well as a positive-energy spectrum up to positive infinity. States between the two continua are discrete bound states. On the other hand, in the spectrum of the inverse Hamiltonian, 1/(h − W ), the two continua come to the middle of the spectrum while the bound states come to the two ends, see Fig. 1(b). The positive (negative)energy bound states come to the top (bottom) of the spectrum when the constant W is set between the Fermi sea and the Dirac sea.
Let us label the energy eigenstates of the Hamiltonian with an integer k according to the energy, so that the energy eigenvalues are denoted as where k > 0 (k < 0) corresponds to the Fermi (Dirac) sea solutions, with |φ k being the eigenstate associated to ǫ k . For simplicity of notations, here we have treated the continuum states as discrete states. Since the spectrum is bound both from the above and the below, we have the rigorous variational principle expressed as  iterative solution of Dirac equations in the coordinatespace representation. In this method, an initial state |ψ (0) , which is not an eigenstate of the Hamiltonian, is evolved with an operator e T /(h−W ) as Note that T has a dimension of energy and has nothing to do with time. The initial state |ψ (0) can be always expanded formally with the true eigenstates of the Hamiltonian h. In the evolution with e T /(h−W ) , all the negative-energy eigenstates contained in the initial state will damp away because the exponent is all negative. On the other hand, all the positive-energy states will grow up exponentially since the exponent is all positive, among which the state closest to W in the Fermi sea, |φ 1 , grows up most rapidly. Therefore, in the limit of T → ∞, the wave function converges to the lowest state in the Fermi sea, i.e., Moreover, if one takes T → −∞, the wave function converges to the highest state in the Dirac sea. Hereafter we only consider a positive T . With this inverse Hamiltonian method, one can obtain exclusively the positiveenergy solutions, which are usually of relevance to nuclear structure calculations. One can obtain the higher-energy states by starting with a set of initial wave functions and orthonormalizing them during the evolution by the Gram-Schmidt method. Alternatively, the iterative solution could also converge to |φ 2 by setting the shift parameter W to be between ǫ 1 and ǫ 2 . This scheme, however, may have a problem since the eigenvalues of h are not known a priori, and 1/(h−W ) becomes singular when W happens to be equal to one of them.
In practical calculations, one may cut T into steps by ∆T , and carry out the evolution iteratively. That is, in every step, the exponential function e ∆T /(h−W ) is expanded to the first order of ∆T . The iterative wave function at T = (n + 1)∆T is then given with the wave function at the previous step as Since this evolution operator is not unitary, one has to normalize the wave function at each step. One can show that the expectation value of the inverse of the Hamiltonian monotonically increases. In order to demonstrate this, let us assume that the step size ∆T is sufficiently small so that its second order is negligible. To the first order of ∆T , one obtains In the last line, the ∆T term is positive definite, since it is proportional to the dispersion of h −1 . One therefore finds, That is, the expectation value of h −1 increases from the step n to the step n + 1. This is a natural and reasonable consequence of the variational principle for the inverse of Hamiltonian.
Notice that, in contrast to h −1 , the behavior of the energy expectation value, h , is not necessarily monotonic: where the sign of the ∆T term depends on the property of ψ (n) , and it can be either positive or negative. Notice also that both in Eqs. (20) and (22), the ∆T term in the last line converges to zero as the iterative wave function converges to an eigenstate of the Hamiltonian.
In Ref. [51], it has been shown that the inverse Hamiltonian method works well for eigenvalue problems of a radial Dirac equations with a given spherical potential.
This method is relatively straightforward to be applied to not only the Dirac equation but also other eigenvalue problems with unbound operators, which is one of its advantages over some other methods for variational collapse [48-50, 75, 76]. In this context, we have successfully applied the inverse Hamiltonian method to non-relativistic Hartree-Fock-Bogoliubov equations with spherical potentials [79].
It is straightforward to extend the inverse Hamiltonian method to the self-consistent calculations. To this end, we simply replace the imaginary time evolution for the non-relativistic mean-field calculations with the Tevolution given by Eq. (17). That is, 1. Prepare an initial set of single-particle wave functions {ψ 3. Construct the single-particle Hamiltonian h (0) , which is a functional of the density ρ (0) .
4. Generate a new set of the single-particle wave functions, {ψ k }, by applying the T -evolution on each wave function:

Orthonormalize the set {ψ
The steps from 2 to 5 are iterated until the convergence is achieved. One may change the energy shift W (n) at each iteration.
Here we closely follow the discussion in Ref. [33] in order to show that the iteration leads to a self-consistent solution. The evolution of the wave functions by a small step ∆T can be approximated by With the Gram-Schmidt orthonormalization, one obtains the wave functions where the matrix elements of the inverse Hamiltonian is denoted as ψ lk . Since the change in the density matrix from the n-th step to the (n + 1)-th step is given by one obtains (27) which leads to Defining , one thus obtains (29) Therefore, the quantity always increases from iteration to iteration, which indicates that {|ψ }. This quantity will eventually converge to its maximum value (for a fixed value of W ), implying that Multiplying h − W from right and left on both sides of the above equation, one finally obtains at the convergence.

A. Dispersion relation and Fermion doubling
In the previous section, we have discussed one of the major problems, variational collapse, in the relativistic mean-field calculations on 3D lattice. In this section we will discuss the other problem, i.e., the fermion doubling. The fermion doubling problem has been well known in the field of lattice QCD, in which the first derivative in the action is replaced by a finite difference by discretizing the space-time within a box with the periodic boundary condition. Here we shall show how the problem arises for a static one-particle Dirac equation.
For simplicity let us consider a Dirac equation for a free particle in 1-dimensional space, where ψ(x) is a two-component spinor and Let us solve this equation by discretizing the coordinate x with a mesh size a within a box of size L. To this end, we approximate the derivative in the kinetic term with the 3-point differential formula, where x i ≡ a · i is the i-th mesh point. If either the periodic ψ(x + L) = ψ(x) or the anti-periodic ψ(x + L) = −ψ(x) boundary condition is imposed, the Dirac equation in the momentum space reads whereψ(k) is the Fourier transform of ψ(x). From Eq. (36) one obtains a dispersion relation which reduces to the ordinary dispersion relation of a relativistic particle in the continuum limit a → 0. In Fig. 2, we compare several approximate dispersion relations of a fermion on lattice with the exact one. The solid curve is the exact dispersion relation in the continuum limit show in Eq. (38), while the dashed curve shows the dispersion relation for a discretized Dirac equation given in Eq. (37). Moreover, the dotted and dash-dotted curves show the dispersion relations for a discretized lattice, obtained with the 5-and 11-point finite differential formulas, respectively. In the continuum limit, energy increases monotonically as a function of momentum as shown by the solid curve. In contrast, for the discretized fermions, a spurious minimum appears at the edge of the Brillouin zone, i.e., at the cut-off momentum k = π/a in the model space. As seen in the figure, this minimum does not disappear even if more accurate point-difference formulas are used, although the lower-momentum behavior is improved. If one solves the discretized Dirac equation, one thus obtains not only physical solutions with low energy and low momentum, but also spurious solutions with low energy and high momentum.
In the 1D case, one spurious state appears for each physical state, while in the 3D cases there are seven spurious states for each physical solution.
In the 3D cases, spurious minima appear on each of three axes of 3D momentum space at the boundary of the Brillouin zone, and the seven spurious solutions correspond to (x phys , y phys , z spur ), (x phys , y spur , z phys ), (x spur , y phys , z phys ), (x phys , y spur , z spur ), (x spur , y phys , z spur ), (x spur , y spur , z phys ), and (x spur , y spur , z spur ), where x phys and x spur denote the physical and spurious solutions in the x-direction, respectively, and similar for the y and z axes.

Wilson fermion
In order to eliminate the spurious fermions, Wilson introduced a term proportional to p 2 into the action, which yields an additional contribution to the dispersion relation. It can separate the energy of the spurious states on the high-momentum side from that of the physical states on the low-momentum side [52,53]. This method, referred to as the Wilson fermion, has been widely used in lattice QCD calculations.
We will demonstrate how it works for a static 1D Dirac equation. With the 3-point formula for the kinetic energy, the action of an 1D Hamiltonian on a wave function ψ(x) for a Wilson fermion is defined as The last term, the so-called Wilson term, has been added to the Hamiltonian. Here R is a dimensionless free parameter, and is referred to as a Wilson parameter. The corresponding Hamiltonian in the continuum limit may be formally written as which can be straightforwardly extended to the 3D space as Notice that, in the continuum limit a → 0, the Wilson term vanishes as a higher order, and the original form of Hamiltonian is recovered. Since the Wilson term, −aRβ∆, is proportional to p 2 , this term lifts up the spurious minimum of the dispersion relation at the edge of the Brillouin zone. Thus, the energy of the spurious states is all pushed upwards. Even though the physical states are also affected by the Wilson term, the effect is much smaller than that for the spurious states, because the expectation values p 2 are much smaller for the physical states. Therefore, the doubling problem can be avoided if one takes the value of R so that all the spurious states are pushed away from the energy region relevant to the calculation, that is, the fermi energy in the case of mean-field calculations. In Appendix, we discuss in detail the behavior of the energy spectrum of spurious single-particle states as a function of R.
For the 1D Hamiltonian given by Eq. (40), the dispersion relation reads when the 3-point formula is adopted for the kinetic term. Notice that this also returns to the exact dispersion relation in Eq. (38) in the continuum limit. Figure 3 shows a comparison of the dispersion relations with and without the Wilson term. To this end, we use R=0.2. For the comparison, the figure also shows the dispersion relation obtained with the 11-point formula (the dash-dotted line) as well as that in the continuum limit (the solid line).
As one can see, the spurious minimum at k = π/a disappears due to the Wilson term. At the same time, the dispersion relation deviates from the exact dispersion relation also on the low momentum side, although the effect vanishes at k = 0. That is, the Wilson term influences not only the spurious states but also the physical states, as may have been expected.
For instance, considering a = 0.8 fm as a typical mesh size for practical 3D calculations, the fermi momentum of nucleons k F ≃ 1.36 fm −1 corresponds to 0.35π/a. This implies that the dispersion relation is required to agree well with the exact one up to k ≃ 0.35π/a. Notice that the fermi energy of nucleons is ǫ F = k 2 F 2 /2m ≃ 38MeV ≃ 0.04m, thus, it is sufficient for the mean-fieldlevel calculations that the upward shift of the unphysical minimum by the Wilson term is 0.04m or larger.
In order to see clearly the effect of the Wilson term in the dispersion relation, we have chosen the value R = 0.2 rather arbitrarily in Fig. 3. This value may be too large given that the shift of the minimum is comparable to m, much larger than a desired value 0.04m. We have, however, confirmed that the physical solutions are still affected by the Wilson term even we take a smaller value of R.

High-order Wilson term
It is desired to find a prescription which affects the physical solutions as small as possible and at the same time affects the spurious solutions as much as possible.
In the field of chemistry and solid-state physics, other prescriptions have also been proposed than the Wilson fermion [56][57][58]. They are, however, not easy to apply when one employs an accurate differential formula, e.g., the 11-point formula, for the kinetic term. Therefore, we here attempt to extend the Wilson fermion, which can be applied in a straightforward way even if the higher-point differential formula is adopted.
To this end, we increase the power of p in the Wilson term as We show in Fig. 4 the dispersion relations for the 1D Hamiltonian obtained with the prescription given by Eq. (43). The dashed curve shows the result with the original Wilson fermion, corresponding to (n, R) = (1, 0.2) with the 3-point differential formula as in Fig. 3. The dotted and dash-dotted curves show the improved dispersion relations obtained with (n, R) = (4, 0.003) together with the 3-point and 11-point differential formulas, respectively. Evidently, the improved (n > 1) Wilson term leads to both a good agreement with the dispersion relation in the continuum limit in the low-k region and an efficient separation of the spurious states from the physical states in the high-k region.

V. SELF-CONSISTENT CALCULATIONS
In the previous sections, we have introduced the two methods to overcome the variational collapse and fermion doubling problems. We now combine them together and carry out the self-consistent relativistic mean-field calculations in the 3D coordinate-space representation. To this end, we first perform a self-consistent calculation with the Wilson term using the inverse Hamiltonian method with an operator e T /(hW −W ) , i.e., After the convergence is achieved, we correct the energy by the first-order perturbation theory, that is, the expectation value of the Wilson term is subtracted from the single-particle energies as well as from the total energy, and whereǫ k andẼ are the converged single-particle and total energies with the Wilson term, respectively. In Eq. (46), the subscript k in β k and ∂ ki means that they act on the k-th single-particle state. For simplicity, we do not correct the single-particle wave functions.
In the following, we present the results of benchmark calculations in order to show the validity of this strategy. To this end, we perform calculations in a cubic numerical box with L = a × (N − 1) long, where a is the mesh size and N is the number of mesh points taken along each direction. We shall vary the box size L and mesh size a, and examine the convergence of the results. In the inverse Hamiltonian method, we need to invert a large sparse matrix of the single-particle Hamiltonians in the coordinate-space representation. For this purpose, we employ an iterative solver, that is, the conjugate residual method [80], for large sparse linear systems. The Wilson parameter is fixed with (n, R) = (5, 0.00015) for all the calculations shown in this section. The first derivative of the kinetic term is approximated by the 11-point formula. We take the box boundary condition, although the discussion of the dispersion relation was made in Sec. IV with the periodic boundary condition. We have confirmed that the two boundary conditions yield almost the same results if the box size is large enough compared to the size of the system.
We first consider the 16 O nucleus without the Coulomb interaction as an example. We shall apply our method also to deformed nuclei. We employ the energy functional of PC-F1 interaction given in Ref. [71]. The selfconsistency of the iterative solutions is judged by a condition that the dispersions of energy of all the occupied single-particle states are smaller than 10 −8 MeV 2 .
Since we discretize the full 3D space without assuming any spatial symmetry, we must constrain the center of mass of the whole nucleus at the origin and also the principal axes aligned to the coordinate axes: x ≡ d 3 r xρ V (r) = 0, y = 0, and z = 0, (47) as well as xy = 0, yz = 0, and zx = 0.
When we draw potential energy surfaces, we in addition impose the constraints on the deformation parameters.
Here we define the nuclear multipole deformation parameters as with R = 1.2 × A 1/3 fm, where X ℓm is a real basis of the spherical harmonics, The ordinary quadrupole deformation parameters, β and γ, are related to α 2m as α 20 = β cos γ and √ 2α 22 = β sin γ. We employ the augmented Lagrangian method [81] for imposing these constraints.

A. Box size dependence
We first discuss the dependence on the box size L. Figure 5 shows the total energy of the 16  various values of L, without the Coulomb interaction and the center-of-mass correction, while the mesh size is fixed to be a = 0.8 fm. The dashed line shows the energyẼ before the correction of the Wilson term, while the solid is obtained after the correction,Ẽ − E W in Eq. (46). Thus the difference between the two lines is nothing but the the expectation value of the Wilson term. One can see that L ≃ 15 fm leads to a well converged result for 16 O, while E W remains less than 0.1 MeV. Figure 6 shows five single-particle energies as functions of L. Even though 16 O is a spherical nucleus, each multiplet of single-particle levels may split in energy in the 3D mesh calculation with the Wilson term for several reasons (see Sec. V D below). Here we have averaged the energies within the same multiplet. For the occupied states, i.e., the 1s 1/2 , 1p 3/2 , and 1p 1/2 states, one sees that L ≃ 15 fm is sufficiently large to have a convergence, which is consistent with the results for the total energy shown in Fig. 5. On the other hand, a larger box size is required for the less bound single-particle states, e.g., 1d 5/2 and 2s 1/2 , because of a long tail in these wave functions.

B. Mesh size dependence
We next discuss the convergence with respect to the mesh size a. To this end, we fix the box size to be L = a × (N − 1) ≃ 25 fm. For each value of a, we tune the number of mesh points N so that L is kept to be close to 25 fm. This value of L yields well converged results for all the bound single-particle energies as shown in Fig. 6, although L ≃ 15 fm may be sufficient for 16 O. Figure 7 shows the total binding energy of 16 O ob- tained as a function of a. For a comparison, the figure also shows by the dotted line the "exact" value obtained with a spherical code, in which the radial Dirac equation is solved by the Runge-Kutta method. One can clearly see that the total energy converges to the exact value as the mesh size decreases. Furthermore, the correction energy for the Wilson term becomes practically negligible for a = 0.6 fm or smaller. This is because the Wilson term is proportional to a 10 in this calculation and vanishes in the continuum limit.
Notice that the total energy increases as the mesh size a becomes smaller. This is opposite to what one would have expected from the variational principle for the non-relativistic systems. This is, however, not surprising since there is no variational principle in the usual sense for the relativistic systems. A similar phenomenon has been observed also in a relativistic mean-field calculation with basis expansion [13]. That is, the total energy of a nucleus increases as the number of basis for the lower component of the Dirac spinor increases. This is due to the variational collapse [76], that is, when a basis for the lower component is increased by one, one state will be added to the negative-energy single-particle spectrum and it will push upwards all the states in the Fermi sea. On the other hand, if a basis for the upper component is increased, a new state will appear on the top of the positive energy spectrum and it will push downwards all the other states. Figure 8 shows the single-particle energies as functions of a. As in Fig. 6, we have taken the average for each multiplet. It can be clearly seen that not only the total energy but also the single-particle energies approach closely to the exact results as the mesh size a decreases.

C. Accuracy of wave functions
We next demonstrate that the wave functions obtained with the present 3D code are as accurate as the total and the single-particle energies. To this end, we compare the total density of 16 O obtained with the 3D code to that with the spherical code. For the 3D mesh calculations, we use the box size of L = 25 fm, the mesh size a = 0.6 fm, and the 11-point formula for the kinetic term. This is the same set up as that for the results shown in Figs. 7 and 8. Figure 9 shows the total density distributions obtained with the two numerical codes, without the Coulomb interaction. The data points in the 3D result is dense for large values of r, because there are many mesh points which correspond to similar values of r. From the figure, one sees that the 3D code yields an almost identical result to the spherical code. In Fig. 10, we make a similar comparison but in the logarithmic scale. Notice that the range of r is wider than in Fig. 9. One can see that the 3D code yields a quite accurate result up to very large r, including the logarithmic tail. This is a particular nice feature of the real-space representation, whereas a ba-  sis representation may not be efficient in describing the asymptotic tail of the wave functions.

D. Artificial violation of the rotational symmetry
In the 3D mesh calculations with the Wilson term, the rotational symmetry is broken for several reasons. Firstly, the high-order Wilson term introduced in Sec. IV B explicitly violates the rotational symmetry for n ≥ 2. The original SO(3) symmetry is broken down to the octahedral symmetry O h . Secondly, the lattice discretization, finite-difference approximation, and the finite volume effect are also the sources of artificial violation of the rotational symmetry. That is, both a cubic lattice and a cubic box reduce the symmetry to O h , similarly to the high-order Wilson term. The effects of the high-order Wilson term, the discretization, and the finite-difference are expected to vanish in the continuum limit a → 0, while the effect of finite volume vanishes with an infinitely large box L → ∞. In the actual calculations, one needs to take a sufficiently small a and a large L so that these artificial symmetry breaking can be neglected.

Hexadecapole deformation
Let us discuss the violation of the rotational symmetry in the actual calculations, by taking again 16 O without the Coulomb interaction as an example. Since the ground state of this nucleus is spherical, all the deformation parameters are expected to vanish in principle. In the actual calculations, however, the octahedral symmetry O h , that is inherent in the 3D mesh calculation, may induce spurious hexadecapole deformation. Indeed, we find that the α 40 and α 44 in Eq. (49) are small but finite in the self-consistent solutions. Table I   L, while the mesh size is kept to be a = 0.8 fm. Meanwhile, the other hexadecapole deformation parameters, as well as the quadrupole and octupole deformation parameters, are all less than 5 × 10 −5 , which are not shown in the Table. The deformation parameters decrease as L increases, and eventually converge to a small value, which is practically negligible. The finite values which still remain even for large L are due to the effects of discretization, finite-difference, and the high-order Wilson term.
In Table II, we show the hexadecapole deformation parameters as a function of the mesh size a, while the box size is kept to be large enough, L ≃ 25 fm. We see that the deformation becomes smaller as a decreases, i.e., as it approaches the continuum limit, where the effects of discretization, finite-difference, and the high-order Wilson term vanish. The hexadecapole deformations due to the artificial symmetry breaking are already negligibly small around a = 0.8 fm, which is a typical mesh size in the 3D calculations.

Splitting of single-particle levels
Other quantities which are affected by the symmetry violation are the single-particle energies. With the group theory, it can be demonstrated that, for j ≥ 5/2, the (2j + 1)-fold degeneracy for the magnetic quantum num-  bers of the single-particle levels splits into several levels due to the reduction of symmetry from SO (3) to O h [82]. Similar level splitting of SO(3) multiplet on 3D lattice was also discussed recently in nuclear physics context by Lu et al. in Ref. [83]. We thus investigate here the splitting of the 1d 5/2 level in 16 O as a function of the mesh size a, while the box size is kept to be L ≃ 25 fm.
In Table III, we show the calculated single-particle energies of the levels corresponding to the 1d 5/2 and 1p 3/2 states in the SO(3) limit for several values of the mesh size a. As expected from the group theory, the 1d 5/2 level splits in energy, whereas the 1p 3/2 level does not within the digits shown in the table. Similarly to the case of the hexadecapole deformation, the level splitting becomes smaller as the mesh size a decreases, approaching to the SO(3) limit.
We note that, with the symmetry violation, the zcomponent of the single-particle angular momentum, j z , is no longer a good quantum number, and its actual expectation values are only approximately equal to the halfintegers. Since we prepare the initial single-particle wave functions on the lattice as spherical spinors, which are the eigenstates of j z , the j z values are approximately conserved during the iteration, in which the axial symmetry is violated only slightly. We also notice that, whereas the exact eigenstates of the mean-field potential with O h symmetry are the irreducible representation of the symmetry group, the resultant single-particle states in our present calculations are not the pure irreducible representations of O h . In this sense, the single-particle states with j ≥ 5/2 in the present calculations cannot completely be the eigenstates of the mean-field Hamiltonian. However, in practice, the convergence and the accuracy of the single-particle energies are satisfactory, as we have seen in this section.

E. Deformed nuclei
In the previous subsections, we have shown that our method yields accurate solutions in the 3D coordinatespace representation for the spherical 16  out the problems of variational collapse and fermion doubling. Let us now apply the method to deformed nuclei, 24 Mg and 28 Si. To this end, we take 20 × 20 × 20 mesh points with a = 0.8 fm. Figure 11 shows the potential energy surface of 24 Mg and 28 Si on the (β, γ) deformation plane obtained with the present 3D CDFT code. One can see that a prolate deformation for 24 Mg as well as an oblate deformation for 28 Si are successfully reproduced in these calculations. We emphasize that this is the first CDFT calculations on 3D lattice with constraints on deformation parameters.

VI. SUMMARY AND PERSPECTIVES
Because of the variational collapse and the fermion doubling, a 3D coordinate-space calculation with covariant density functionals has been impossible for a long time. In order to realize such calculations for the first time, we have proposed a novel and practical method to solve Dirac equations in the 3D coordinate space. To this end, we have introduced the two different prescriptions and combined them to overcome the variational collapse and the fermion doubling. For the variational collapse, we have employed a method based on the varia-tional principle for the inverse of a single-particle Hamiltonian, while for the fermion doubling, we have extended the method of Wilson fermion, which has been widely employed in lattice QCD calculations.
Using 16 O as an example, we have confirmed that our strategy provides accurate solutions for self-consistent mean-field calculations without the influence of the negative-energy spectrum and the spurious solutions of a discretized Dirac equation. We have shown with 24 Mg and 28 Si that this method is also applicable to deformed solutions in the (β, γ) deformation plane. Such calculations are important for discussions of deformation properties as well as fission properties of heavy nuclei. They can also be used as inputs for the generator coordinate method and five-dimensional collective Hamiltonian calculations in order to investigate spectroscopic properties of nuclei by taking into account the quantum fluctuation of shape degrees of freedom [11,84,85].
There will be many possible applications of our new code. In particular, it will enable one to i) study any complicated structure of nuclei with a single numerical code, ii) compare directly the results of the relativistic models to those of 3D mesh calculations with the non-relativistic models, and iii) provide reliable theoretical predictions with the relativistic models for unknown nuclei allowing symmetry-breaking solutions. We emphasize that our new relativistic 3D code allows one to study arbitrary shape of nuclei such as exotic deformations, halo structure, complicated shape relevant to nuclear fission, and even a cluster structure, without any restriction on the spatial symmetry and without significantly increasing the numerical cost. It also allows a straightforward extension of the finite-amplitude method [86,87] within the relativistic framework [88,89] for a study of nuclear excitations in deformed nuclei. We therefore believe that the method proposed in this work makes an important step to drastically extend the flexibility of the CDFT calculations in near future. In this appendix, we will show an example of the fermion doubling problem and examine our prescription given in Sec. IV in detail.
For this purpose, we consider a Dirac equation with the scalar S(x) and vector V (x) potentials in a 1-dimensional space, where ψ(x) = (ψ 1 (x), iψ 2 (x)) T is a two-component spinor, see Eq. (33). The mass is set to be m = 939 MeV/c 2 , while we take a Woods-Saxon type for the potentials V (x) and S(x), in which we use the parameters corresponding to 40 Ca given in Ref. [90]. We discretize the coordinate x with mesh size a = 0.1 fm with N = 400 mesh points and impose the box boundary condition. We use the 3-point differential formula for the kinetic energy term. The energy eigenvalues and expectation values of p 2 /2m for bound state solutions obtained with the inverse Hamiltonian method are summarized in Table IV. We have checked that the dispersion of the Hamiltonian, h 2 − h 2 , is close to zero for all the states shown in the Table. We find 5 pairs of bound states, and in each pair the two states have exactly the same energy. The expectation value of p 2 is extremely large for one of the states in each pair, implying that it is a spurious solution.  Table IV. Panel (a) to (f) correspond to k = 1 to 6 in the Table, while the ψ1 and ψ2 correspond to the upper and lower components of the wave functions, respectively.
In Fig. 12, we show the wave functions for the three lowest pairs. The left and right panels correspond to the physical and the spurious solutions, respectively. It is seen that the spurious states has the same amplitude as the corresponding physical states, but the sign of their wave functions alter at every mesh point. More precisely, if we denote a physical state by ψ (p) (x j ) = ψ It can be shown that, in the specific case of 3-point formula f ′ j = (f j+1 − f j−1 )/2a, the expectation values of kinetic term for the physical state and the corresponding spurious state are exactly the same. Their expectation values of potential term are also the same.
We have shown that each pair of physical and spurious states is degenerate in energy if one uses the 3point formula. What happens to the spurious states if one approximates the derivative in the kinetic term with more accurate differential formula? Figure 13  ingly the expectation value of p 2 /2m is quite large compared to a typical value of the kinetic energy of a nucleon ( 38 MeV). With the higher-order formulas, the spurious states are no longer degenerate to the corresponding physical states. For instance, with the 5-point formula, the spurious states in the 4-th and 5-th pairs are already pushed away to the continuum region, while the lower spurious states are still in the bound region. With the 7-point formula, the third unphysical state also goes up to the continuum region. The energy shift for the higher spurious states is larger than that for the lower spurious states. Therefore, it is important that the higher-point differential formulas for the kinetic term resolve the degeneracy and lift up the spurious states, although the energy shift may not be large enough to remove all the spurious states from the energy region of interest.
In Fig. 14, we show the wave functions of the lowest pair of states obtained with the 3-, 7-, and 11-point formulas. One can see that the amplitude of the spurious state deviates more from the physical state as the differential formula becomes more accurate.
Let us now switch on the Wilson term and investigate the change of the spectrum as a function of the Wilson parameter R. In Fig. 15, we show the energy spectrum of the Dirac equation as a function of R. Figure 15(a) is obtained with the original Wilson term (n = 1). The derivatives in the Wilson term and in the kinetic term are both approximated by the 3-point formula. Figure 15(b), on the other hand, shows the result with an high-order Wilson term (n = 4) evaluated with the 9-point formula, and the kinetic term is computed with the 11-point formula. One can see that the energy shift for the unphysical states is proportional to the Wilson parameter R in both cases. In the case of n = 1, all the unphysical states are pushed up to the continuum region around R = 0.008, while all of them are pushed up to the continuum already at about R = 0.0001 for n = 4. The high-order Wilson term is thus much more powerful than the original one. Figure 16 shows a comparison of the wave function for the lowest single-particle state obtained with the Wilson term to the exact one. The upper panel shows the case with the original normal Wilson term (n = 1 and R = 0.01) with the 3-point formula. On the other hand, the lower panel shows the result with an high-order Wilson term (n = 4 and R = 0.00015) computed by the 9-point formula. In both cases, the kinetic term is approximated by the 11-point formula. The "exact" wave function is obtained by solving the Dirac equation without the Wilson term by the inverse Hamiltonian method. The energy of this state with the Wilson term after the corresponding correction is ǫ =ǫ − ǫ W = −65.8726 MeV and ǫ = −65.8918 MeV for the cases (a) and (b), respectively, which can be compared with the exact energy ǫ = −65.8918 MeV evaluated with the 11-point formula for the kinetic term. Although the energies coincide between (a) and (b) only up to the first three digits, the difference in the wave functions is almost invisible in the scale shown in the figure. In order to quantify the deviation in the wave functions, one can compute the overlap probability | ψ ex |ψ W | 2 , where ψ ex is the exact wave function and ψ W is the approximate wave function obtained with the Wilson term. The deviation from unity is less than 10 −14 for both cases.
In short, the properties of the Wilson fermion and a spectrum of spurious states are summarized as follows, • Without the Wilson term, degenerate physical and unphysical states appear in pairs, if the first derivative in the kinetic term is approximated by the 3point formula. When a more accurate differential formula is used in evaluating the kinetic term, the degeneracies are resolved, i.e., the spurious states are pushed upwards while the physical states stay unchanged. However, the energy shifts for the unphysical states are not always large enough to remove all the spurious states up to the continuum.
• By switching on the Wilson term, all the spurious states are moved upwards by a similar amount of energy. The energy shifts are nearly proportional to the Wilson parameter R. • An increase of n in the high-order Wilson term makes the shifts of the spurious states drastically large while the shifts of the physical states remain small. Thus a Wilson term with n ≥ 2 has more effects on the spurious states and less effects on the physical states.
• The solutions obtained with the high-order Wilson fermion are close to the exact solutions to a sufficient accuracy both in the energy eigenvalues and in the wave functions.