Chebyshev-expanded covariant density functional from relativistic Brueckner theory

Introduction

Covariant density functional (CDF) theory provides a successful relativistic framework for finite nuclei, uniform nuclear matter, and neutron-star matter [1, 2, 3, 4]. In density-dependent CDFs, many-body correlations are encoded in effective meson–nucleon vertices that depend on the baryon density. Microscopic relativistic Brueckner–Hartree–Fock (RBHF) calculations provide guidance for this density dependence. A particularly relevant precedent is the work of Typel and Wolter [5], who used density-dependent $\sigma$, $\omega$, and $\rho$ vertices for nuclear matter and finite nuclei and compared their behavior with couplings inferred from Dirac–Brueckner self-energies. Subsequent work further developed the microscopic-to-local connection [6, 7, 8, 9]. Conventional exponential and rational parametrizations provide economical alternatives for organizing the density dependence.

The purpose of the present work is to construct a new covariant density functional, CE-ME1, through the reduction of modern full-Dirac-space RBHF information to local Hartree fields. The availability of all three Lorentz components $U_S$, $U_{0,0}$, and $U_V$ permits us to quantify the role of the spatial-vector self-energy in this reduction, to construct meson-resolved density-dependent Hartree couplings, and to test simultaneously their consistency with the microscopic RBHF equation of state (EoS). For a compact representation of the resulting density dependence of the vertices $g_i(\rho)$, we use the Chebyshev expansion

where $T_n$ are Chebyshev polynomials of the first kind, $i = \sigma,\omega,\delta,\rho$ labels the meson channel, and $m_i$ is its meson mass. In uniform matter the Hartree approximation implies that only the combinations $C_i$, not the individual $g_i$ and $m_i$, appear. For a fixed truncation order this spans the same polynomial space as an ordinary power series; the practical advantages are bounded basis functions on a finite interval, a hierarchical density-mode organization, and analytic density derivatives.

The microscopic input was taken from full-Dirac-space RBHF calculations [10, 11, 12] via a numerical reconstruction (digitization) of published results. The single-particle potential of neutrons $n$ and protons $p$ is decomposed as

The full-Dirac-space treatment determines all three Lorentz components without the ambiguities associated with positive-energy-only projections. A local Hartree functional has no explicit $\gamma\cdot\hat{p}$ term. A central point of this work is therefore to use the microscopic $U_V$ information explicitly before performing the Hartree reduction, rather than discarding it.

The self-energy data are available at the density grid $\rho= 0.16, 0.24, 0.32, 0.40,$ and $0.48\ \mathrm{fm}^{-3}$ as functions of the isospin asymmetry, whereas the RBHF energy-per-particle data extend down to about $0.08\ \mathrm{fm}^{-3}$. We use the 199 energy points inside $0.08 \leq\rho\leq0.48\ \mathrm{fm}^{-3}$. Thirteen additional energy points above $0.48\ \mathrm{fm}^{-3}$, extending to $0.50054\ \mathrm{fm}^{-3}$, are used only as near-boundary extrapolation diagnostics because no meson-resolved self-energy input is available there. Below $0.16\ \mathrm{fm}^{-3}$, the coupling shapes are constrained by the energy data and the smooth parametrization, not by meson-resolved self-energies.

Density-dependent Hartree framework

We consider nucleons coupled to effective isoscalar-scalar, isoscalar-vector, isovector-scalar, and isovector-vector meson fields. The interaction Lagrangian is

where the total baryon density is $\rho=\rho_n+\rho_p$. With $C_i(\rho)$ defined by Eq. (1), the direct Hartree self-energies in uniform matter are

where $\rho_s=\rho_{s,n}+\rho_{s,p}$, $\rho_{s3}=\rho_{s,n}-\rho_{s,p}$, and $\rho_3=\rho_n-\rho_p$. Introducing the isospin asymmetry $\alpha=(\rho_n-\rho_p)/\rho$, the neutron/proton densities are given by $\rho_{n,p}=[(1\pm\alpha)/2]\rho$ and the respective Fermi momenta are $k_{F,\tau}=(3\pi^2\rho_\tau)^{1/3}$. The scalar density is

where $M_\tau^*=M+\Sigma_{S,\tau}$ and $M$ is the bare nucleon mass. Because the vertices depend on the total baryon density, the common rearrangement self-energy is

Here a prime denotes differentiation with respect to $\rho$. It contributes to the chemical potentials and to the pressure through $\rho\Sigma_R$ term, but not explicitly to the energy density. Its inclusion is required for thermodynamic consistency.

