{'fraction_non_alphanumeric': 0.03775546934367876, 'fraction_numerical': 0.008398992120945487, 'mean_word_length': 4.800046392948272, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nuclear energy density functionals from empirical ground-state densities

A model is developed, based on the density functional perturbation theory and the inverse Kohn-Sham method, that can be used to improve relativistic nuclear energy density functionals towards an exact but unknown Kohn-Sham exchange-correlation functional. The improved functional is determined by empirical exact ground-state densities of finite systems. A test of the model and an illustrative calculation are performed, starting from two different approximate functionals, to reproduce the parameters and density dependence of a target functional, using exact ground-state densities of symmetric N=Z systems.


I. INTRODUCTION
Nuclear energy density functionals (EDFs) have been developed, more or less systematically, over the last two decades into a universal theoretical framework for the analysis of ground-state properties, low-energy collective excitations, and reaction dynamics of mediumheavy and heavy nuclei [1,2]. At present no other theoretical approach can be used to consistently describe diverse emergent nuclear phenomena with the same level of simplicity and accuracy, and at a comparable computational cost. Based on universal EDFs, a number of microscopic models, such as the random phase approximation [3][4][5], the interacting boson model [6,7], and the generator coordinate method [8][9][10], have also been designed to explore low-energy spectroscopy and large-amplitude dynamics on a quantitative level, and to calculate various details for astrophysical applications.
The exact but unknown nuclear EDF must be approximated by functionals of powers and gradients of ground-state nucleon densities and currents. Even though a great number of studies have been devoted to the microscopic formulation of a universal EDF framework (for a recent review, see Ref. [11]), the most successful nuclear EDFs are either semi-empirical or completely phenomenological. In a semi-empirical approach, one starts from a microscopically motivated ansatz for the nucleonic density dependence of the energy of a system of protons and neutrons. Part of the parameters of such a functional can be determined, at least qualitatively, by microscopic calculations of the energy of isospin symmetric and asymmetric infinite nuclear matter as a function of the nucleonic density (or Fermi momentum). The remaining parameters are usually adjusted to selected ground-state data, e.g., masses and charge radii, of an arbitrarily large set of nuclei. Fully phenomenological EDFs, for instance Skyrme, Gogny, and relativistic ones, usually take into account some empirical properties of nuclear matter at saturation, but all parameters are adjusted to ground-state data of finite nuclei.
The question we address in this work is how to improve a given functional towards the exact but unknown nuclear EDF. One could start, for instance, from a general expansion in powers of densities and currents and retain all terms allowed by symmetries up to a given order. Such a functional could be derived, in principle, from a microscopic theory (lowenergy QCD) that describes the underlying many-body dynamics. The problem, however, is that available low-energy nuclear data can only constrain a relatively small subset of terms and determine the corresponding parameters. Moreover, nuclear EDFs are "sloppy", that is, they generally exhibit an exponential range of sensitivity to parameter variations, and one finds many soft linear combinations of bare model parameters that are poorly constrained by data. This often indicates the presence of low-dimensional effective functionals associated with the relevant (stiff ) parameter combinations. In Refs. [12,13], we considered, in the context of nuclear structure, the interesting problem of a systematic construction of reduced low-dimensional functionals from a more complete but sloppy framework. Using methods of information geometry, it has been shown how to systematically construct effective EDFs of successively lower dimension in parameter space until sloppiness is eventually eliminated and the resulting functional contains only stiff combinations of parameters.
Instead of using low-energy data to reduce the complexity of a very general functional, one could also start from a relatively simple functional form and improve it towards the exact but unknown EDF. Such an expansion must, of course, be constrained by available data. In the spirit of density functional theory (DFT) [14][15][16], the empirical (exact) groundstate densities should determine the improved EDF. In fact, the inverse problem of DFT is formulated as a density-to-potential inversion that, starting from a given exact ground-state density, determines the true Kohn-Sham (KS) exchange-correlation potential [17][18][19][20][21][22][23][24][25]. The first application of the inverse KS method to nuclear EDFs has recently been reported in Ref. [26]. The inverse Kohn-Sham (IKS) method is often used to benchmark and test the accuracy of various approximate exchange-correlation functionals but, since its implementation depends on the exact density of a specific system, it does not provide direct information on the universal EDF.
The problem of improving the functional starting from exact ground-state densities has recently been addressed in Ref. [27], where a model was introduced based on a combination of the IKS method and the density functional perturbation theory (DFPT). The idea is to consider the difference, supposedly small, between the known functional and the exact but unknown EDF as a first-order perturbation. In Ref. [27], the method was successfully tested in benchmark calculations of systems of noble-gas atoms.
In this work, we apply the novel approach to atomic nuclei and, in particular, to relativistic nuclear EDFs. The principal reason for considering relativistic functionals is the fact that they automatically take into account the spin-orbit potential. The strong coupling between the orbital angular momentum and nucleon spin accounts for the empirical magic numbers and shell gaps. While in atomic systems the effect of the spin-orbit coupling is perturbative, in nuclei the energy spacings between spin-orbit partner single-nucleon states can be as large as the gaps between major shells. However, because the spin-orbit potential is a completely phenomenological addition to the non-relativistic KS potential, it cannot be determined from the ground-state density using the IKS method [26]. For relativistic EDFs, the spin-orbit potential emerges automatically with the proper strength as the constructive combination of the scalar and vector nucleon self-energies.
In Sec. II, we introduce the formalism of DFPT and IKS, and develop the model that will be used in this work. In Sec. III, a test case is discussed, and Sec. IV presents an illustrative calculation. Section V contains a brief summary and outlook for future studies.

A. Density functional perturbation theory
Here, we apply the method of Ref. [27] to improve a relativistic nuclear EDF, starting from a given empirical ground-state density. It is assumed that the unknown exact Hartreeexchange-correlation (interaction) functional can be written in the following form: where E (0) int [ρ] denotes the known functional that we wish to make better and λ is a small parameter. The main premise of this method, therefore, is that the difference between the exact functional E int [ρ] and the starting functional E (0) int [ρ] is small enough to be treated perturbatively. The correction will be determined by the exact ground-state density, using the IKS approach.
The Dirac KS single-nucleon equation, derived from the relativistic EDF, reads with the scalar and vector KS potentials, respectively, The corresponding scalar and vector densities ρ S,gs (r) = j∈occ ψ † j (r) βψ j (r) and ρ V,gs (r) = j∈occ ψ † j (r) ψ j (r) (4) are obtained from the self-consistent solutions of the single-nucleon Dirac KS equation in the no-sea approximation, which omits explicit contributions of negative-energy states to densities and currents [31,33,34] and, thus, the sums run only over the occupied positiveenergy single-nucleon orbitals.
The exact Dirac spinors, that is, the solutions of Eq. (2) for the exact EDF, can also be expanded to the first order in λ: where the first-order correction is orthogonal to the zeroth-order spinor The exact ground-state scalar and vector densities take the form: where the first term on the right-hand side denotes the densities derived by the known functional. We will assume that, given the exact densities ρ V(S),gs (r), one can use the IKS method to calculate the exact single-nucleon Dirac spinors ψ j (r) and energies j .
The ground-state energy can also be decomposed as follows: where the KS kinetic energy reads andt = α · p + βm. By expanding the Dirac spinors as in Eq. (5) and retaining only terms linear in λ, we obtain the following expression for the kinetic energy The second and third terms on the right-hand side of Eq. (9) denote the interaction (Hartreeexchange-correlation) contribution to the total energy. Using the expansion [Eq. (8)] for the ground-state densities up to terms linear in λ, one derives: S,gs , respectively, into the expansion for the interaction energy Eqs. (12a) and (12b), together with the expression for the kinetic energy Eq. (11), the ground-state energy reads One notices that the expression in square brackets represent the zeroth-order (unperturbed) Dirac Hamiltonian, that is The corresponding terms in Eq. (13) vanish because of the orthogonality relation [Eq. (6)] and, thus, a much simpler relation for the ground-state energy reads On the other hand, the ground state energy E gs = E kin + E int [ρ gs ], can be written in the following form where the Dirac KS equation has been used to eliminate the explicit contribution of the kinetic energy term, and j are the exact single-particle energies with the summation over occupied states. If we separate the zeroth-order and first-order terms of the exact interaction From Eq. (15), we can express the first-order correction to the interaction energy as a function of the zeroth-order ground-state densities where obviously E gs denotes the ground-state energy calculated with the known functional. Next, Eq. (17) for E gs is inserted in this expression and the following relation is obtained, interaction functional) are collected on the left-hand side The right-hand side of this equation depends only on the exact ground-state densities and the known functional E int . For given ground-state empirical densities, therefore, we can calculate all terms on the right-hand side, except the first term which is a sum of the exact single-particle energies. These energies are, of course, implicit functionals of the exact ground-state densities. One can, therefore, use the IKS method to calculate the singleparticle energies starting from given ground-state densities, a procedure that we describe in the following section.
In practical terms, one must assume a certain ansatz for the functional E , that will also include parameters to be determined from Eq. (20) for a choice of empirical ground-state densities. There is no guarantee, especially in the case of several undetermined parameters for the first-order correction, that the improved functional will reproduce the exact densities to a desired level of accuracy. The solution is an iterative procedure [27], in which the functional improved in the first iteration step is considered as the known functional for the next iteration, i.e., and the operation is repeated until the exact densities are reproduced by the solutions of the resulting n-th iteration Dirac KS equation to a desired accuracy.

B. Inverse Kohn-Sham method
Determining the KS potential for a given density presents an inverse problem. According to the Hohenberg-Kohn theorem, this inverse problem has a solution and the KS exchangecorrelation potential for a given system of interacting particles can be calculated starting from its ground-state density. Here, we perform the density-to-potential inversion in order to determine the single-particle energies that appear on the right-hand side Eq. (20). Starting from the single-nucleon Dirac KS equation (2), we rewrite the KS potentials: By multiplying Eq. (22) with ψ † j (r) from the left and summing over the occupied positiveenergy states, one obtains The scalar and vector densities that appear in this expression can also be combined in the following form: ρ + (r) = ρ V,gs (r) + ρ S,gs (r) and ρ − (r) = ρ V,gs (r) − ρ S,gs (r), so that If Eq. (22) is multiplied withψ j = ψ † j β, the following expression is obtained: (26) Finally, by adding and subtracting Eqs. (25) and (26), we derive from which the KS potentials V + and V − are expressed The set of IKS equations (28a) and (28b) can be solved iteratively. If we assume that the densities in the denominator are the exact (target) densities, and use Eqs. (27a) and (27b) in the numerator to define the densities and potentials for the n-th step, the resulting potentials in the (n + 1)-th step read In actual IKS calculations, we have modified an algorithm proposed in Ref. [21], and replaced Eqs. (29) and (30) with This algorithm was also used in the first nuclear IKS calculation with non-relativistic EDFs [26], and justified by the following argument. Equation (29) for the potential V + (r) = V (r) + S (r), which is the equivalent of the non-relativistic KS potential, has a simple interpretation: the potential is enhanced in absolute value in those regions where the density is larger than the target density, and reduced in regions where the density is smaller than the target density. However, this is what one expects for repulsive potentials (e.g., the Coulomb potential for electrons), whereas in the case of attractive potentials (e.g., the nuclear potential for nucleons) the opposite should happen. In Ref. [26], it has been shown that this issue can be avoided by adopting the modified algorithm of Eqs. (31) and (32).
Following Ref. [26], here we use the value of 1 MeV for both parameters γ + and γ − . In actual calculations we have encountered some stability issues for large values of the radial coordinate, due to small values of denominators in Eqs. (31) and (32) beyond the nuclear radius. This problem can be simply solved by introducing a cut-off radius r cut , and setting the potentials to zero for r > r cut . For the initial KS potential, a realistic Woods-Saxon potential [36] is used, and the Broyden mixing procedure [37] is employed to solve Eqs. (31) and (32). The convergence criterion used to halt the iterative IKS algorithm is defined in terms of the absolute variation of the potential, i.e.