Reduction of the full-Dirac-space self-energies

At the Fermi surface the microscopic RBHF dispersion relation can be written as

Introducing $\eta_\tau=U_{V,\tau}(k_{F,\tau})/k_{F,\tau}$, the reduced self-energies at this momentum can be written as

Equation (9) is an exact rewriting at the Fermi surface and identifies the quantities to be compared with local Hartree fields. It does not eliminate the momentum dependence of the RBHF self-energies or define a unique mapping of the microscopic energy functional onto a local Hartree functional.

The four Hartree coupling combinations are then extracted as

The microscopic RBHF self-energies are matched here to the direct Hartree self-energies. We exclude $|\alpha|<0.05$ from the pointwise isovector extraction, because $\rho_3,\rho_{s3}\to0$. The rearrangement contribution in Eq. (7) is not added to $U^{\mathrm{red}}_{0,\tau}$ in this mapping. In a density-dependent hadron field theory it is generated by functional differentiation and enters the effective-theory chemical potentials and pressure [13]. This direct-field convention does not imply equality of the full quasiparticle energies: the final self-consistent rearrangement field reaches $-25.6$ MeV at $\rho=0.40\ \mathrm{fm}^{-3}$ for $\alpha=0.4$.

At each density we use the median over the available asymmetry points as the density-only microscopic baseline. The central 68% interval over $\alpha$ is used as a diagnostic of the information lost in reducing the momentum- and asymmetry-dependent RBHF self-energy to four density-only Hartree channels.

The role of $U_V$ is numerically important. If it is omitted, the density-only least-squares residual of the scalar self-energy rises from about 1 MeV at $\rho=0.16\ \mathrm{fm}^{-3}$ to about 84 MeV at $0.48\ \mathrm{fm}^{-3}$. After the reduction (9), the corresponding residual (RMS over neutron/proton scalar self-energies over the sampled asymmetries at fixed density) remains below 4 MeV at all five densities. This estimate refers only to the microscopic projection, before thermodynamic calibration: final CE-ME1 scalar-field residuals can reach about 20 MeV. The strong high-density growth of the directly projected scalar coupling is largely an artifact of discarding the space-like vector component.

Chebyshev representation and microscopic baseline

We map the calibrated density interval onto $[-1,1]$ according to

where $\rho_{\min}=0.08\ \mathrm{fm}^{-3}$ and $\rho_{\max}=0.48\ \mathrm{fm}^{-3}$. The final functional uses the cubic representation

with $T_1=z$, $T_2=2z^2-1$, and $T_3=4z^3-3z$. The continuous Chebyshev orthogonality relation,

is used to organize the spectral expansion, but it does not make the discrete, nonlinear RBHF fit orthogonal: correlations among fitted coefficients depend on the data, imposed constraints, and the nonlinear form of the functional itself. We use the Chebyshev basis for its boundedness, hierarchical form, and analytic derivatives, not as an argument for statistical independence. The $UV$-reduced microscopic baseline is listed in Table 1. The couplings are smooth: $C_\sigma$ and $C_\omega$ decrease with density, $C_\delta$ increases moderately, and $C_\rho$ decreases. The asymmetry-induced relative half-width is below 1.2% for the isoscalar channels, while it reaches about 18% for $C_\delta$ at the lowest density and remains at the few-percent level over most of the interval.

$\rho$ [fm$^{-3}$]	$C_\sigma$	$C_\omega$	$C_\delta$	$C_\rho$
0.16	298.78	210.51	15.56	30.80
0.24	280.56	199.75	19.67	24.64
0.32	268.40	186.99	22.17	22.07
0.40	261.92	174.76	25.05	20.28
0.48	252.65	162.42	25.85	18.72

Thermodynamic calibration

The reduced self-energy projection provides the meson-resolved microscopic baseline, but matching self-energies at the Fermi surface is not identical to matching the RBHF energy functional. We therefore apply a second, deliberately low-dimensional thermodynamic calibration. The cubic functions are required to remain positive and monotonic in the same sense as the microscopic input over the full interval. This prevents the energy fit from exploiting weakly identified high-order scalar–vector mutual compensations.

The 16 cubic coefficients are fitted sequentially, with each of the 199 energy observations carrying equal weight. Five equalities and positivity/monotonicity conditions restrict the fit; the latter are not assumed statistically independent. We first calibrate $C_\sigma$ and $C_\omega$ to the symmetric-matter energy per particle $E/A$, imposing the microscopic saturation parameters

while pressure $P(\rho_{\mathrm{sat}}) = 0$ and incompressibility $K_\infty= 258$ MeV. Second, these isoscalar functions are held fixed and $C_\delta$ and $C_\rho$ are determined from the five asymmetric data sets. We define the quadratic symmetry energy by

where $S_2(\rho) = \frac{1}{2}\left.\partial^2(E/A)/\partial\alpha^2\right|_{\alpha=0}$. Its saturation values $J = S_2(\rho_{\mathrm{sat}})$ and $L = 3\rho_{\mathrm{sat}}S_2'(\rho_{\mathrm{sat}})$ are fixed to

These values are taken from the published Bonn-A full-Dirac-space RBHF calculation [12]. The curvature $K_{\mathrm{sym}} = 9\rho_{\mathrm{sat}}^2S_2''(\rho_{\mathrm{sat}})$ is not used in calibration. The energy residuals

and proximity to the smooth microscopic baseline select the remaining isovector freedom under the positivity and monotonicity constraints.

The final CE-ME1 coefficients are given in Table 2. They reproduce the fitted saturation quantities $E_{\mathrm{sat}}$, $K_\infty$, $J$, and $L$, and predict $K_{\mathrm{sym}} = 33.14$ MeV. The combined RMS deviation from the $N = 199$ RBHF energy points is 0.6127 MeV, with equal weight assigned to each point. Since the RBHF data carry no pointwise uncertainties, this RMS measures the discrepancy between the functional and the microscopic input rather than a statistical goodness of fit. The six asymmetry sets contain $N_\alpha= (36,39,30,35,30,29)$ points, with individual RMS deviations increasing from 0.413 MeV at $\alpha= 0$ to 0.745 MeV at $\alpha= 1$. For the sampled fixed-asymmetry sequences above saturation, $dP/d\rho$ remains positive and $c_s^2 \leq0.797$, indicating positive compressibility and subluminal sound speeds throughout the calibrated density interval.

$i$	$a_{i0}$	$a_{i1}$	$a_{i2}$	$a_{i3}$
$\sigma$	272.600466	-38.723320	2.119386	0.745273
$\omega$	195.097250	-30.905766	0.173863	2.168564
$\delta$	20.602444	7.137289	-1.594963	-0.079499
$\rho$	30.501194	-8.418124	3.878273	-0.862753

Figure 1 shows the RBHF energies and the final CE-ME1 curves. The residuals have a systematic density dependence rather than noise-like scatter; the largest fitted absolute deviation is about 1.76 MeV.

Figure 2 compares the final functions with the $UV$-reduced microscopic projection. At the five microscopic densities the final $C_\sigma$ and $C_\omega$ differ from the baseline by up to 6.3% and 2.6%, respectively; for the scalar field this corresponds to differences of order 20 MeV at the highest density. The final $C_\alpha$ remains within about 7.3% of the microscopic values. The largest adjustment occurs in $C_\rho$, reaching about 34% at $0.48\ \mathrm{fm}^{-3}$. The near-plateau of $C_\rho$ is sensitive to the selected shape constraints, but their quantitative influence has not been isolated by a relaxed refit.

Beta-equilibrated matter and comparison with DD-ME2

We next test the parametrization in cold charge-neutral $npe\mu$ matter, focusing on its equation of state and proton fraction, see Figure 3. As a phenomenological benchmark we include the same results obtained using the DD-ME2 parametrization [14]. Our implementation of the DD-ME2 model reproduces its standard saturation properties, $\rho_0 = 0.152\ \mathrm{fm}^{-3}$, $E/A = -16.14\ \mathrm{MeV}$, $K \simeq251\ \mathrm{MeV}$, $J \simeq32.3\ \mathrm{MeV}$, and $L \simeq51.3\ \mathrm{MeV}$ [4].

For CE-ME1 the proton fraction increases from $Y_p = 0.052$ at the RBHF saturation density to 0.152 at $0.48\ \mathrm{fm}^{-3}$. Along this sequence we find $\min(dP/d\rho) = 64.0\ \mathrm{MeV}$, $\max c_s^2 = 0.753$, where $c_s$ is the speed of sound. The baryonic curvature matrix $\partial\mu_i/\partial\rho_j$ is also positive definite on a grid spanning $0 \leq\alpha\leq0.9$ above saturation, with a minimum eigenvalue of approximately $305\ \mathrm{MeV\,fm^3}$. We do not impose mechanical stability below saturation, where a liquid–gas spinodal is physical, and no crust construction is attempted here.