III. A TEST CASE
The atomic nucleus is a complex quantum mechanical system with two types of constituent particles of spin one-half and, therefore, a general EDF will be a functional of isoscalar, isovector, and spin densities, as well as corresponding currents. Even in the simplest case of spin-saturated even-even nuclei, the EDF will depend on isoscalar and isovector densities. The problem, of course, is that accurate data exists only on charge (proton) densities and, thus, adjusting a general functional to empirical densities is not a straightforward procedure. In the particular case of relativistic EDFs that we consider here, the functional depends also on the Lorentz scalar single-nucleon density, which is not an observable. In the final section, we will discuss a possible approach that can be used to determine the scalar and isovector densities in an indirect way, but, for the purpose of testing the proposed method, here only N = Z systems without Coulomb interaction are considered. For such artificial nuclei, to demonstrate the relativistic IKS method with DFPT, we will use an existing relativistic EDF as the exact target functional, and apply the method developed in the previous section to improve different approximate functionals towards the target functional. In real nuclei, the exact functional is, of course, unknown and we will need more than ground-state data to determine the functional dependence on various nuclear densities.
For the exact target EDF, we will use the relativistic EDF DD-PC1 [35], for which the single-nucleon Hamiltonian readŝ where the scalar potential, vector potential, and rearrangement terms are respectively defined by In these expressions, m is the nucleon mass, α S (ρ), α V (ρ), and α TV (ρ) are density-dependent couplings for different space-isospace channels, δ S is the coupling constant of the derivative term, e is the electric charge, 1−τ 3 2 A 0 corresponds to the Coulomb interaction, and the single-nucleon densities ρ S (scalar-isoscalar density), ρ V (time-like component of the isoscalar current), and ρ TV (time-like component of the isovector current).
In addition to contributions from the isoscalar-vector four-fermion interaction and the electromagnetic interaction, the isoscalar-vector self-energy includes the rearrangement terms in Σ R that arise from the variation of the vertex functionals α S , α V , and α TV with respect to the nucleon fields in the vector density operatorρ V .
In a phenomenological construction of a relativistic EDF, one starts from an assumed ansatz for the medium dependence of the mean-field nucleon self-energies, and adjusts the free parameters directly to ground-state data of finite nuclei. Guided by the microscopic density dependence of the vector and scalar self-energies, the following practical ansatz for the functional form of the couplings was adopted in Ref. [35]: In the simplified case of N = Z doubly closed-shell nuclei without Coulomb interaction, there is no contribution of the isovector channel either. We will use four N = Z systems:  In the test case, we will assume for the known functional E (0) int [ρ V , ρ S ] a simple form that is actually a part of the DD-PC1 functional The values for the a  For the remaining unknown part of the functional E where α (1) s (ρ V ) and α (1) v (ρ V ) have the functional form of the density-dependent parts of the DD-PC1 couplings: with x = ρ V /ρ sat , and ρ sat = 0.152 fm −3 . As explained above, the parameters of the functional DD-PC1 were adjusted to reproduce the nuclear matter equation of state and the s , and b (1) v ) by using density functional perturbation theory and the IKS scheme, that is, using Eq. (20). Because the right-hand side of this equation is just a number that can be evaluated provided the exact single-particle energies and vector and scalar densities are known, a different finite system is needed for each parameter of the unknown functional. Since this is an illustrative test, we will employ three N = Z systems: 8 8 16, 28 28 56 and 50 50 100 to determine the constants b s and b (1) v , while d s and d v are again fixed to the DD-PC1 values. Note that, even though the problem has been simplified to a certain extent, nevertheless the test is far from being trivial. Namely, only three artificial systems are used to reproduce the values of parameters that were originally adjusted to the experimental masses of a large number of nuclei. Moreover, since our choice for the unperturbed functional is obviously not close to the exact target functional, it is far from obvious that a first-order perturbation method will determine the unknown parameters with sufficient accuracy. Hence, we repeat the calculation in several iterative steps, as described in the previous section [Eq. (21)], and at each step improve the values of b holds. Hence, hereafter, b v , the parameters converge to the values that correspond to the DD-PC1 target functional, denoted by the horizontal lines in Fig. 3. The results of this test demonstrate not only the feasibility of the IKS + DFPT method for nuclear densities, but also the convergence and accuracy of the iteration scheme.