For a direct microscopic comparison we also construct beta-equilibrated matter from the same RBHF energy-per-particle data used in the calibration. To ensure an accurate density derivative needed for the pressure, we represent the 199 RBHF energies in the calibration domain by a smooth tensor-product surface, using Chebyshev polynomials through order five in density and a cubic polynomial in $\alpha^2$. The representation has an RMS deviation of 0.097 MeV from the original RBHF energies. The sensitivity band in Figure 3 is obtained by varying the density order from four to seven and the $\alpha^2$ order from two to four. This band should be read as an interpolation/differentiation sensitivity, not as a statistical uncertainty.

Figure 3 shows the pressure and proton fraction for CE-ME1, DD-ME2, and the reconstructed RBHF EoS. This is an in-sample derivative/composition consistency check, not a prediction against independent RBHF data. The pressure agrees reasonably up to about $0.4\ \mathrm{fm}^{-3}$: at that density the pressures are 44.4 and 45.3 $\mathrm{MeV\,fm^{-3}}$, respectively. At the upper edge of the calibrated interval CE-ME1 gives $P = 94.6\ \mathrm{MeV\,fm^{-3}}$, compared with the central smooth RBHF value 87.3 $\mathrm{MeV\,fm^{-3}}$; the latter varies between approximately 82 and 90 $\mathrm{MeV\,fm^{-3}}$ under the interpolation changes described above. The CE-ME1 edge pressure is about 8.4% above the central reconstructed RBHF value and lies outside this interpolation-sensitivity band. The corresponding proton fractions are 0.152 and 0.155. DD-ME2 is appreciably stiffer over the same interval, reaching $P = 128.1\ \mathrm{MeV\,fm^{-3}}$ at $0.48\ \mathrm{fm}^{-3}$, and has a smaller high-density proton fraction, $Y_p = 0.111$.

Discussion

The main improvement of the CE-ME1 construction over a direct $U_S$–$U_0$ projection is the explicit use of the space-like RBHF self-energy. At $\rho= 0.48\ \mathrm{fm}^{-3}$ its magnitude is comparable with a sizable fraction of the Fermi momentum, so omitting it strongly distorts the scalar field inferred for a momentum-independent Hartree theory. Once $U_V$ is absorbed at the Fermi surface, all four microscopic coupling combinations become smooth and can be represented at cubic order.

A residual difference between the Fermi-surface self-energy mapping and the RBHF energy per particle is nevertheless expected. The full RBHF self-energies retain momentum dependence over the occupied Fermi sea, while the local Hartree functional is momentum independent. Consequently, matching the Lorentz components at one momentum is not an energy-functional identity. The thermodynamic calibration should therefore be interpreted as a controlled correction to the microscopic local projection, not as a new independent determination of the meson channels.

The cubic Chebyshev representation spans the ordinary cubic-polynomial space but offers bounded basis functions and analytic rearrangement derivatives. It does not remove correlations among the nonlinear fit coefficients. Higher-order constrained fits and convergence tests remain for future work; the reported asymmetry spreads and residuals are model-reduction diagnostics only.

The beta-equilibrium comparison (Figure 3) checks derivative and composition consistency of the thermodynamic reduction using the same microscopic energies. The increasing pressure difference near $0.48\ \mathrm{fm}^{-3}$ is consistent with the larger derivative sensitivity of the energy fit. The 13 unused near-boundary points, $0.48 < \rho\le0.50054\ \mathrm{fm}^{-3}$, have an RMS residual of 0.837 MeV and a maximum absolute residual of 1.710 MeV. Note that the relatively high $\rho_{\mathrm{sat}} = 0.188\ \mathrm{fm}^{-3}$ is inherited from Bonn-A RBHF and is retained in our analysis, i.e., we do not fit phenomenological saturation point. The DD-ME2 CDF, as a functional widely used in astrophysical applications, provides a phenomenological reference with a different saturation point and a markedly stiffer high-density equation of state.

The calibration is restricted to the density range $0.08 \le\rho\le0.48\ \mathrm{fm}^{-3}$. Extrapolating to nuclear surfaces, neutron-star crusts, or the central densities of massive stars would require additional microscopic input at low or high density. Several effects could improve the calibration of the functional: many-body contributions beyond the ladder summation, inelastic scattering at high energies, and short-range correlations. Momentum-dependent or Hartree–Fock extensions of the CDF developed here would also help, since they retain more of the microscopic input. Another route to a better functional is optimization against finite-nucleus data, which we leave for future work.