IV. AN ILLUSTRATIVE CALCULATION
In the second example, we again use DD-PC1 as the unknown target functional, and the corresponding exact single-particle energies that appear on the right-hand side of Eq. (20) are obtained by the IKS method as described in Sec. II B. Also for the known functional int [ρ V , ρ S ] the simple form of Eq. (37) is adopted, that is we choose a polynomial form of the couplings α (1) s (ρ V ) and α (1) v (ρ V ): with x = ρ V /ρ sat , and ρ sat = 0.152 fm −3 . Therefore, we will examine whether the known functional E or expressed in terms of x = ρ V /ρ sat :

V. SUMMARY
In this study we have considered an interesting problem in the framework of nuclear energy density functionals, namely, how to improve a given functional towards an exact but unknown Kohn-Sham exchange-correlation functional. Based on the density functional perturbation theory and inverse Kohn-Sham method, a model has been developed that can be used to improve an approximate relativistic EDF.
Using the method introduced in Ref. [27] for non-relativistic functionals, and based on the density functional perturbation theory, we have derived Eq. (20) that is used to compute the first-order correction to an approximate zeroth-order functional. The input to this equation are the exact ground-state densities and the known functional, as well as the exact single- particle energies that are also implicit functionals of the densities. We then use the inverse Kohn-Sham (IKS) method to calculate these single-particle energies starting from given ground-state densities.
In practice, one must assume a certain ansatz for the first-order correction to the Kohn-Sham exchange-correlation functional, and use empirical exact ground-state densities of finite systems to determine the corresponding phenomenological parameters. In case the first-order functional does not reproduce the exact densities to a desired level of accuracy, the functional can be further improved in an iterative procedure, in which the first-order functional obtained in each iteration is considered as the known functional for the next iteration.
The model has been tested using the relativistic functional DD-PC1 as the exact target functional. A simplified form of DD-PC1 has been employed for the known functional.
Assuming for the first-order correction the same functional form as in the remaining part of DD-PC1, the method described above has been applied to determine three parameters. By employing only three finite N = Z systems, and with less than ten iterations, the resulting parameters of the first-order correction are found in excellent agreement with the original parameters of the functional. In a further illustrative calculation the target functional has been approximated by a different functional form, namely, a quadratic polynomial in the densities, determined by four parameters of the scalar and vector KS potentials. Even though the assumed density dependence differs from that of the target functional DD-PC1, nevertheless the model accurately reproduces the density-dependent coupling functions, as well as the target densities of four N = Z systems.
As noted in the introduction, the reason for considering relativistic functionals is that they automatically take into account the nuclear spin-orbit potential. The inclusion of this term in the nuclear KS potential is crucial to reproduce the empirical magic numbers and shell gaps, and yet in the non-relativistic case the spin-orbit term cannot be determined by the IKS method because there is no information on the corresponding density. The relativistic formulation does not provide a direct solution though. The reason is that the spin-orbital potential emerges as a constructive combination of the scalar and vector nucleon potentials, but the corresponding scalar density does not represent an observable. This brings us to the fact that accurate data exists only for charge (proton) densities, while in the IKS construction of the potential we need not only the isoscalar vector and scalar, but also the isovector densities. A possible solution would be to combine the model developed in this work, which utilizes empirical exact densities of finite nuclei, with the equations of state of isospin symmetric and isospin asymmetric nuclear matter. Namely, data on the proton vector densities in finite nuclei can be used together with the (microscopic) equations of state of nuclear matter to determine the isoscalar-scalar and isovector channels of the Kohn-Sham potential. Work along these lines is in progress.

Acknowledgments
This work has been supported in part by the Croatian Science Foundation under the project Uncertainty quantification within the nuclear energy density framework (IP-2018-01-5987). It has also been supported by the QuantiXLie Centre of Excellence, a project co-financed by the Croatian Government and European Union through the European Regional Development Fund -the Competitiveness and Cohesion Operational Programme