Conclusions

We have constructed the CE-ME1 density-dependent covariant functional from full-Dirac-space RBHF input. The microscopic $U_S$, $U_0$, and $U_V$ components are first reduced at the Fermi surface to effective scalar and time-like self-energies and then mapped onto the four-meson $\sigma$-$\omega$-$\delta$-$\rho$ Hartree channels. Including $U_V$ is essential: it removes the strong high-density growth of the scalar coupling found when the space-like component is neglected, and it reduces the residual of the local self-energy mapping to a few MeV. The density dependence of the couplings is represented by cubic Chebyshev polynomials on $0.08 \le\rho\le0.48\ \mathrm{fm}^{-3}$. A sequential isoscalar/isovector calibration to 199 RBHF energy-per-particle values reproduces the microscopic saturation point exactly and gives $K_\infty= 258.0$ MeV, $J = 33.10$ MeV, and $L = 65.20$ MeV, with a combined energy RMS of 0.613 MeV. The couplings remain positive and monotonic throughout the interval, and uniform matter is mechanically stable and causal above saturation. In beta equilibrium the CE-ME1 pressure and proton fraction agree reasonably with the same-data RBHF interpolation over most of the range, while the pressure discrepancy grows near $0.48\ \mathrm{fm}^{-3}$. The DD-ME2 CDF turns out to be substantially stiffer than the CE-ME1 counterpart at high density.

The curvature $K_{\mathrm{sym}} = 9\rho_{\mathrm{sat}}^2 S''(\rho_{\mathrm{sat}})$ is not used in the calibration. CE-ME1 predicts 33.14 MeV, whereas the published full-Dirac-space Bonn-A RBHF value is $-67.3$ MeV [12]. The value cannot reliably be recovered from second derivatives of our data interpolation, which are strongly order dependent. The discrepancy of $\simeq100$ MeV illustrates how sensitively this second derivative probes the low-order local representation; the restricted isovector shapes and the discarded momentum dependence may both play roles, without a separate quantitative attribution. We therefore regard $K_{\mathrm{sym}}$ as a sensitivity test of the reduction.

The present construction shows that a low-order Chebyshev representation can retain much of the microscopic Lorentz structure while providing a compact, thermodynamically consistent CDF. Additional full-Dirac-space self-energy points, especially above $0.48\ \mathrm{fm}^{-3}$, will be needed before CE-ME1 can be extended to high-density neutron-star applications. Such calculations may, however, fail to converge or to yield reliable results at these densities. In that case, one would have to resort to extrapolation schemes based on polytropes or speed-of-sound parametrizations.

In closing, it is instructive to compare the construction with that of Alford et al. [15], who also calibrate a relativistic mean-field functional using microscopic uniform-matter input. The two approaches differ in their microscopic input, functional form, and density range. Ref. [15] fits a density-independent $\sigma$-$\omega$-$\rho$ Lagrangian with nonlinear meson interactions to chiral-EFT pure neutron matter over $0.5 \le n/n_0 \le2$, supplemented by empirical saturation properties. Its extension to higher densities is determined by the adopted functional form and constrained by astrophysical observations. CE-ME1 instead employs density-dependent $\sigma$-$\omega$-$\delta$-$\rho$ couplings derived from full-Dirac-space RBHF self-energies, including $U_V$, and energies at six isospin asymmetries. Its microscopic calibration extends to $\rho= 0.48\ \mathrm{fm}^{-3}$, corresponding to about 2.55 times the RBHF saturation density, but no extrapolation beyond this range is attempted. The two approaches are therefore complementary: CE-ME1 is anchored in the Lorentz-resolved self-energies of relativistic RBHF theory [10, 11, 12], whereas the functional of Ref. [15] is constrained by nonrelativistic chiral-EFT calculations of neutron matter. The latter is extended beyond the microscopic density range within its relativistic mean-field framework, subject to additional astrophysical constraints. A comparable extension of CE-ME1 to higher densities will require confronting it with neutron-star observations, which we leave for future work.

Acknowledgments

This work was funded by Deutsche Forschungsgemeinschaft Grant No. SE 1836/6-1 and Polish NCN Grant No. 2023/51/B/ST9/02798.

Declaration of competing interest

The author declares that there are no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The $\alpha$-resolved RBHF $U_S$, $U_0$, $U_V$, and energy-per-particle input tables are numerical digitizations from previously published figures of Ref. [12]. The reconstructed tables, CE-ME1 coefficients, EoS and validation tables, and Python scripts are available from author upon request.