The Interrelated Roles of Correlations in the Nuclear Equation of State and in Response Functions: Application to a Chiral Confining Theory

Abstract

We study the role of short-range correlations, as well as pion and rho loops governing long-range RPA correlations, in nuclear matter properties and response functions. We use an adapted formulation of the Brueckner G-matrix approach to generate a pair correlation function satisfying the Beg–Agassi–Gal theorem, providing a natural cutoff to the loop integrals. We present results for the case of a relativistic chiral theory, including the effects of quark confinement and of the chirally broken vacuum in a version where parameters are directly connected to QCD observables or constrained by well-established hadron phenomenology. This provides a unified and coherent view of the nuclear matter equation of state and the effect of correlations on neutrino–nucleus scattering.

1. Historical Introduction

More than thirty years ago, the observed depletion of the nuclear structure function with respect to the free nucleon one, with the maximum effect observed when the Bjorken variable had a value of $x\simeq 0.6$ [1,2], was tentatively attributed to a possible modification of the nucleon structure in the nuclear medium. The reason was the inability of the standard mean-field Hartree–Fock approach [3] to reproduce this celebrated EMC effect, as discussed in [4]. This depletion is governed by the separation energy, $〈ϵ〉$, a negative quantity that is the opposite of the mean energy needed to remove a nucleon from the nucleus. For typical nuclei from ${}^{12}$C to ${}^{56}$Fe, the value $〈ϵ〉=-(40÷50)$ MeV gives a correct fit to the EMC effect. Hence, according to the Koltun sum rule, $-〈ϵ〉=〈t〉+2B$ [5], this necessitates a nucleon mean kinetic energy as large as $〈t〉\sim 35$ MeV, in strong disagreement with the conventional HF calculation giving ${〈t〉}_{HF}\sim 20$ MeV. The reason for this large value of kinetic energy is well known and is linked to correlations. Realistic calculations based on correlated wave functions have shown that the bare nucleon momentum distribution acquires a long tail beyond the Fermi momentum since the effect of short-range correlations is the depopulation of the Fermi sea. For instance, in the evaluation by Ciofi degli Atti et al. [6,7,8], for nuclei ranging from ${}^{12}$C to ${}^{56}$Fe, the kinetic energy per nucleon becomes larger, namely, $〈t〉\simeq 35$ MeV, and the absolute value of the separation energy also becomes larger, $-〈ϵ〉\simeq 50$ MeV. Thus, conventional nuclear effects are indeed able to reproduce the major part of the EMC effect. One can conclude that the independent particle picture is a very approximate view of the nucleus or that the Hartree–Fock scheme with bare nucleons is a poor approximation of the nuclear ground state. This conclusion seems surprising at first since we know that independent particle models based on the HF scheme accurately reproduce many basic properties of nuclei. The reason for the apparent discrepancy is due to the implicit identification of bare and dressed nucleons. In reality, Hartree–Fock nucleons should not be treated as bare nucleons but as nucleons surrounded by a polarization cloud, constituting quasi-particles (in the Landau sense), moving independently in a smooth mean-field potential. As pointed out in [9], such dressed objects couple differently from bare ones to an external (electroweak) probe.

A similar story occurred twenty years later in 2009/2010, when the measurement of the charged-current quasielastic-like (CCQE-like) neutrino–nucleus cross-section performed by the MiniBooNE collaboration [10,11] turned out to be in strong disagreement with the predictions based on the Fermi gas model (the model implemented in all Monte Carlo at that time; see, e.g., [11]), suggesting that scattering between a neutrino and independent nucleons strongly underestimates the measured cross-section. It was first proposed that this effect was due to a possible in-medium modification of the nucleon axial form factor, i.e., to an intrinsic modification of the bound nucleon structure. However, very soon, the explanation of this disagreement was given by the Lyon group [12]: nucleons are correlated via short-range correlations (SRCs) and meson exchange currents (MECs), which implies the possible ejection of a pair of nucleons from the scattering of the neutrino. This process is referred as the two-particle–two-hole (2p-2h) or n-particle–n-hole (np-nh) contribution to the cross-section. As an elaboration of the approach outlined in [9], we recently proposed in [13] a general technique based on a unitary transformation mapping of the Fermion operators relative to bare nucleons into quasi-particle operators relative to dressed nucleons. In short, this unitary transform is governed by the G matrix, which is the effective interaction felt by the dressed HF nucleon, with this G matrix $G(r)$ being very schematically related to the bare nucleon interaction $V(r)$ by $G(r)=V(r)f(r)$, where $f(r)$ is the pair correlation function.

In the following, we elucidate the influence of short-range correlations, which depopulate the Fermi sea and enhance the response functions, on the nuclear matter equation of state. The discussion is specifically illustrated on the basis of a microscopic model, inspired by QCD, describing the nucleon–nucleon interaction and its modification in the nuclear medium.

In the second section, we slightly reformulate the approach of [13], which allows the generation of the G matrix governing the interaction between dressed nucleons. We then propose a way to cause the resulting correlation function to vanish at zero distance, in agreement with the Beg–Agassi–Gal theorem [14,15,16]. The parameter $q_C$ appearing in this correlation function represents the inverse of the correlation hole radius induced by short-range correlations (SRCs). Its value can be fixed from the bare $NN$ interaction derived from a relativistic chiral theory, including the effect of quark confinement, developed in the recent years [17,18,19,20,21,22,23,24,25].

This approach, which we now call the chiral confining model, is the subject of Section 3. The original phenomenological version [18,20,21] described in Section 3.1 is based on a linear sigma model (L$\sigma$M) Lagrangian with standard Yukawa couplings of the various meson fields (scalar s, $\pi ,\rho ,\omega ,\dots$) to the nucleon but with three important physical pillars. First, the scalar field s, associated with the radial fluctuation of the chiral condensate, is identified with the sigma meson of relativistic theories [17,19]. Second, the effect of the quark substructure of the nucleon is reflected by its polarizability in the presence of the nuclear scalar field generating a repulsive three-nucleon force, providing an efficient saturation mechanism [18,20,21,25]. Third, the associated response parameters, namely, the nucleon scalar coupling constant, $g_S$, and scalar susceptibility, ${\kappa }_{NS}$, can be related to two chiral properties of the nucleon given by lattice QCD simulations, imposing severe constraints on their values [20,21,25]. This model has been applied in the past to the equation of state of nuclear matter and neutron stars, as well as to the study of the chiral properties of nuclear matter at different levels of approximation in the treatment of the many-body problem (RMF, Relativistic Hartree–Fock, or RHF, pion loop correlation energy).

However, it has been realized that the correct description of saturation properties compatible with lattice QCD constraints requires a large value of the scalar susceptibility, ${\kappa }_{NS}$, which is not compatible with realistic confining models of the nucleon. Hence, in Section 4.2, we present an enriched version recently proposed in [26], the NJL chiral confining model, which is based on the semi-bosonization of the Nambu–Jona-Lasinio model, as described in [23]. The main new feature of this improved framework is the replacement of the $L\sigma M$ effective potential, $V_{\chi ,L\sigma M}$, by a more repulsive NJL one, $V_{\chi ,NJL}$, considerably reducing the tension between the expected model values of the response parameters and the lattice QCD constraints.

A further step is accomplished in Section 3.3, where we build an effective Hamiltonian inspired by QCD that allows the simultaneous introduction, within some ansatz prescriptions, of a confining model for the nucleon, allowing the calculation of the response parameters, and an equivalent NJL model. This is the foundation of the QCD-connected chiral confining model, which basically depends on two genuine QCD quantities, namely, the string tension, $\sigma$, and the gluon condensate, ${\mathcal{G}}_2$.

Section 4 is devoted to nuclear matter calculations. We first summarize the values and the origin of the various parameters entering the (in-medium) nuclear interaction, trying to distinguish between those coming from the QCD-connected model (essentially the scalar and pionic sectors) and those coming from hadron phenomenology (essentially the vector sector), and we explain how the effect of short-range correlations is implemented in this context through the adapted G-matrix approach. We thus finally discuss the combined effect of short-range correlations and long-range RPA-like correlations (pion–rho loops) on nuclear matter properties and the depopulation of the Fermi sea, which is an indicator of the magnitude of the 2p-2h contribution to the neutrino–nucleus scattering cross-section.

2. Connecting Different Microscopic Models for Correlations

2.1. From Bare to Dressed Nucleons: Landau Quasi-Particles

Let us consider a nuclear system described by a Hamiltonian where the bare nucleon–nucleon interaction reduces to a two-body potential, written with standard notation as

$$H=\sum _{k_1,k_2}<k_1|T|k_2>a_{k_1}^{†}{a}_{k_2}+\frac{1}{4}\sum _{k_1{k}_2{k}_3{k}_4}<k_1{k}_2|\overline{V}|k_3{k}_4>a_{k_1}^{†}{a}_{k_2}^{†}{a}_{k_4}{a}_{k_3},$$

(1)

where a ($a^{†}$) is the annihilation (creation) operator relative to the bare nucleons. According to Goldstone [27], the correlated ground state of the nucleus, $|0>$, can be obtained from the state of uncorrelated bare nucleons, $|0{>}_{uncorr}$, by a unitary transformation as follows:

$$\begin{array}{ccc}\hfill |0>& =& U|0{>}_{uncorr}=U(a)\left(\prod _h{a}_h^{†}\right)|vac>=U(a)\left(\prod _h{a}_h^{†}\right)U^{†}(a)|vac>\hfill \\ & \equiv & \left(\prod _h{A}_h^{†}\right)|vac>,\hfill \end{array}$$

(2)

where $|vac>$ is the true vacuum state with a baryonic number equal to zero. The label h stands for hole states with a momentum below the Fermi momentum $k_F$. Hence, the correlated ground state $|0>$ can be obtained as a Slater determinant of dressed nucleons created by the $A_h^{†}$ operators. These operators are related to the bare creation operators, $a_h^{†}$, by a unitary transformation preserving the anticommutation relation:

$$A_h^{†}=U(a)a_h^{†}{U}^{†}(a)\iff a_h^{†}=U^{†}(A)A_h^{†}U(A).$$

(3)

Starting, as in references [28,29], with a unitary operator $U(A)=exp(S(A))$, with $S(A)$ being an anticommutation operator truncated at 2p-2h excitations,

$$S(A)=\frac{1}{4}(s_{p_1{p}_2{h}_1{h}_2}{A}_{p_1}^{†}{A}_{p_2}^{†}{A}_{h_1}{A}_{h_2}-hc),$$

(4)

one obtains, to the leading order in $S(A)$, the following result for hole (h) and particle (p) creation and annihilation operators:

$$\begin{array}{ccc}\hfill a_h& =& \sqrt{Z_h}{A}_h-\frac{1}{2}\sum _{h_2{p}_3{p}_4}{s}_{p_3{p}_4{h}_{2}h}{A}_{h_2}^{†}{A}_{p_4}{A}_{p_3}\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill a_p& =& \sqrt{Z_p}{A}_p-\frac{1}{2}\sum _{p_2{h}_3{h}_4}{s}_{p{p}_2{h}_3{h}_4}{A}_{p_2}^{†}{A}_{h_3}{A}_{h_4},\hfill \end{array}$$

(6)

where $Z_h$ and $Z_p$ are normalization factors that are fixed to ensure that $<0|\{a_k,a_k^{†}\}|0>=1$:

$$Z_h=1-\frac{1}{2}\sum _{h_2{p}_3{p}_4}{\left|s_{p_3{p}_4{h}_{2}h}\right|}^2,Z_p=1-\frac{1}{2}\sum _{p_2{h}_3{h}_4}{\left|s_{p{p}_2{h}_3{h}_4}\right|}^2.$$

(7)

To the same order, the ground state energy can obtained by a straightforward although tedious calculation as a functional of $s_{p_1{p}_2{h}_1{h}_2}$. These coefficients can be obtained variationally, i.e., by minimizing the ground-state energy, yielding:

$$\begin{array}{ccc}& & \left({ϵ}_{h_1}+{ϵ}_{h_2}-{ϵ}_{p_1}-{ϵ}_{p_2}\right)s_{p_1{p}_2{h}_1{h}_2}=<p_1{p}_2|\overline{V}|h_1{h}_2>\hfill \\ & & +\frac{1}{4}\sum _{p_1^{\prime },p_2^{\prime }}<p_1{p}_2|\overline{V}|p_1^{\prime }{p}_2^{\prime }>s_{p_1^{\prime }{p}_2^{\prime }{h}_1{h}_2}+\frac{1}{4}\sum _{h_1^{\prime },h_2^{\prime }}{s}_{p_1{p}_2{h}_1^{\prime }{h}_2^{\prime }}<h_1^{\prime }{h}_2^{\prime }|\overline{V}|h_1{h}_2>\hfill \\ & & \mathrm{with}\mathrm{single}-\mathrm{particle}\mathrm{energy}\text{:}{ϵ}_k=<k|T|k>+\sum _h<kh|\overline{V}|kh>\equiv <k|h_0|k>.\hfill \end{array}$$

(8)

The solution of this equation can be found as follows:

$$s_{p_1{p}_2{h}_1{h}_2}=\frac{<p_1{p}_2|\overline{G}(E={ϵ}_{h_1}+{ϵ}_{h_2})|h_1{h}_2>}{{ϵ}_{h_1}+{ϵ}_{h_2}-{ϵ}_{p_1}-{ϵ}_{p_2}},$$

(9)

where $G(E)$, the G matrix, is a two-body operator satisfying

$$\begin{array}{ccc}& & G(E)=V+V\frac{Q}{E-h_0}G(E)⟺\hfill \\ & & <k_1{k}_2|G(E)|k_1^{\prime }{k}_2^{\prime }>=<k_1{k}_2|V|k_1^{\prime }{k}_2^{\prime }>\hfill \\ & & +\sum _{p_1,p_2>k_F}<k_1{k}_2|V|p_1{p}_2>\frac{1}{E-{ϵ}_{p_1}-{ϵ}_{p_2}}<p_1{p}_2|G(E)|k_1^{\prime }{k}_2^{\prime }>,\hfill \end{array}$$

(10)

and Q is the Pauli-blocking operator projecting on particle states above the Fermi sea. A different derivation of this result, not based on a variational approach, is given in references [28,29]. It follows that the binding energy of nuclear matter is

$$E_0=\sum _{h<k_F}<h|T|h>+\frac{1}{4}\sum _{h_1{h}_2<k_F}<h_1{h}_2|\overline{G}(E={ϵ}_{h_1}+{ϵ}_{h_2})|h_1{h}_2>,$$

(11)

reflecting the fact that the dynamics of dressed HF nucleons is governed by the G matrix. The important point is that this G matrix, often identified with the effective interaction in the Hartree–Fock calculation, refers to objects, the dressed or renormalized nucleons, which can be seen as Landau quasi-particles that are definitely different from bare nucleons. If an energy-independent mean field is assumed, then the single-particle energy difference entering, for example, Equations (5), (6) and (9) involves only the kinetic energy difference between particle and hole states.

2.2. Response Function and Occupation Numbers

Using the correspondence between bare and dressed nucleon operators (Equations (5) and (6)), an explicit form of the genuine particle and hole spectral functions (i.e., related to the bare nucleons) can be obtained (see [13]):

$$\begin{array}{ccc}\hfill S_{k,k^{\prime }}^h(E)& =& \sum _n<0|a_k^{†}|n><n|a_{k^{\prime }}|0>\delta \left(E+E_n^{A-1}-E_0^A\right)=Z_k{\delta }_{k{k}^{\prime }}\delta \left(E-{ϵ}_k\right){\rho }_k\hfill \\ & +& \frac{1}{2}\sum _{p_2{h}_3{h}_4}\frac{<k{p}_2|\overline{G}|h_3{h}_4><k^{\prime }{p}_2|\overline{G}|h_3{h}_4>}{\left({ϵ}_k+{ϵ}_{p_2}-{ϵ}_{h_3}-{ϵ}_{h_4}\right)\left({ϵ}_{k^{\prime }}+{ϵ}_{p_2}-{ϵ}_{h_3}-{ϵ}_{h_4}\right)}\hfill \\ & & \delta \left(E-\left[{ϵ}_{h_3}+{ϵ}_{h_4}-{ϵ}_{p_2}\right]\right){\overline{\rho }}_k{\overline{\rho }}_{k^{\prime }}\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill S_{k^{\prime },k}^p(E)& =& \sum _n<0|a_{k^{\prime }}|n><n|a_k^{†}|0>\delta \left(E-E_n^{A+1}+E_0^A\right)=Z_k{\delta }_{k{k}^{\prime }}\delta \left(E-{ϵ}_k\right){\overline{\rho }}_k\hfill \\ & +& \frac{1}{2}\sum _{h_2{p}_3{p}_4}\frac{<p_3{p}_4|\overline{G}|h_2{k}^{\prime }><p_3{p}_4|\overline{G}|h_{2}k>}{\left({ϵ}_k+{ϵ}_{h_2}-{ϵ}_{p_3}-{ϵ}_{p_4}\right)\left({ϵ}_{k^{\prime }}+{ϵ}_{h_2}-{ϵ}_{p_3}-{ϵ}_{p_4}\right)}\hfill \\ & & \delta \left(E-\left[{ϵ}_{p_3}+{ϵ}_{p_4}-{ϵ}_{h_2}\right]\right){\rho }_k{\rho }_{k^{\prime }}.\hfill \end{array}$$

(13)

By integrating the spectral function, the occupation numbers are also obtained:

$$\begin{array}{ccc}\hfill n_h\equiv 1-\mathrm{\Delta }{n}^h=& =& 1-\frac{1}{2}\sum _{h_2{p}_3{p}_4}{\left|\frac{<p_3{p}_4|\overline{G}|h_{2}h>}{{ϵ}_h+{ϵ}_{h_2}-{ϵ}_{p_3}-{ϵ}_{p_4}}\right|}^2\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill n_p\equiv \mathrm{\Delta }{n}^p& =& \frac{1}{2}\sum _{p_2{h}_3{h}_4}{\left|\frac{<p{p}_2|\overline{G}|h_3{h}_4>}{{ϵ}_p+{ϵ}_{p_2}-{ϵ}_{h_3}-{ϵ}_{h_4}}\right|}^2,\hfill \end{array}$$

(15)

which is a well-known form, quoted, for instance, in ref. [30], Equations (2.31 and 2.32).

We can now obtain the response function to a probe (such as the vector boson related to charged current neutrino interaction), which transfers an energy-momentum to the nucleus $(\omega ,\mathbf{q})$ and couples to individual nucleons with time-independent operators: $\mathcal{O}(j)={\tau }_j^{\pm },({\mathbf{\sigma }}_j·\widehat{q}){\tau }_j^{\pm },{({\mathbf{\sigma }}_j\times \widehat{q})}^i{\tau }_j^{\pm }$1. This response is defined as

$$R(\omega ,\mathbf{q})=\sum _n|<n|\sum _{j=1}^A\mathcal{O}(j)e^{i\mathbf{q}·{\mathbf{x}}_j}{|0>|}^2\delta (\omega -E_n+E_0),$$

(16)

where $|n>$ and $E_n$ are the eigenstates and eigenvalues of the full nuclear Hamiltonian. In the factorization approximation, this response can be expressed in terms of the above spectral functions according to:

$$\begin{array}{ccc}\hfill R(\omega ,\mathbf{q})& =& \sum _{k_1^{\prime },k_2^{\prime },k_1,k_2}<k_1^{\prime }|{\mathcal{O}}^{†}{e}^{-i\mathbf{q}·\mathbf{x}}|k_2^{\prime }><k_2|\mathcal{O}{e}^{i\mathbf{q}·\mathbf{x}}|k_1>\hfill \\ & & {\int }_{-\infty }^{{\epsilon }_F}d{E}_1{S}_{k_1,k_1^{\prime }}^h(E_1){\int }_{{\epsilon }_F}^{\infty }d{E}_2{S}_{k_2^{\prime },k_2}^p(E_2)\delta (\omega -E_2+E_1).\hfill \end{array}$$

(17)

The indices k ($k=(\mathbf{k},s,t)$) stand for the basis of single-particle states, and $a_k$ and $a_k^{†}$ are the associated destruction and creation operators with discrete normalization $\left\{a_k,a_{k^{\prime }}^{†}\right\}={\delta }_{k{k}^{\prime }}$. The expression of the response function involves the known coupling of the external (electroweak) probe to the bare nucleon. The leading term (Equation (34) in [13]) is the pure p-h response quenched by Z factors. Correlations generate two additional contributions, the so-called correlation piece (Equation (35) in [13]) and the so-called polarization piece (Equation (36) in [13]) involving 2p-2h excitations at the origin of the enhancement of the neutrino–nucleus quasielastic cross-section seen by the MiniBooNE collaboration [10,11]. In addition, notice that, as explained in detail in [13], the delta resonance states can be incorporated as extra particle states in all parts of this approach.

An excellent indicator of the role of the correlations is provided by the sum rule:

$$S(\mathbf{q})=\frac{1}{A}\int d\omega R(\omega ,\mathbf{q})=\frac{4}{A}\sum _{\mathbf{k}}{n}_{\mathbf{k}}\left(1-n_{\mathbf{k}+\mathbf{q}}\right).$$

(18)

In particular, at zero momentum, this sum rule does not vanish. Indeed, we established the result in [13]:

$$S(0)\approx 2\mathrm{\Delta }n,$$

(19)

where $\mathrm{\Delta }n$ represents the mean depopulation (i.e., per nucleon) of the Fermi sea,

$$\mathrm{\Delta }n=\frac{1}{2A}\sum _{h_1,h_2,p_1,p_2}{\left|\frac{<p_1{p}_2|\overline{G}|h_1{h}_2>}{{ϵ}_{h_1}+{ϵ}_{h_2}-{ϵ}_{p_1}-{ϵ}_{p_2}}\right|}^2\equiv \kappa ,$$

(20)

which just coincides with the dimensionless wound integral of the Brueckner theory [31,32], which is the norm of the defect wave function averaged over all NN pairs. The defect wave function is the difference between the correlated and uncorrelated NN wave functions, and therefore, $\kappa$ is a measure of the importance of NN correlations, whose net effect is the depopulation of the Fermi sea. To show the effect of the correlation on the response function with a concrete example, we display in Figure 1 the result of the calculation (see [13]) for the sum rule in a toy model adjusted to reproduce the values obtained in references [6,7,8] for light nuclei such as carbon:

$$\mathrm{\Delta }n=0.15,<t>=\frac{3}{k_F^3}{\int }_0^{\infty }dk{k}^2\frac{k^2}{2M}{n}_{\mathbf{k}}=35\mathrm{MeV}.$$

(21)

with $k_F$ = 245 MeV or $\rho =0.8{\rho }_0$. We can check that $S(0)\approx 2\mathrm{\Delta }n=0.3$ with good accuracy and that the enhancement due to the correlations persists up to a momentum $q=350\mathrm{MeV}$.

2.3. The G Matrix and the Pair Correlation Function

To simplify the matter described below, we first omit the tensor force that we will reintroduce later perturbatively. We separate the total CM momentum and the relative momentum of the two-nucleon states and write the G-matrix elements in each given spin–isospin ($S,T$) channel as

$$\begin{array}{ccc}& & <{\mathit{k}}_1{\mathit{k}}_2|G^{ST}(E)|{\mathit{t}}_1{\mathit{t}}_2>=\frac{1}{\mathcal{V}}<\mathit{k}|G_K^{ST}(E)|\mathit{t}>{\delta }_{\mathbf{K},\mathbf{T}}\equiv \frac{1}{\mathcal{V}}{G}_K^{ST}(\mathit{k},\mathit{t};E){\delta }_{\mathbf{K},\mathbf{T}}\hfill \end{array}$$

(22)

$$\begin{array}{ccc}& & \mathrm{with}\text{:}{\mathit{k}}_1=\frac{\mathit{K}}{2}+\mathit{k},{\mathit{k}}_2=\frac{\mathit{K}}{2}-\mathit{k},{\mathit{t}}_1=\frac{\mathit{T}}{2}+\mathit{t},{\mathit{t}}_2=\frac{\mathit{T}}{2}-\mathit{t},\hfill \end{array}$$

(23)

where all the plane-wave states are normalized to one in a box with a volume $\mathcal{V}$. In view of the nuclear matter calculation, the entry energy will be the total energy of the two initial nucleons in the Fermi sea:

$$E={ϵ}_{\mathbf{K}/\mathbf{2}+\mathbf{t}}+{ϵ}_{\mathbf{K}/\mathbf{2}-\mathbf{t}}<2{ϵ}_F.$$

(24)

Similarly, the matrix elements of the bare two-body potential are written as:

$$\begin{array}{ccc}& & <{\mathit{k}}_1{\mathit{k}}_2|V^{ST}|{\mathit{t}}_1{\mathit{t}}_2>=\frac{1}{\mathcal{V}}<\mathit{k}|V^{ST}|\mathit{t}>{\delta }_{\mathbf{K},\mathbf{T}}\hfill \end{array}$$

(25)

$$\begin{array}{ccc}& & <\mathit{k}|V^{ST}|\mathit{t}>\equiv V^{ST}(\mathit{k}-\mathit{t})=\int d\mathit{r}{e}^{-i(\mathbf{k}-\mathbf{t})·\mathbf{r}}{V}^{ST}(\mathit{r}).\hfill \end{array}$$

(26)

Just to fix the idea, let us first consider the case of a potential only involving the sigma, omega, pion and rho exchange channels. In the simplest non-relativistic scheme, ignoring recoil correction, one has (with standard notations) for the most important spin–isospin states (i.e., the ones surviving in the orbital s state after antisymmetrization):

$$\begin{array}{ccc}& & S=0,T=1\mathrm{and}S=1,T=0\mathrm{states}:\hfill \\ & & V^{01}(q)=V_{\sigma }(q)+V_{\omega }(q)-\left(V_{\pi }(q)+2V_{\rho }(q)\right)=V^{10}(q)\hfill \\ & & V_{\sigma }(q)=-g_{\sigma }^2\frac{1}{q^2+m_{\sigma }^2}{\mathrm{\Gamma }}_S^2(q),V_{\omega }(q)=g_V^2\frac{1}{q^2+m_V^2}{\mathrm{\Gamma }}_V^2(q),\hfill \\ & & V_{\pi }(q)=-{\left(\frac{g_A}{2F_{\pi }}\right)}^2\frac{q^2}{q^2+m_{\pi }^2}{\mathrm{\Gamma }}_{\pi }^2(q),V_{\rho }(q)=-{\left(\frac{g_A}{2F_{\pi }}\right)}^2{C}_{\rho }\frac{q^2}{q^2+m_{\rho }^2}{\mathrm{\Gamma }}_{\rho }^2(q),\hfill \\ & & \mathrm{with}{C}_{\rho }={\left(\frac{F_{\pi }{g}_{\rho }(1+{\kappa }_{\rho })}{g_A{M}_N}\right)}^2,{\mathrm{\Gamma }}_{\alpha }(q)=\frac{{\mathrm{\Lambda }}_{\alpha }^2}{{\mathrm{\Lambda }}_{\alpha }^2+q^2},\hfill \end{array}$$

(27)

where q represents the modulus of the exchanged three-vector momentum. The form factors that mathematically regularize the UV (i.e., the short-distance) behavior of the potential are physically motivated by the compositeness and the finite size of the nucleon source of the emitted mesonic fields. In particular, they remove the contact piece of the central part of the (spin–isospin) $\pi$ and $\rho$ exchange and make the globally repulsive interaction at $r=0$ finite, transforming the hard core into a soft core. Of course, they do not affect the globally attractive long-range part of the interaction, dominated by the pion exchange and rho exchange in the strong rho scenario (${\kappa }_{\rho }\sim 6$) [33]:

$${\left[V^{01}\right]}_{longrange}(r)=-g_S^2\frac{e^{-m_{S}r}}{4\pi r}+g_V^2\frac{e^{-m_{V}r}}{4\pi r}-{\left(\frac{g_A}{2F_{\pi }}\right)}^2\left(m_{\pi }^2\frac{e^{-m_{\pi }r}}{4\pi r}+2C_{\rho }{m}_{\rho }^2\frac{e^{-m_{\rho }r}}{4\pi r}\right).$$

(28)

The combination of the repulsive short-range piece and the attractive long-range piece produces a pocket (see Figure 2) able to generate an approximately zero-energy bound state. Of course, distinguishing the singlet channel ($S=0,T=1$) from the triplet deuteron channel ($S=1,T=0$) necessitates the introduction of the tensor potential related to the difference between pion and rho exchange potentials and the time component of the Lorentz vector part of the rho exchange, with these two effects making the triplet interaction significantly more attractive than the singlet interaction.

In order to isolate the role of the short-range correlations, we introduce a mixed representation of the G matrix (see ref. [34]) according to:

$$\begin{array}{ccc}\hfill G_K^{ST}(\mathit{r},\mathit{t})& =& <\mathit{r}|G_K^{ST}(E)|\mathit{t}>=\int \frac{d\mathit{q}}{{(2\pi )}^3}{e}^{i\mathbf{q}·\mathbf{r}}{G}_K^{ST}(\mathit{q},\mathit{t};E)\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill V^{ST}(\mathit{r},\mathit{t})& =& <\mathit{r}|V^{ST}|\mathit{t}>=\int \frac{d\mathit{q}}{{(2\pi )}^3}{e}^{i\mathbf{q}·\mathbf{r}}<\mathit{q}|V^{ST}|\mathit{t}>=V^{ST}(\mathit{r})e^{i\mathbf{t}·\mathbf{r}}.\hfill \end{array}$$

(30)

This suggests the need to search for a solution for the G matrix in the form:

$$G_K^{ST}(\mathit{r},\mathit{t})=V^{ST}(\mathit{r})f_K^{ST}(\mathit{r},\mathit{t}).$$

(31)

It follows that the G-matrix equation is equivalent to an equation for the correlation function $f_K^{ST}(\mathit{r},\mathit{t};E)$:

$$\begin{array}{ccc}\hfill f_K^{ST}(\mathit{r},\mathit{t})& =& e^{i\mathbf{t}·\mathbf{r}}+\int d\mathit{y}{\left[\frac{Q}{e}\right]}_K(\mathit{r}-\mathit{y})V^{ST}(\mathit{y})f_K^{ST}(\mathit{y},\mathit{t})\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill {\left[\frac{Q}{e}\right]}_K(\mathit{r}-\mathit{y})& =& \int \frac{d\mathit{q}}{{(2\pi )}^3}{e}^{i\mathbf{q}·(\mathbf{r}-\mathbf{y})}\frac{P_{\mathbf{K}/\mathbf{2}+\mathbf{q}}{P}_{\mathbf{K}/\mathbf{2}-\mathbf{q}}}{{ϵ}_{\mathbf{K}/\mathbf{2}+\mathbf{t}}+{ϵ}_{\mathbf{K}/\mathbf{2}-\mathbf{t}}-{ϵ}_{\mathbf{K}/\mathbf{2}+\mathbf{q}}-{ϵ}_{\mathbf{K}/\mathbf{2}-\mathbf{q}}}.\hfill \end{array}$$

(33)

All these equations can be rewritten in an operator form,

$$\begin{array}{ccc}\hfill G_K^{ST}(E)& =& V^{ST}{F}_K^{ST}(E),G_K^{ST}(E)=V^{ST}+V^{ST}{\left[\frac{Q}{e}\right]}_K(E)G_K^{ST}(E)\hfill \\ \hfill F_K^{ST}(E)& =& 1+{\left[\frac{Q}{e}\right]}_K(E)V^{ST}{F}_K^{ST}(E)=1+{\left[\frac{Q}{e}\right]}_K(E)G_K^{ST}(E),\hfill \end{array}$$

(34)

in terms of operators acting in the space of relative motion and defined by their matrix elements:

$$\begin{array}{ccc}& & <\mathit{r}|1|\mathit{t}>=e^{i\mathbf{t}·\mathbf{r}},<\mathit{r}|V^{ST}|\mathit{t}>=e^{i\mathbf{t}·\mathbf{r}}{V}^{ST}(\mathit{r})<\mathit{r}|F_K^{ST}(E)|\mathit{t}>=f_K^{ST}(\mathit{r},\mathit{t})\hfill \\ & & <\mathit{r}|G_K^{ST}(E)|\mathit{t}>=G_K^{ST}(\mathit{q},\mathit{t})=F_K(\mathit{q},\mathit{t})V^{ST}(\mathit{r}),\hfill \\ & & <\mathit{x}|{\left[\frac{Q}{e}\right]}_K(E)|\mathit{y}>={\left[\frac{Q}{e}\right]}_K(\mathit{x}-\mathit{y};E).\hfill \end{array}$$

(35)

The correlation function can be expanded in partial waves, and one can define a correlation function for each two-body channel identified by the quantum numbers $STl$:

$$\begin{array}{ccc}\hfill f_K^{ST}(\mathit{r},\mathit{t})& =& \sum _{lm}4\pi {(i)}^l{g}_K^{STl}(r,t)Y_{lm}(\widehat{\mathit{r}})Y_{lm}^{\ast }(\widehat{\mathit{t}}).\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill \mathrm{with}:& & g_K^{STl}(r,t)\equiv f_K^{STl}(r,t)j_l(tr).\hfill \end{array}$$

(37)

The resulting integral equation in the channel $STl$ reads:

$$\begin{array}{ccc}& & g_K^{STl}(r,t)=j_l(tr)+{\int }_0^{\infty }4\pi y^{2}dy{\left[\frac{Q}{e}\right]}_K(r,y)V(y)g_K^{STl}(y,t)\hfill \end{array}$$

(38)

$$\begin{array}{ccc}& & \mathrm{with}:{\left[\frac{Q}{e}\right]}_K(r,y)=\int \frac{4\pi q^{2}dq}{{(2\pi )}^3}{j}_l(qr){\overline{\left[\frac{Q}{e}\right]}}_K(q)j_l(qy),\hfill \end{array}$$

(39)

$$\begin{array}{ccc}& & \mathrm{with}:{\overline{\left[\frac{Q}{e}\right]}}_K(q)=\int \frac{d\widehat{\mathit{q}}}{(4\pi )}\frac{P_{\mathbf{K}/\mathbf{2}+\mathbf{q}}{P}_{\mathbf{K}/\mathbf{2}-\mathbf{q}}}{{ϵ}_{\mathbf{K}/\mathbf{2}+\mathbf{t}}+{ϵ}_{\mathbf{K}/\mathbf{2}-\mathbf{t}}-{ϵ}_{\mathbf{K}/\mathbf{2}+\mathbf{q}}-{ϵ}_{\mathbf{K}/\mathbf{2}-\mathbf{q}}}.\hfill \end{array}$$

(40)

This equation is approximate in the sense that the Pauli-blocking operator divided by the energy denominator in relative momentum space has been replaced by its angle average. It nevertheless becomes exact for zero CM momentum. The functions $f_K^{STl}(r,t)=g_K^{STl}(r,t)/j_l(tr)$ play the role of correlation functions to be compared with the state-independent correlation function $f(r)$ of the variational LOCV approaches (see [34] and references therein). At variance with the variational framework, these Brueckner correlation functions depend on the two-body channel and on the initial state with total CM momentum K (P in [34]) and relative momentum t (q in [34]). However, according to the results of ref. [34], for the dominant $l=0$ Fourier component at normal nuclear density, one finds that the correlation functions for the ${}^1{S}_0(S=0,T=1)$ channel and the ${}^3{S}_1(S=1,T=0)$ channel (which is sensitive to the tensor force) calculated at given t and K values are very similar. In addition, it is found that these correlation functions are also very similar to the state-channel-independent correlation function of the variational LOCV method. Moreover, using the $v18$ Argonne potential [35], this calculation demonstrates that the dependence on the CM momentum and on the relative momentum is also weak. Hence, for convenience and practical reasons, but not only these reasons (see below), we propose representing the effect of short-range correlations by a unique state and channel-independent Jastrow function as:

$$f_{\alpha }(r)=1-\alpha j_0(q_{c}r).$$

(41)

Indeed, it was demonstrated in ref. [36] that, starting from the Bonn potential, this simple Jastrow ansatz can rather accurately reproduce the Brueckner reaction matrix in the dominant channels (but with state-dependent $q_c$ and $\alpha$ parameters). Of course, in the Brueckner theory, the correlation function has no reason to vanish (i.e, $\alpha =1$) at the origin. For instance, one systematically has $f(r=0)\sim 0.2$ (or $\alpha =0.8$ for the corresponding Jastrow ansatz) for various values of the relative momentum, of the CM momentum and of the density in the results reported in [34]. The fact that two nucleons may coexist at the same point is, mathematically, a consequence of modeling through the presence of the form factor transforming the hard core potential into a soft core potential; i.e., $V(r=0)$ remains finite. One can very well consider that this modeling approach is not appropriate at very short distances if we state that two composite nucleons cannot coexist at the same place as a six-quark cluster, which has never been observed. For this reason, one can take the following ansatz for the Jastrow function:

$$f_K^{STl}(r,t)=f_c(r)=1-j_0(q_{c}r)$$

(42)

satisfying $f_c(r=0)=0$. An immediate consequence is that the effective interaction automatically vanishes at high momentum, providing the natural regularization of the loop integral entering, for instance, the two-pion (or two-rho) exchange contribution to the energy per nucleon of nuclear matter. As discussed in [16], this property can be seen as a consequence of the Beg–Agassi theorem, which states that virtual mesons propagate freely inside the correlation hole. However, the function $g_K^{01l}(r,t)=f_c(r)j_0(tr)$ cannot be an exact solution of the integral Equation (38). To find the optimal choice, we proceed as follows: we enter the function $f_c(r)=1-j_0(q_{c}r)$ into the r.h.s of Equation (38) and calculate the output function ${\tilde{f}}_c(r)$. We thus fix the value of $q_c$ to obtain ${\tilde{f}}_c(r=0)=0$. Using typical parameters of the Bonn potential (potential A in [37]), we obtain a value of the order of $q_c\sim 600$ MeV, which changes very little if we vary the input parameters. We show the result of the calculation for a Bonn-like non-relativistic OBE potential but with the parameters taken from or inspired by the chiral confining model in its QCD-connected version, as described below in Section 3:

$$\begin{array}{ccc}& & g_{\sigma }=8.37,m_{\sigma }=919\mathrm{MeV},{\mathrm{\Lambda }}_S=1\mathrm{GeV},\hfill \\ & & g_V=7.5,m_V=783\mathrm{MeV},{\mathrm{\Lambda }}_V=1\mathrm{GeV},\hfill \\ & & g_A=1.26,F_{\pi }=92.4\mathrm{MeV},m_{\pi }=140\mathrm{MeV},{\mathrm{\Lambda }}_{\pi }=1\mathrm{GeV},\hfill \\ & & g_{\rho }=g_V/3,{\kappa }_{\rho }=6,C_{\rho }=1.865,m_{\rho }=770\mathrm{MeV},{\mathrm{\Lambda }}_{\rho }=2\mathrm{GeV}.\hfill \end{array}$$

(43)

In this version of the model, one actually introduces two scalar potentials: one associated with the QCD-connected scalar field, s, describing the fluctuations of the chiral condensate (see next section), and the other, ${\sigma }^{\prime }$, simulating the two-pion (and two-rho) exchange potential with the delta resonance in the intermediate state. In the traditional OBE model, the correlated and uncorrelated two-pion exchanges correspond to all the diagrams in Figure 3 in ref. [38] (apart from the iterative process with NN intermediate states). The parameters are:

$$g_{2\pi }=4.8\equiv g_{{\sigma }^{\prime }},m_{2\pi }=550\mathrm{MeV}\equiv m_{{\sigma }^{\prime }},{\mathrm{\Lambda }}_{2\pi }=1\mathrm{GeV}\equiv {\mathrm{\Lambda }}_{{\sigma }^{\prime }}.$$

(44)

We show the various contributions to this simple non-relativistic $NN$ potential in Figure 2 and note the existence of an attractive pocket in the full potential, which is expected to be slightly deepened when the tensor force (needed to produce the deuteron bound state) and the time component of the rho exchange are incorporated. We recall in passing that such a strong rho scenario (i.e., a strong (Lorentz) tensor coupling, $K_{\rho }=6$) is required to decrease the (Wigner) tensor force without obtaining too large a D-state probability in the deuteron.

In Figure 3, we show the output correlation function $\tilde{f}(r)$ after one G-matrix iteration for a two-nucleon system in the $l=0$ state at normal nuclear density for either the singlet channel, ${}^1{S}_0(S=0,T=1)$, or the triplet channel, ${}^3{S}_1(S=1,T=0)$. We fix $q_c=670$ MeV to have a vanishing correlation function for a relative momentum $t=100$ MeV. As can be seen in this figure, varying t induces a very moderate change in the correlation function, which is reminiscent of the very weak relative momentum dependence of the correlation function observed in ref. [34]. In addition, varying $g_{{\sigma }^{\prime }}$ and $m_{{\sigma }^{\prime }}$ by ∼20% or varying the density between zero and twice the nuclear matter density induces a change in $q_c$ of only a few tens of MeV.

If we now accept that the correlation function is well represented by this Jastrow function, the effective tensor force can be obtained very simply as $G_T(r)=V_T(r)f_c(r)$. For the total spin–isospin interaction written in momentum space, this explicitly gives:

$$\begin{array}{ccc}\hfill G_{\sigma \tau }(\mathit{q})& =& \int d\mathit{r}{e}^{-i\mathbf{q}·\mathbf{r}}[\int \frac{d\mathit{k}}{{(2\pi )}^3}{e}^{i\mathbf{k}·\mathbf{r}}(\frac{1}{3}\left(V_{\pi }(k)+2V_{\rho }(k)\right){\mathit{\sigma }}_1·{\mathit{\sigma }}_2\hfill \\ & & +\frac{1}{3}\left(V_{\pi }(k)-V_{\rho }(k)\right)\left(3{\mathit{\sigma }}_1·\widehat{\mathbf{k}}{\mathit{\sigma }}_2·\widehat{\mathbf{k}}-{\mathit{\sigma }}_1·{\mathit{\sigma }}_2\right)\left)\right]\left(1-j_0(q_{c}r)\right).\hfill \end{array}$$

(45)

One can adopt two different decompositions for the effective spin–isospin interaction: central and tensor or longitudinal and transverse,

$$\begin{array}{ccc}\hfill G_{\sigma \tau }(\mathit{q})& \hfill =& {\left(\frac{g_A}{2F_{\pi }}\right)}^2\left(v_c(q){\mathit{\sigma }}_1·{\mathit{\sigma }}_2+v_t(q)\left({\mathbf{3}\mathit{\sigma }}_1·\widehat{\mathbf{q}}{\mathit{\sigma }}_2·\widehat{\mathbf{q}}-{\mathit{\sigma }}_2·{\mathit{\sigma }}_2\right)\right)\hfill \end{array}$$

(46)

$$\begin{array}{ccc}& \hfill =& {\left(\frac{g_A}{2F_{\pi }}\right)}^2\left(v_L(q){\mathit{\sigma }}_1·\widehat{\mathbf{q}}{\mathit{\sigma }}_2·\widehat{\mathbf{q}}+v_T(q){\mathit{\sigma }}_1\times \widehat{\mathbf{q}}·{\mathit{\sigma }}_2\times \widehat{\mathbf{q}}\right),\hfill \end{array}$$

(47)

with

$$\begin{array}{ccc}\hfill v_c(q)& =& \frac{1}{3}\left(v_{\pi }(q)+2v_{\rho }(q)+3g^{\prime }(q)\right),v_t(q)=\frac{1}{3}\left(v_{\pi }(q)-v_{\rho }(q)+3h^{\prime }(q)\right)\hfill \\ \hfill v_L(q)& =& v_{\pi }(q)+g^{\prime }(q)+2h^{\prime }(q),v_T(q)=v_{\rho }(q)+g^{\prime }(q)-h^{\prime }(q),\hfill \end{array}$$

(48)

and

$$\begin{array}{ccc}\hfill v_{\pi }(q)& =& -{\mathrm{\Gamma }}_{\pi }^2(q^2)\frac{q^2}{q^2+m_{\pi }^2}\hfill \\ \hfill v_{\rho }(q)& =& -C_{\rho }{\mathrm{\Gamma }}_{\rho }^2(q^2)\frac{q^2}{q^2+m_{\rho }^2}\hfill \\ \hfill g^{\prime }(q)& =& -\frac{1}{3}\int \frac{d\mathit{k}}{{(2\pi )}^3}\left(v_{\pi }(k)+2v_{\rho }(k)\right)\frac{2{\pi }^2}{q_c}\delta \left(|\mathit{q}-\mathit{k}|-q_c\right)\hfill \\ & \simeq & \frac{1}{3}{\mathrm{\Gamma }}_{\pi }^2(q^2+q_c^2)\frac{q^2+q_c^2}{q^2+q_c^2+m_{\pi }^2}+\frac{2C_{\rho }}{3}{\mathrm{\Gamma }}_{\rho }^2(q^2+q_c^2)\frac{q^2+q_c^2}{q^2+q_c^2+m_{\rho }^2}\hfill \\ \hfill h^{\prime }(q)& =& -\frac{1}{3}\int \frac{d\mathit{k}}{{(2\pi )}^3}\left(v_{\pi }(k)-v_{\rho }(k)\right)\left(3{(\mathit{k}·\widehat{\mathit{q}})}^2-1\right)\frac{2{\pi }^2}{q_c}\delta \left(|\mathbf{q}-\mathbf{k}|-q_c\right)\hfill \\ & \simeq & \frac{1}{3}{\mathrm{\Gamma }}_{\pi }^2(q^2+q_c^2)\frac{q^2}{q^2+q_c^2+m_{\pi }^2}-\frac{C_{\rho }}{3}{\mathrm{\Gamma }}_{\rho }^2(q^2)\frac{q^2}{q^2+q_c^2+m_{\rho }^2}.\hfill \end{array}$$

(49)

The analytical forms of $g^{\prime }(q)$ and $h^{\prime }(q)$, which are exact for vanishing and large momenta, provide an accurate interpolation between these two limits. This is the well-known result quoted in [39]. In the Landau limit, we obtain the Landau–Migdal parameter $g^{\prime }=g^{\prime }(0)=0.59$, which is compatible with the range of accepted values [40]. We can also remark that a large value of the cutoff for the $\rho NN$ (Lorentz) tensor coupling is needed to obtain a value for g’ that is not too small. Finally, we also note that the whole G-matrix interaction is only made of terms in the form $F(q^2)-F(q^2+q_c^2)$, which vanishes in the UV for $q>>q_C$ as a consequence of the Beg–Agassi–Gal theorem.

It is well known that nuclear matter calculations based on a pure two-body interaction fail to simultaneously reproduce the binding energy and the saturation density (the famous Coester band problem) even if relativistic effects (Walecka mechanism) [41,42] included in the Dirac–Brueckner–Hartree–Fock (DBHF) approach may improve the situation, hence providing a guide for parameterization for an in-medium RMF with density-dependent coupling constants simulating many-body forces [43]. In the following, we follow another route, where effects linked to the quark substructure and the broken chiral symmetry of the QCD vacuum generate the needed three-body forces. This is performed in the chiral confining model that we describe in the next section.

3. The Chiral Confining Model

3.1. The Phenomenological Model

The very early motivation of this chiral model for dense nuclear matter (and neutron star interior) was to find a firm theoretical basis of the very successful relativistic theories initiated by Walecka and collaborators [41,42]. This type of approach indeed provides a very economical saturation mechanism, and a spectacular well-known success is the correct magnitude of the spin–orbit potential since the nucleons move in an attractive background scalar field and in a repulsive vector background field, which contribute in an additive way. Although the origin of the repulsive vector field can be safely identified as associated with the omega vector-meson exchange, the real nature of the attractive Lorentz scalar field has been a controversial subject since there is no sharp scalar resonance that would lead to a simple scalar particle exchange. More fundamentally, the question of the very nature of these background fields has to be elucidated; in other words, it is highly desirable to clarify their relationship with the QCD condensates and, in particular, with the chiral quark condensate and, more generally, with the low-energy realization of chiral symmetry, which is spontaneously broken in the QCD vacuum and is expected to be progressively restored when the density increases.

To bridge the gap between relativistic theories of the Walecka type and approaches insisting on chiral symmetry, it was proposed in [17] (see also the detailed discussion given in ref. [19]) that the “nuclear physics scalar sigma meson” of the relativistic Walecka model at the origin of the nuclear binding be identified; let us call it ${\sigma }_W$, with the chiral invariant $s=S-F_{\pi }$ field associated with the radial fluctuation of the chiral condensate S around the “chiral radius” $F_{\pi }$, identified to the pion decay constant. Said differently, we take the point of view that the effective theory has to be formulated, as a starting point, in terms of the fields associated with the fluctuations of the chiral quark condensate parameterized in a matrix form:

$$\begin{array}{ccc}\hfill W& =& \sigma +i\overrightarrow{\tau }·\overrightarrow{\pi }\equiv SU\equiv (s+F_{\pi })U\equiv ({\sigma }_W+F_{\pi })U\hfill \\ & & \mathrm{with}U(x)=e^{i\overrightarrow{\tau }·\overrightarrow{\varphi }(x)/F_{\pi }}.\hfill \end{array}$$

(50)

The scalar field $\sigma$ (S) and pseudoscalar fields $\overrightarrow{\pi }$ ($\overrightarrow{\varphi }$) written in Cartesian (polar) coordinates appear as dynamical degrees of freedom and may deviate from the vacuum value, ${〈\sigma 〉}_{vac}={〈S〉}_{vac}=f_{\pi }\propto {〈\overline{q}q〉}_{vac}$. The sigma and pion fields, associated with the amplitude $s\equiv {\sigma }_W$ and phase fluctuations $\varphi$ of this condensate, are promoted to the rank of effective degrees of freedom. Their dynamics are governed by an effective potential, $V_{\chi }\left(\sigma ,\overrightarrow{\pi })\right)$, having a typical Mexican hat shape associated with a broken (chiral) symmetry of the QCD vacuum.

There is, however, a well-identified problem concerning nuclear saturation with the usual chiral effective theories [44,45,46,47]. Independently of the particular chiral model, in the nuclear medium, one can move away from the minimum of the vacuum effective potential (Mexican hat potential), i.e., into a region of smaller curvature. This single-effect equivalent to the lowering of the sigma mass destroys the stability, creating problems for the applicability of such effective theories in the nuclear context. The effect can be associated with an $s^3$ tadpole diagram, generating attractive three-body forces and destroying saturation, even if the repulsive three-body force from the Walecka mechanism is present. One possible way to solve this problem is to introduce the nucleonic response to the scalar field, ${\kappa }_{\mathrm{NS}}$, which is the central ingredient of the quark–meson coupling model (QMC), introduced in the original pioneering work of P. Guichon [48] and successfully applied to finite nuclei with an explicit connection to the Skyrme force [49,50,51]. This effect, which is associated with the polarization of the quark substructure in the presence of the nuclear scalar field, will unavoidably generate three-body forces, which may provide the desired repulsion. In practice, this response or, more precisely, the nucleon scalar susceptibility ${\kappa }_{\mathrm{NS}}$ generates a non-linear coupling of the scalar field to the nucleon or, equivalently, a decrease in the scalar coupling constant with increasing density. Hence, to achieve saturation, in a set of successive works [18,20,21,22,25], we have complemented the relativistic chiral approach in such a way that the effect of the nucleon response is able to counterbalance the attractive chiral tadpole diagram to obtain good saturation properties, especially the correct curvature coefficient—the empirical incompressibility parameter.

Phenomenologically, the model is described with standard notations by the Lagrangian

$$\mathcal{L}=\overline{\mathrm{\Psi }}i{\gamma }^{\mu }{\partial }_{\mu }\mathrm{\Psi }+{\mathcal{L}}_s+{\mathcal{L}}_{\omega }+{\mathcal{L}}_{\rho }+{\mathcal{L}}_{\pi },$$

(51)

with

$$\begin{array}{ccc}\hfill {\mathcal{L}}_s& =& -M_N^{\ast }(s)\overline{\mathrm{\Psi }}\mathrm{\Psi }-V(s)+\frac{1}{2}{\partial }^{\mu }s{\partial }_{\mu }s\hfill \\ \hfill {\mathcal{L}}_{\omega }& =& -g_V{\omega }_{\mu }\overline{\mathrm{\Psi }}{\gamma }^{\mu }\mathrm{\Psi }+\frac{1}{2}{m}_V^2{\omega }^{\mu }{\omega }_{\mu }-\frac{1}{4}{F}^{\mu \nu }{F}_{\mu \nu }\hfill \\ \hfill {\mathcal{L}}_{\rho }& =& -g_{\rho }{\rho }_{a\mu }\overline{\mathrm{\Psi }}{\gamma }^{\mu }{\tau }_a\mathrm{\Psi }-g_{\rho }\frac{{\kappa }_{\rho }}{2M_N}{\partial }_{\nu }{\rho }_{a\mu }\mathrm{\Psi }{\overline{\sigma }}_{}^{\mu \nu }{\tau }_a\mathrm{\Psi }+\frac{1}{2}{m}_{\rho }^2{\rho }_{a\mu }{\rho }_a^{\mu }-\frac{1}{4}{G}_a^{\mu \nu }{G}_{a\mu \nu }\hfill \\ \hfill {\mathcal{L}}_{\pi }& =& \frac{g_A}{2F_{\pi }}{\partial }_{\mu }{\phi }_{a\pi }\overline{\mathrm{\Psi }}{\gamma }^{\mu }{\gamma }^5{\tau }_a\mathrm{\Psi }-\frac{1}{2}{m}_{\pi }^2{\phi }_{a\pi }^2+\frac{1}{2}{\partial }^{\mu }{\phi }_{a\pi }{\partial }_{\mu }{\phi }_{a\pi }.\hfill \end{array}$$

(52)

It involves the scalar field s (the “nuclear physics sigma meson” ${\sigma }_W$), the pion field ${\phi }_{a\pi }$, and the vector fields associated with the omega meson channel (${\omega }^{\mu }$) and with the rho meson channel (${\rho }_a^{\mu }$). Each meson-nucleon vertex is regularized by a meson-nucleon form factor mainly originating from the compositeness of the nucleon, hence generating a bare OBE Bonn-like NN interaction, as given in Equation (28). The specific crucial ingredients beyond the simplest approach are the presence of a chiral effective potential associated with the chirally broken vacuum and an in-medium-modified nucleon mass, which is supposed to embed its quark substructure.

In the majority of our previous works [18,19,20,21,22,24,25], the chiral effective potential had the simplest linear sigma model (L$\sigma$M) form:

$$V_{\chi ,L\sigma M}(s)=\frac{1}{2}{M}_{\sigma }^2{s}^2+\frac{1}{2}\frac{M_{\sigma }^2-M_{\pi }^2}{F_{\pi }}{s}^3+\frac{1}{8}\frac{M_{\sigma }^2-M_{\pi }^2}{F_{\pi }^2}{s}^4.$$

(53)

The effective Dirac nucleon mass $M_N^{\ast }(s)$ deviates from the bare nucleon mass in the presence of the nuclear scalar field s:

$$\begin{array}{ccc}\hfill M_N^{\ast }(s)& =& M_N+g_{S}s+\frac{1}{2}{\kappa }_{NS}{s}^2+\mathcal{O}(s^3).\hfill \end{array}$$

(54)

where $g_S$ is the scalar coupling constant of the model. In [18,19,20,21,22,24], we took the pure linear sigma model value $g_S=M_N/F_{\pi }$, but in the most recent work [25], $g_S$ was allowed to deviate from L$\sigma$M and was fixed by performing a Bayesian analysis. This quantity actually corresponds to the first-order response of the nucleon to an external scalar field and can be obtained in an underlying microscopic model of the nucleon. The nucleon scalar susceptibility ${\kappa }_{\mathrm{NS}}$ is another response parameter and reflects the polarization of the nucleon, i.e., the self-consistent readjustment of the quark wave function in the presence of the scalar field. Very generally, the scalar coupling constant, $g_S$, and the nucleon response parameter, ${\kappa }_{NS}$, depend on the subquark structure and the confinement mechanism, as well as the effect of spontaneous chiral symmetry breaking. In our previous works [18,20,21,22,24,25], we introduced a dimensionless parameter,

$$C\equiv \frac{{\kappa }_{\mathrm{NS}}{F}_{\pi }^2}{2M_N},$$

(55)

which is expected to be of the order $C\sim 0.5$ as in the MIT bag used in the QMC framework. The physical motivation to introduce this nucleonic response is the observation that nucleons experience huge fields at finite density, e.g., the scalar field is of the order of a few hundred MeV at saturation density. Nucleons, being composite objects in reality, will react to the nuclear environment (i.e., the background nuclear scalar field) through the (self-consistent) modification of quark wave functions. This effect may generate a three-body force that provides the desired repulsion if confinement dominates spontaneous chiral symmetry breaking in the nucleon mass origin, as discussed in ref. [23] within particular models. Indeed, it is possible to show that the scalar sector generates a three-nucleon repulsive contribution to the energy per nucleon:

$$\frac{E^{(3)}}{A}\simeq \frac{g_S^2}{2M_{\sigma }^4}\left({\kappa }_{\mathrm{NS}}-\frac{g_S}{F_{\pi }}\right){\rho }_s^2=\frac{g_S^3}{2M_{\sigma }^4{F}_{\pi }}\left(2{\tilde{C}}_3-1\right){\rho }_s^2\mathrm{with}{\tilde{C}}_3\simeq \frac{M_N}{g_S{F}_{\pi }}C.$$

(56)

This is the result previously quoted in Equation (44) in ref. [20], but without the factor $M_N/g_S{F}_{\pi }$, which appears when the scalar coupling constant is allowed to deviate from its pure L$\sigma$M value, as in ref. [25]. This three-body effect provides a very natural mechanism for the saturation mechanism but is at variance with the QMC model, which ignores the attractive tadpole diagram present in the chiral approach, this requires a C parameter that is close to or even larger than one [18,20,21,22,24,25].

One very important point is that it is possible to relate these scalar response parameters to the chiral properties of the nucleon, namely, the first derivative of the nucleon mass with respect to to the current quark mass (or to the pion mass), which is related to the light-quark sigma commutator and the second derivative, linked to the chiral susceptibility of the nucleon. The nucleon mass and other intrinsic properties of the nucleon (sigma term, chiral susceptibilities) are indeed QCD quantities, which are, in principle, obtainable from lattice simulations. The problem is that lattice calculations of this kind are difficult for small current quark mass, m, or equivalently, small pion mass, $M_{\pi }^2$. Here, $M_{\pi }$ represents the pion mass to the leading order in the quark mass (i.e., ignoring the NLO chiral logarithm correction), $M_{\pi }^2=2mB=-2m{\langle \overline{q}q\rangle }_{\chi L}/F^2$ (Gell–Mann–Oakes–Renner (GOR) relation). The quantities B and F, the pion decay constant in the chiral limit, are two low-energy parameters appearing in chiral perturbation theory [52]. The difficulty of the extrapolation is linked to the nonanalytical behavior of the nucleon mass as a function of m (or equivalently, $M_{\pi }^2$), which comes from the pion cloud contribution. The idea of the Adelaide group [53,54,55,56] was to separate the pion cloud self-energy, ${\Sigma }_{\pi }(M_{\pi }^2,\mathrm{\Lambda })$, from the rest of the nucleon mass and to calculate it with just one adjustable cutoff parameter $\mathrm{\Lambda }$ entering the $\pi NN,\pi N\mathrm{\Delta }$ form factor regularizing the pion loops. Actually, different cutoff forms for the pion loops (Gaussian, dipole, monopole, sharp) were used with the adjustable parameter $\mathrm{\Lambda }$. This formulation of ChiPT is thus called the Finite-Range Regulator (FFR) method. The remaining nonpionic part is expanded in terms of powers of $M_{\pi }^2$ as follows:

$$M_N(M_{\pi }^2)=a_0+a_2{M}_{\pi }^2+a_4{M}_{\pi }^4+\cdots +{\Sigma }_{\pi }(M_{\pi },\mathrm{\Lambda }),$$

(57)

Depending of the details of the chiral extrapolation, the extracted values of the parameters are within the range of $a_2\simeq 1.5{\mathrm{GeV}}^{-1},a_4\simeq -0.5{\mathrm{GeV}}^{-3}$ [54] and $a_2\simeq 1.0{\mathrm{GeV}}^{-1},a_4\simeq -0.25{\mathrm{GeV}}^{-3}$ [56]. The explicit connection between lattice QCD parameters $a_2$ and $a_4$ and the response parameters $g_S$ and C in the L$\sigma$M case was obtained in [20,21]:

$$a_2=\frac{F_{\pi }{g}_S}{M_{\sigma }^2},a_4=-\frac{F_{\pi }{g}_S}{2M_{\sigma }^4}\left(3-2{\tilde{C}}_3\right).$$

(58)

Notice that in the expression of $a_4$, the factor $M_N/F_{\pi }{g}_s$ in ${\tilde{C}}_3$, present in our recent paper [25], was absent in [20,21] since the nucleon mass was fixed to be $M_N=F_{\pi }{g}_s$. The quantity $a_2{M}_{\pi }^2\sim$ 20–30 MeV represents the nonpionic piece of the sigma commutator directly associated with the scalar field, s (see the detailed discussion of this quantity in ref. [20]). One very robust conclusion is that the lattice result for $a_4$ is much smaller than that obtained in the simplest linear sigma model, ignoring the nucleonic response ($C=0$), for which $a_4\simeq 3.5{\mathrm{GeV}}^{-3}$. Hence, lattice data require a strong compensation from the effects governing the three-body repulsive force needed for the saturation mechanism: compare Equations (58) and (56). Hence, both lattice data constraints and nuclear matter phenomenology require quite a large value of the dimensionless response parameter, C, which must be at least larger than one. Moreover, in a recent work based on a Bayesian analysis with lattice data as an input [25], we found that the response parameter is strongly constrained to the value $C\sim 1.4$, very close to the value where the scalar susceptibility changes in sign: $C=1.5$.

3.2. The NJL Chiral Confining Model

The problem that one has to face is that it seems impossible to find a realistic confining model for the nucleon able to generate C larger than one. For instance, as mentioned above, in the MIT bag model used in the QMC scheme, one has $C_{\mathrm{MIT}}\simeq 0.5$. One possible reason for this discrepancy between models and phenomenological values of C lies in the use of the linear sigma model (L$\sigma$M), which is probably too naive. Hence, one should certainly use an enriched chiral effective potential from a model able to give a correct description of the low-energy realization of chiral symmetry in the hadronic world. A good, easily tractable candidate is the Nambu–Jona–Lasinio (NJL) model. Indeed, in ref. [23], an explicit construction of the background scalar field was performed in the NJL model using a bosonization technique based on an improved derivative expansion valid at low (space-like) momenta [57]. Various confining interactions have been incorporated (quark–diquark string interaction, linear and quadratic confining interaction) on top of the NJL model, which seem to be sufficient to generate saturation, although the response parameter C remains relatively small, of the order of $C\sim 0.5$. The reason is that, for a given scalar mass, the NJL chiral effective potential generates a significantly smaller $s^3$ cubic term (hence generating a smaller attractive tadpole diagram) than the simplistic linear sigma model. We discuss this point below in more detail and demonstrate that the repulsive three-body force generating saturation is determined not only by the nucleon response C but also by the cubic term of the NJL potential, hereafter governed by the parameter $C_{\chi }$, with the parameters C and $C_{\chi }$ combining together in the three-body interaction. Below, we briefly outline this approach, which is described in detail in a parallel paper [26].

Let us now consider the NJL model defined by the Lagrangian:

$$\begin{array}{ccc}\hfill \mathcal{L}& =& \overline{\psi }\left(i{\gamma }^{\mu }{\partial }_{\mu }-m\right)\psi +\frac{G_1}{2}\left[{\left(\overline{\psi }\psi \right)}^2+{\left(\overline{\psi }i{\gamma }_5\overrightarrow{\tau }\psi \right)}^2\right]\hfill \\ & & -\frac{G_2}{2}\left[{\left(\overline{\psi }{\gamma }^{\mu }\overrightarrow{\tau }\psi \right)}^2+{\left(\overline{\psi }{\gamma }^{\mu }{\gamma }_5\overrightarrow{\tau }\psi \right)}^2+{\left(\overline{\psi }{\gamma }^{\mu }\psi \right)}^2\right].\hfill \end{array}$$

(59)

It depends on four parameters: the coupling constants $G_1$ (scalar) and $G_2$ (vector), the current quark mass m and a (noncovariant) cutoff parameter $\mathrm{\Lambda }$. Three of these parameters ($G_1$, m and $\mathrm{\Lambda }$) are adjusted to reproduce the pion mass, the pion decay constant and the quark condensate. We refer the reader to [23] for more details. Using path integral techniques (i.e., integrating out $q\overline{q}$ pairs in the Dirac sea) and after the chiral rotation of the quark field, it can be equivalently written in a semi-bosonized form involving a pion field $\overrightarrow{\varphi }$ embedded in the unitary operator $U={\xi }^2=exp(i\overrightarrow{\tau }·\overrightarrow{\varphi }(x)/F_{\pi })$, a scalar field, $\mathcal{S}$, a vector field, $V^{\mu }$, and an axial-vector field, $A^{\mu }$, with all these meson fields being coupled to the constituent valence quarks. It has the explicit form given in Equations (2) and (7)–(11) in ref. [23]. Subtracting vacuum expectation values, the chiral effective potential can be expressed as:

$$V_{\chi ,\mathrm{NJL}}(s)=-2N_c{N}_f\left(I_0(\mathcal{S})-I_0(M_0)\right)+\frac{{\left(\mathcal{S}-m\right)}^2-{\left(M_0-m\right)}^2}{2G_1}.$$

(60)

The quantity $-2N_c{N}_f{I}_0(\mathcal{S})$ is nothing but the total (in-medium) energy of the Dirac sea of constituent quarks with the NJL loop integral $I_0(\mathcal{S})$ given in ref. [26]. The $\mathcal{S}$ fields whose vacuum expectation value coincides with the vacuum constituent quark mass $M_0$ is related to the “nuclear physics sigma meson” field s according to:

$$\mathcal{S}\equiv \frac{M_0}{F_{\pi }}\left(s+F_{\pi }\right).$$

(61)

For a comparison with the usual RMF model using the L$\sigma$M chiral effective potentials in Equation (53) or, equivalently, non-linear sigma couplings, we expand the effective potential to the third order in s as:

$$V_{\chi ,\mathrm{NJL}}(s)=V_{\chi }(0)+V_{\chi }^{\prime }(0)s+\frac{1}{2}{V}_{\chi }^{\prime \prime }(0)s^2+\frac{1}{6}{V}_{\chi }^{\prime \prime \prime }(0)s^3+\dots$$

(62)

An explicit calculation of the derivatives of the potential yields [26]:

$$V_{\chi ,\mathrm{NJL}}(s)=\frac{1}{2}{M}_{\sigma }^2{s}^2+\frac{1}{2}\frac{M_{\sigma }^2-M_{\pi }^2}{F_{\pi }}{s}^3\left(1-C_{\chi ,\mathrm{NJL}}\right)+\dots ,$$

(63)

The effective sigma mass $M_{\sigma }\sim 2M_0$ and the pion mass $M_{\pi }$ and $F_{\pi }$ are calculated within the model. The quantity $C_{\chi ,\mathrm{NJL}}$ is another NJL parameter, whose expression in terms of the NJL loop integral is given in [26]. For a large cutoff value, the $C_{\chi ,\mathrm{NJL}}$ parameter goes to zero, but for typical values of the NJL parameters, its value is in the range 0.4–0.5. Figure 4 shows that the approximate expansion (63) reproduces the exact NJL potential (60) very well. In this figure, the values of the NJL parameters are taken from the QCD-connected model presented in the next subsection. Comparing L$\sigma$M with the NJL scalar potential, one sees that the attractive tadpole term is larger in the case of L$\sigma$M for the same value of the effective sigma mass, $M_{\sigma }$. The effect of the parameter $C_{\chi ,\mathrm{NJL}}$ is then the reduction in the attractive tadpole diagram and the greater repulsion of the chiral potential. In the same figure, we also show the potential of the original Walecka model (again for the same effective sigma mass) that is limited to the quadratic term.

As derived in [26], to the leading order in density, the scalar sector generates a three-nucleon contribution contribution to the energy per nucleon, which is modified according to (compare with Equation (56)):

$$\begin{array}{ccc}& & \frac{E^{(3b)}}{A}\simeq \frac{g_S^2}{2M_{\sigma }^4}\left({\kappa }_{\mathrm{NS}}-\frac{g_S{V}_{\chi }^{\prime \prime \prime }(0)}{3M_{\sigma }^2}\right){\rho }_s^2=\frac{g_S^3}{2M_{\sigma }^4{F}_{\pi }}\left(2{\tilde{C}}_3-1\right){\rho }_s^2\hfill \\ & & \mathrm{with}{\tilde{C}}_3\simeq \frac{M_N}{g_S{F}_{\pi }}C+\frac{1}{2}{C}_{\chi }.\hfill \end{array}$$

(64)

One very important point is that, as established in [26], the relationship between lattice QCD parameters is also modified according to (compare with Equation (58)):

$$a_2=\frac{F_{\pi }{g}_S}{M_{\sigma }^2},a_4=-\frac{F_{\pi }{g}_S}{2M_{\sigma }^4}\left(3-2{\tilde{C}}_L\right)\mathrm{with}{\tilde{C}}_L=\frac{M_N}{g_S{F}_{\pi }}C+\frac{3}{2}{C}_{\chi }.$$

(65)

The important conclusion discussed in detail in ref. [26] is that for a given three-nucleon force allowing saturation, one needs a lower value of the dimensionless parameter C to obtain a small value of the lattice parameter $a_4$ compatible with lattice data. To provide insight into the effect of an enriched chiral effective potential, we return to our our original paper [18]. In this paper, where L$\sigma$M was used ($g_S=M_N/F_{\pi }$), we obtained the correct saturation properties with $C=1$ (see Figure 4 in this paper). We can retrospectively calculate the $a_2$ and $a_4$ parameters: we find $a_2=1.67$ GeV${}^{-1}$ and $a_4=-1.48$ GeV${}^{-3}$. If the value obtained for $a_2$ is not very far from the lattice values, $a_4$ is at a magnitude three times larger than the upper value compatible with the lattice calculation. We can now play a little game by just incorporating the NJL-like $(1-C_{\chi })$ correction in the cubic term term of the L$\sigma$M chiral effective potential with $C_{\chi }=0.44$. Keeping all the other parameters at their original value, we take $C=0.78$ so as to keep the same value of the repulsive three-body force, i.e., ${\tilde{C}}_3=C+C_{\chi }/2=1$ (64). The new saturation point remains very close to the original one, but the $a_4$ parameter becomes very close to zero, $a_4=-0.1$ GeV${}^{-3}$, in much better agreement with the lattice data mentioned above, namely, $a_4\simeq -0.5{\mathrm{GeV}}^{-3}$ [54] and $a_4\simeq -0.25{\mathrm{GeV}}^{-3}$ [56]. See also the discussion given in [26]. The conclusion is that the use of an enriched chiral effective potential of the NJL type significantly increases the agreement with lattice data, together with the expected model values of the nucleonic response parameter C. In the last section, we present the calculation of the nuclear matter equation of state with parameters taken from the QCD-connected model. We show explicitly that the inclusion of the $C_{\chi }$ parameter will generate a correct value of the $a_4$ parameter, together with a moderate value of the C parameter, while preserving the saturation properties.

3.3. The QCD-Connected Chiral Confining Model

The general picture underlying our approach can be summarized as follows. Nuclear matter is made of nucleons, which are themselves built from quarks and gluons and look like Y-shaped strings generated by a nonperturbative confining force, with constituent quarks at the ends. These quarks acquire a large mass from a quark condensate, which is the order parameter associated with the spontaneous breaking of chiral symmetry in the QCD vacuum. When the density $\rho$ of nuclear matter increases, the QCD vacuum is modified by the presence of the nucleons, yielding a decrease in the quark condensate, indicating the progressive restoration of chiral symmetry. Hence, what is usually called “the nuclear medium” can be seen as a “shifted vacuum” with a lower value of the order parameter. The mass of the constituent quarks coincides with the in-medium expectation value, $M=\overline{\mathcal{S}}(\rho )$, of the chiral-invariant scalar field $\mathcal{S}$, associated with the radial fluctuation mode of the chiral condensate. To implement such a physical picture, in this section, we propose an effective model of low-energy QCD that incorporates both chiral symmetry breaking, i.e., the condensation of quark–antiquark pairs in the QCD vacuum, and a confining string force that binds the massive constituent quarks inside the nucleon. It is based on the field correlator method (FCM) developed by Y. Simonov and collaborators [58,59,60,61,62,63,64,65,66], where the starting point is the Euclidean QCD Lagrangian and partition function written with conventional notation:

$$\begin{array}{ccc}\hfill L_{\mathrm{QCD}}& =& L_0+L_1+L_G,Z=\int d\psi d\overline{\psi }d{A}_{\mu }{e}^{-\left[L_0+L_1+L_G\right]},\hfill \\ \hfill \mathrm{with}{L}_0& =& \int d^{4}x\overline{\psi }(x)\left({\gamma }_{\mu }{\partial }_{\mu }+m\right)\psi (x),\hfill \\ \hfill L_1& =& \int d^{4}x\overline{\psi }(x)ig{\gamma }_{\mu }{A}_{\mu }\psi (x),L_G=-\frac{1}{4}\int d^{4}x{F}_{\mu \nu }^a{F}_{\mu \nu }^a.\hfill \end{array}$$

(66)

In the first step, we integrate over the colored gluon fields ($A_{\mu }\equiv A_{\mu }^a{t}_a$) to generate a pure quark effective action. This is the gluon field averaging briefly described below. However, before carrying out this gluon integration, one makes a gauge choice (modified Fock–Schwinger gauge [58,59,60,61,62,63,64,65]) such that

$$A_{\mu }\left(x_4,{\mathbf{r}}_0\right)=0,{\left(\mathbf{x}-{\mathbf{r}}_0\right)}_j·A_j\left(x_4,\mathbf{x}\right)=0,$$

(67)

where ${\mathbf{r}}_0$ is an arbitrary point, which will be identified later with a string junction coordinate. In this particular gauge, it is possible to express $A_{\mu }$ in terms of the field strength tensor $F_{\mu \nu }$ ([58,59,60,61,62,63,64,65]):

$$\begin{array}{ccc}\hfill A_4\left(x_4,\mathbf{x}\right)& =& {\int }_0^{1}dv{\left(\mathbf{x}-{\mathbf{r}}_0\right)}_k{F}_{k4}\left(x_4,{\mathbf{r}}_0+v\left(\mathbf{x}-{\mathbf{r}}_0\right)\right)\hfill \\ \hfill A_j\left(x_4,\mathbf{x}\right)& =& {\int }_0^{1}dvv{\left(\mathbf{x}-{\mathbf{r}}_0\right)}_k{F}_{kj}\left(x_4,{\mathbf{x}}_0+v\left(\mathbf{x}-{\mathbf{r}}_0\right)\right).\hfill \end{array}$$

(68)

Hence, the vector field appears as a contour integral, with the contour $C(\mathbf{x})$ reducing to a straight line $\mathbf{z}(v)={\mathbf{r}}_0+v\left(\mathbf{x}-{\mathbf{r}}_0\right)$ between the $x_4$ axis and the point where the gluon field is calculated. The partition function can be written in successive forms,

$$\begin{array}{ccc}\hfill Z& =& \int d\psi d\overline{\psi }d{A}_{\mu }{e}^{-\left[L_0+L_1+L_G\right]}\equiv \int d\psi d\overline{\psi }{e}^{-L_0}{〈e^{-L_1}〉}_G\hfill \\ & \equiv & \int d\psi d\overline{\psi }{e}^{-\left[L_0+L_1^{\mathrm{eff}}\right]}\equiv \int d\psi d\overline{\psi }{e}^{-L^{\mathrm{eff}}},\hfill \end{array}$$

(69)

where, in the second form, gluon field averaging has been performed, which subsequently defines the effective quark Lagrangian $L^{\mathrm{eff}}$ once the gluon field has been integrated out. In the field correlator method, the latter quantity can be obtained as an expansion of the exponential ${〈e^{-L_1}〉}_G$. Using cluster expansion [58,59,61,62], $L^{\mathrm{eff}}$ can be written as an infinite sum containing averages such as $\langle {(A_{\mu })}^k\rangle$. According to ref. [62], one can make a Gaussian approximation, neglecting all correlators $\langle {(A_{\mu })}^k\rangle$ of degrees higher than $k=2$. The numerical accuracy has been studied on the lattice, demonstrating that the corrections to the Gaussian approximation are not larger than $3\%$ [62,67]. Hence, neglecting all higher correlators, this effective Lagrangian (or more precisely, this effective action) reads

$$\begin{array}{ccc}\hfill L^{\mathrm{eff}}& \simeq & L_0+\frac{1}{2}{〈L_1{L}_1〉}_G=\int d^{4}x\overline{\psi }(x)\left({\gamma }_{\mu }{\partial }_{\mu }+m\right)\psi (x)\hfill \\ & & +\frac{1}{2}\int d^{4}x{d}^{4}y\overline{\psi }(x){\gamma }_{\mu }{t}_a\psi (x)\overline{\psi }(y){\gamma }_{\nu }{t}_b\psi (y)\frac{{\delta }_{ab}}{C_F}{J}^{\mu \nu }(x,y),\hfill \end{array}$$

(70)

with $C_F=(N_c^2-1)/2N_c=4/3$ and

$$\begin{array}{ccc}& & J^{\mu \nu }(x,y)=\frac{g^2}{2N_c}〈A_a^{\mu }(x)A_a^{\nu }(y)〉\hfill \\ & =& {\int }_0^{1}dv{\int }_0^{1}dw{\alpha }_{\mu }(v){\alpha }_{\nu }(w){\left(\mathbf{x}-{\mathbf{r}}_0\right)}_k{\left(\mathbf{y}-{\mathbf{r}}_0\right)}_q\frac{g^2}{2N_c}〈F_a^{k\mu }\left(z(v)\right)F_a^{q\nu }\left(z^{\prime }(w)\right)〉,\hfill \end{array}$$

(71)

with $z_4=x_4$, $\mathbf{z}(v)={\mathbf{r}}_0+v\left(\mathbf{x}-{\mathbf{r}}_0\right)$, $z_4^{\prime }=y_4$, ${\mathbf{z}}^{\prime }(w)={\mathbf{r}}_0+w\left(\mathbf{y}-{\mathbf{r}}_0\right)$, ${\alpha }_4(v)=1$, ${\alpha }_k(v)=v$. This expression of the kernel can be shown to be gauge-invariant, the reason being that the parallel transporters on the contour $C(x,y)$ are identically equal to unity in this gauge. The whole nonperturbative physics is contained in the nonlocal gluon condensate $g^2〈F_a^{k\mu }\left(z(v)\right)F_a^{q\nu }\left(z^{\prime }(w)\right)〉$, and parameterization is known from lattice measurements [66]. As in many previous works by Simonov and collaborators, we keep only the nonperturbative confining piece,

$$\frac{g^2}{2N_c}〈F_a^{\rho \mu }(z)F_a^{\lambda \nu }(z^{\prime })〉=\left({\delta }_{\rho \lambda }{\delta }_{\mu \nu }-{\delta }_{\rho \nu }{\delta }_{\mu \lambda }\right)D(z-z^{\prime }),$$

(72)

where $D(x)$ decreases in all directions and describes the profile of the bilocal correlator of the nonperturbative gluonic fields in the QCD vacuum. It depends on two QCD quantities, the gluon condensate ${\mathcal{G}}_2\sim 0.015$ GeV${}^4$ and the gluon correlation length $T_g$, which physically corresponds to the string width. It may also be parameterized in terms of two QCD parameters, the string tension $\sigma =0.18Ge{V}^2$ and $T_g$. We adopt a convenient Gaussian parameterization approach [61,62]:

$$D(x)=D(0)e^{-x^2/4T_g^2}\mathrm{with}D(0)=\frac{\sigma }{2\pi T_g^2}=\frac{{\pi }^2}{18}{\mathcal{G}}_2,\mathrm{hence}\sigma =\frac{{\pi }^3}{9}{\mathcal{G}}_2{T}_g^2.$$

(73)

This form of the bilocal correlator, which has the advantage of simplifying the calculations, was justified in [61,62]. In short, since all the observables are integrals of $D(x)$, its explicit form is not essential at large distances, provided that it has a finite range $T_g$ and that the string tension, i.e., the coefficient in the area law of the Wilson loop, is equal to $\sigma =\int d^{2}uD(u)/2$. Numerically, the gluon correlation length can be estimated as:

$$T_g=0.336{\left(\frac{\sigma (({\mathrm{GeV}}^2)}{0.2}\right)}^{1/2}{\left(\frac{0.02}{{\mathcal{G}}_2({\mathrm{GeV}}^4)}\right)}^{1/2}\mathrm{fm}.$$

(74)

A “small” parameter is:

$${\eta }^2=\sigma T_g^2=0.2\left(\frac{\sigma ({\mathrm{GeV}}^2)}{0.2}\right){\left(\frac{T_g(\mathrm{fm})}{0.1973}\right)}^2=0.462\left(\frac{\sigma ({\mathrm{GeV}}^2)}{0.2}\right){\left(\frac{T_g(\mathrm{fm})}{0.3}\right)}^2.$$

(75)

As in references [61,62], we will ignore the magnetic piece of the kernel, keeping only the dominant electric piece $J_{44}$. In addition, we make a static approximation, ignoring retardation effects:

$$e^{-x_4^2/4T_g^2}\simeq 2\sqrt{\pi }{T}_g\delta (x_4).$$

(76)

This approximation, already made in ref. [68] in the context of the heavy–light quark system, can be valid if the energy scale $T_g^{-1}\sim 700–800$ MeV is bigger than the other scale $\sigma \sim 400$ MeV of the problem or, equivalently, if the parameter $\eta$ is smaller than one. Hence, we end up with an effective static Lagrangian and, consequently, a more tractable effective Hamiltonian governing the dynamics of light quarks in the presence of a three-quark junction (or a static heavy quark) placed in ${\mathbf{r}}_0$ in the QCD vacuum. Said differently, this effective Hamiltonian has to be seen as a way to derive baryonic Green’s function and bound-state equations [69,70,71,72], already given in Refs. [58,59,62]. In principle, the point ${\mathbf{r}}_0$ is arbitrary and should be chosen as the one minimizing the string lengths (Torricelli point) joining the three quarks (see the discussion after Equation (28) in [62]). In practice, we take it as a constant parameter coinciding with the three-quark string junction, and the three-quark wave function is simply expressed in a factorized form in terms of single-quark orbitals centered in ${\mathbf{r}}_0$ [62]. Moreover, as discussed below, a strong point of this approach is that vacuum chiral symmetry breaking will show up as being intimately related to the confinement mechanism. Physically, this means that the three-quark string keeping the quarks together is constructed on top of a chirally broken vacuum, building a constituent quark at the end of the strings [68]. Indeed, the effective Hamiltonian reads:

$$H=\int d^{3}x{\psi }^{†}(\mathbf{x})\left(-i\overrightarrow{\alpha }·\overrightarrow{\nabla }+m\right)\psi (\mathbf{x})+\frac{1}{2}\int d^{3}x\int d^{3}y{\psi }^{†}(\mathbf{x})t_a\psi (\mathbf{x})V\left(\mathbf{x},\mathbf{y}\right){\psi }^{†}(\mathbf{y})t_a\psi (\mathbf{y}).$$

(77)

Introducing $\mathbf{R}=(\mathbf{x}+\mathbf{y})/2-{\mathbf{r}}_0$, $\mathbf{r}=\mathbf{x}-\mathbf{y}$, $\mathbf{X}=\mathbf{x}-{\mathbf{r}}_0$ and $\mathbf{Y}=\mathbf{y}-{\mathbf{r}}_0$, the potential V, derived from $J_{44}$ [61,62], can be written in various forms:

$$\begin{array}{ccc}\hfill V\left(\mathbf{x},\mathbf{y}\right)& =& \frac{1}{C_F}\frac{\sigma }{\sqrt{\pi }{T}_g}\left(\mathbf{x}-{\mathbf{r}}_0\right)·\left(\mathbf{y}-{\mathbf{r}}_0\right){\int }_0^{1}dv{\int }_0^{1}dw{e}^{-\frac{{\left(v\left(\mathbf{x}-{\mathbf{r}}_0\right)-w\left(\mathbf{y}-{\mathbf{r}}_0\right)\right)}^2}{4T_g^2}}\equiv \hfill \\ \hfill V\left(\mathbf{R},\mathbf{r}\right)& =& \frac{1}{C_F}\frac{\sigma }{\sqrt{\pi }{T}_g}\left({\mathbf{R}}^2-\frac{{\mathbf{r}}^2}{4}\right){\int }_0^{1}dv{\int }_0^{1}dw{e}^{-{\left(\frac{\left(v+w\right)\mathbf{r}}{4T_g}+\frac{\left(v-w\right)\mathbf{R}}{2T_g}\right)}^2}\equiv \hfill \\ \hfill V\left(\mathbf{X},\mathbf{Y};\mathbf{r}\right)& =& \frac{1}{C_F}\frac{\sigma }{\sqrt{\pi }{T}_g}\left(\frac{{\mathbf{X}}^2+{\mathbf{Y}}^2}{2}-\frac{{\mathbf{r}}^2}{2}\right)I\left(\mathbf{X},\mathbf{Y};\mathbf{r}\right)\hfill \\ \hfill \mathrm{with}& & I\left(\mathbf{X},\mathbf{Y};\mathbf{r}\right)={\int }_0^{1}dv{\int }_0^{1}dw{e}^{-{\left(\frac{\left(v+w\right)\mathbf{r}}{4T_g}+\frac{\left(v-w\right)}{2T_g}\frac{\mathbf{X}+\mathbf{Y}}{2}\right)}^2}.\hfill \end{array}$$

(78)

The relative variable $\mathbf{r}$ is the one associated with the self-interaction of the quark. Keeping only the $\mathbf{r}$ dependence in the second form of Equation (79), namely, taking $\mathbf{R}=0$, i.e, ${\mathbf{r}}_0$ at the mid-point between $\mathbf{x}$ and $\mathbf{y}$, the $V(\mathbf{r})$ interaction allows the generation of a BCS-like chirally broken vacuum and a mass gap equation for the constituent quark. The light $q\overline{q}$ mesons and, in particular, the Goldstone pion mode can be generated using this approach in a way equivalent to the RPA in the many-body nuclear problem [73]. It is also equivalent to the quark Dyson–Schwinger (gap) and bound-state Bethe–Salpeter equations [69,74,75]. In the FCM approach, this program is implemented using a Fierz transformation and a bosonization procedure, taking the origin ${\mathbf{r}}_0$ at the mid-point between the quark and the antiquark [63,64,65]; see also Section VII-D of ref. [76] for a short summary of this method. The variable $\mathbf{R}$ physically corresponds to the length of the confining string. Indeed, at a large value of R, one can analytically check that the integral of the Gaussian factor in Equation (79) has a $1/R$ behavior, which, together with the $R^2$ prefactor, yields a linear behavior: $C_{F}V(R,r)\sim \sigma R$. We see (Equation (79)) that, in practice, the self-interacting part and the string-like confining piece of the kernel might be mixed up in a very complicated way in the case of the presence of an “external” source, either a string junction or a static heavy (anti)quark. The treatment of chiral symmetry breaking versus confinement entanglement is certainly the most prominent problem of QCD, which already manifests at the level of the bare nucleon, as well as at the level of dense and hot nuclear or hadronic matter. However, some approximation schemes have been developed in the literature to partially disentangle these two aspects [68,77]. Inspired by the third writing of the potential in Equation (79), we propose a specific ansatz interpolating between the ansatz used in Refs. [26,68]:

$$V\left(\mathbf{X},\mathbf{Y};\mathbf{r}\right)={\tilde{V}}_C(\mathbf{X})+{\tilde{V}}_C(\mathbf{Y})+{\tilde{V}}_{CSB}(\mathbf{r}),$$

(79)

with

$${\tilde{V}}_{CSB}(\mathbf{r})=-\frac{1}{C_F}\frac{\sigma }{2\sqrt{\pi }{T}_g}{\mathbf{r}}^{2}I\left(0,0;\mathbf{r}\right)=-\frac{1}{C_F}\frac{\sigma }{2\sqrt{\pi }{T}_g}{\mathbf{r}}^2{\int }_0^{1}dv{\int }_0^{1}dw{e}^{-\left({\left(v+w\right)}^2\frac{{\mathbf{r}}^2}{16T_g^2}\right)}$$

(80)

$${\tilde{V}}_C(\mathbf{X})=\frac{1}{C_F}\frac{\sigma }{2\sqrt{\pi }{T}_g}{\mathbf{X}}^{2}I\left(\mathbf{X},\mathbf{X};0\right)=\frac{1}{C_F}\frac{\sigma }{2\sqrt{\pi }{T}_g}{\mathbf{X}}^2{\int }_0^{1}dv{\int }_0^{1}dw{e}^{-\left({\left(v-w\right)}^2\frac{{\mathbf{X}}^2}{4T_g^2}\right)}$$

(81)

$${\tilde{V}}_C(\mathbf{Y})=\frac{1}{C_F}\frac{\sigma }{2\sqrt{\pi }{T}_g}{\mathbf{Y}}^{2}I\left(\mathbf{Y},\mathbf{Y};0\right)=\frac{1}{C_F}\frac{\sigma }{2\sqrt{\pi }{T}_g}{\mathbf{Y}}^2{\int }_0^{1}dv{\int }_0^{1}dw{e}^{-\left({\left(v-w\right)}^2\frac{{\mathbf{Y}}^2}{4T_g^2}\right)}.$$

(82)

It turns out that ${\tilde{V}}_{CSB}(\mathbf{r})$ has the following form:

$${\tilde{V}}_{CSB}(\mathbf{r})=\frac{1}{C_F}\frac{4\sigma T_g}{\sqrt{\pi }}\left(v_l(r/T_g)-1\right),$$

(83)

where $v_l(r)$ is a rapidly decreasing function of r with $v_l(0)=1$ and $v_l(\infty )=0$. This suggests decomposing the potential according to

$$V\left(\mathbf{X},\mathbf{Y};\mathbf{r}\right)=V_C(\mathbf{X})+V_C(\mathbf{Y})+V_{CSB}(\mathbf{r}),$$

(84)

with

$$V_{CSB}\left(\mathbf{r}\right)=\frac{1}{C_F}\frac{2\sigma T_g}{\sqrt{\pi }}\left(2-\frac{{\mathbf{r}}^2}{4T_g^2}{\int }_0^{1}dv{\int }_0^{1}dw{e}^{-\left({\left(v+w\right)}^2\frac{{\mathbf{r}}^2}{16T_g^2}\right)}\right)\equiv \frac{1}{C_F}\frac{4\sigma T_g}{\sqrt{\pi }}{v}_l(r/T_g)$$

(85)

$$V_C\left(\mathbf{R}\right)={\tilde{V}}_C(\mathbf{R})-\frac{1}{C_F}\frac{2\sigma T_g}{\sqrt{\pi }}=\frac{1}{C_F}\frac{\sigma T_g}{2\sqrt{\pi }}\frac{{\mathbf{R}}^2}{T_g^2}{\int }_0^{1}dv{\int }_0^{1}dw{e}^{-\left({\left(v-w\right)}^2\frac{{\mathbf{R}}^2}{4T_g^2}\right)}-\frac{1}{C_F}\frac{2\sigma T_g}{\sqrt{\pi }}.$$

(86)

Hence, the effective interaction Hamiltonian can be split according to:

$$\begin{array}{ccc}\hfill H_{int}& =& H_{CSB}+H_C\hfill \end{array}$$

(87)

$$\begin{array}{ccc}\hfill H_{CSB}& =& \frac{1}{2}\int d^{3}x\int d^{3}y{\psi }^{†}(\mathbf{x})t_a\psi (\mathbf{x})V_{CSB}\left(\mathbf{r}\right){\psi }^{†}(\mathbf{y})t_a\psi (\mathbf{y})\hfill \end{array}$$

(88)

$$\begin{array}{ccc}\hfill H_C& =& \frac{1}{2}\int d^{3}x\int d^{3}y{\psi }^{†}(\mathbf{x})t_a\psi (\mathbf{x})K_C\left(\mathbf{X},\mathbf{Y}\right){\psi }^{†}(\mathbf{y})t_a\psi (\mathbf{y})\hfill \\ & & \mathrm{with}{K}_C\left(\mathbf{X},\mathbf{Y}\right)=V_C\left(\mathbf{X}\right)+V_C\left(\mathbf{Y}\right).\hfill \end{array}$$

(89)

$H_{CSB}$ represents a short-range interaction generating chiral symmetry breaking and the condensation of light $q\overline{q}$ pairs in the QCD vacuum, whereas $H_C$ represents a long-range interaction confining quarks inside the nucleon. In the following, the QCD vacuum and the “shifted” vacuum at finite density will be determined variationally as a BCS-like state $\left|\phi (\rho )\right.〉$. It is a simple matter to check that $\langle \phi (\rho )|H_C|\phi (\rho )\rangle =0$. Hence, the confining interaction in the presence of a string junction does not contribute to the chiral effective potential. Since this confining interaction exists only in the presence of an external source, i.e., a string junction, we can consider it to act only on valence quarks and not on the Dirac sea quarks. The same mechanism was described in ref. [68] for the heavy–light quark system, with a heavy quark seen as an external source located in ${\mathbf{r}}_0$.

3.3.1. The Self-Energy Kernel

The short-range $\mathbf{r}$-dependent piece $H_{CSB}$ will generate dynamical chiral symmetry breaking. One numerically finds that $\int d^{3}r{v}_l(r/T_g)={\pi }^{3/2}{T}_g^3{\mu }^3$ with ${\mu }^3=30$. In Figure 5, we compare the exact form with the approximate Gaussian form:

$$v_l^G(r/T_g)=e^{-\frac{r^2}{{\mu }^2{T}_g^2}}=e^{-{\left(\frac{r}{3.1T_g}\right)}^2}.$$

(90)

This interaction can be rewritten as:

$$\begin{array}{ccc}\hfill V_{CSB}\left(\mathbf{r}\right)& =& \frac{1}{C_F}\frac{4\sigma T_g}{\sqrt{\pi }}{\pi }^{3/2}{T}_g^3{\mu }^3\frac{v_l(r/T_g)}{\int d^3{r}^{\prime }{v}_l(r^{\prime }/T_g)}=\frac{1}{C_F}4\pi \sigma T_g^4{\mu }^3\frac{v_l(r/T_g)}{\int d^3{r}^{\prime }{v}_l(r^{\prime }/T_g)}\hfill \\ & \simeq & \frac{1}{C_F}4N_c{N}_f{G}_1\delta (\mathbf{r})\mathrm{with}4N_c{N}_f{G}_1\equiv 4\pi \sigma T_g^4{\mu }^3.\hfill \end{array}$$

(91)

In the second line of Equation (92), we have approximated the kernel by a point-like interaction with the same total strength. By performing a Fierz transform [60,63,64,65] of the Lagrangian or Hamiltonian of this type, $({\psi }^{†}{t}_a\psi )({\psi }^{†}{t}_a\psi )$, one can check that this interaction is equivalent to a Nambu–Jona–Lasinio (NJL) model with a scalar sector coupling constant $G_1$ (we take $N_f=2,N_c=3)$:

$$G_1=\frac{4\pi \sigma T_g^4{\mu }^3}{4N_c{N}_f}=3.1416\left(\frac{\sigma ({\mathrm{GeV}}^2)}{0.2}\right){\left(\frac{T_g(\mathrm{fm})}{0.1973}\right)}^4{\mathrm{GeV}}^{-2}.$$

(92)

The reason for introducing a point-like interaction is mainly practical. Keeping the finite-range CSB interaction would imply the use of bilocal meson fields resulting from the bosonization procedure, which was used in our previous paper [23]. This cumbersome formalism would have required completely reformulating and extending our approach and would also lack its simplicity. Conversely, one can also say that this FCM approach justifies, on firm grounds, the NJL phenomenology, which is itself an important result, despite the approximations performed. The momentum space representation of this interaction can be written as:

$$\begin{array}{ccc}\hfill {\widehat{V}}_{CSB}\left(\mathbf{q}\right)& =& \frac{1}{C_F}4\pi \sigma T_g^4{\mu }^3\mathrm{\Gamma }(\mathbf{q})=\frac{1}{C_F}4N_c{N}_f{G}_1\mathrm{\Gamma }(\mathbf{q})\hfill \\ \hfill \mathrm{\Gamma }(\mathbf{q})& =& \frac{\int d^{3}rexp\left(-i\mathbf{q}·\mathbf{r}\right)v_l(r/T_g)}{\int d^{3}r{v}_l(r/T_g)}\simeq exp\left(-\frac{{\mu }^2{T}_g^2{q}^2}{4}\right)\simeq \mathrm{\Theta }\left(\mathrm{\Lambda }-q)\right).\hfill \end{array}$$

(93)

In the second line of Equation (94), we have approximated the interaction in momentum space by the Gaussian form and then by a theta function with the cutoff $\mathrm{\Lambda }$ which regularizes the equivalent NJL model. Since this correspondence is not exact, fixing the equivalent NJL cutoff is somewhat arbitrary. A possible prescription is:

$$V_{CSB}\left(r=0\right)=\frac{1}{C_F}\frac{4\sigma T_g}{\sqrt{\pi }}=\int \frac{d^{3}q}{{(2\pi )}^3}{\widehat{V}}_{CSB}\left(\mathbf{q}\right)\simeq \frac{1}{C_F}4\pi \sigma T_g^4{\mu }^3\int \frac{d^{3}q}{{(2\pi )}^3}\mathrm{\Theta }\left(\mathrm{\Lambda }-q\right).$$

(94)

This yields:

$$\mathrm{\Lambda }={\left(\frac{6\sqrt{\pi }}{{\mu }^3}\right)}^{1/3}\frac{1}{T_g}=0.71\left(\frac{0.1973}{T_g(\mathrm{fm})}\right)\mathrm{GeV}.$$

(95)

One relevant dimensionless parameter to appreciate the strength of the NJL equivalent interaction is:

$$G_1{\mathrm{\Lambda }}^2=\frac{{(36{\pi }^4{\mu }^3)}^{1/3}}{N_c{N}_f}\sigma T_g^2=7.86\sigma T_g^2=7.86{\eta }^2.$$

(96)

By considering the NJL gap equation [26] in the chiral limit,

$$1=24G_1{\int }_0^{\mathrm{\Lambda }}\frac{d^{3}p}{{(2\pi )}^3}\frac{1}{\sqrt{p^2+M_0^2}},$$

it is easy to establish that one must have $G_1>{\pi }^2/3\simeq 3.29$ to have a nontrivial solution of the gap equation, i.e., a nonvanishing constituent quark mass, $M_0$. In practice, this necessitates a sizable ${\eta }^2>0.5$.

Starting with accepted values of the string tension and gluon condensate, for orientation, we take $\sigma =0.18{\mathrm{GeV}}^2$ and ${\mathcal{G}}_2=0.025{\mathrm{GeV}}^4$, which implies that $T_g=0.286\mathrm{fm}$, and hence, $\eta =\sqrt{\sigma }{T}_g=0.615$. Consequently, the NJL parameters are $G_1=12.514{\mathrm{GeV}}^{-2}$ and $\mathrm{\Lambda }=0.488\mathrm{GeV}$. It follows that $G_1{\mathrm{\Lambda }}^2=3$ is just below the critical value for chiral symmetry breaking. Hence, as in many QCD-inspired approaches, our model, at this level of approximation, is not able to properly generate the physical broken vacuum. We thus decided to increase the cutoff parameter to $\mathrm{\Lambda }=0.604$ MeV to obtain the correct phenomenology.

We solve the NJL gap equation, $V_{\chi ,NJL}^{\prime }(\mathcal{S})=0$, with a light-quark mass $m=5.8$ MeV to obtain the vacuum constituent quark mass $\mathcal{S}=M_0$. With this set of parameters, we ignore the vector interaction and, consequently, the $\pi$-axial mixing and calculate the pion decay constant, the pion mass and the quark condensate from the GOR relation, with all formal details being given in [23]. We obtain the following results: $M_0=356.7\mathrm{MeV},F_{\pi }=91.9\mathrm{MeV},m_{\pi }=140\mathrm{MeV},\langle \overline{q}q\rangle =-{\left(241.1\mathrm{MeV}\right)}^3$. It is remarkable that by taking accepted values of the two basic QCD parameters, namely, the string tension and the gluon condensate, one recovers the NJL parameters, yielding good NJL phenomenology, or at least their order of magnitude.

We also find that the NJL parameter entering the cubic term of the chiral effective potential is significant, $C_{\chi }=0.488$, and the effective sigma mass parameter is $M_{\sigma }=716.4$ MeV. We will show below with the FCM confining potential that the scalar coupling constant (ignoring the effect of the pion cloud) is $g_S=6.52$, and the response parameter is $C=0.32$.

3.3.2. The Confining Kernel

The confining interaction is given by Equation (86). We see immediately that it has a quadratic behavior at a short distance ($R<<T_g$), whereas at a large distance, it is possible to show that it has a linear behavior:

$$\begin{array}{ccc}\hfill R<<T_g& :& {(V_C)}_S\left(\mathbf{R}\right)+\frac{1}{C_F}\frac{2\sigma T_g}{\sqrt{\pi }}=\frac{1}{C_F}\frac{\sigma }{2\sqrt{\pi }{T}_g}{R}^2\hfill \\ \hfill R>>T_g& :& {(V_C)}_L\left(\mathbf{R}\right)+\frac{1}{C_F}\frac{2\sigma T_g}{\sqrt{\pi }}=\frac{1}{C_F}\sigma \left(R-\frac{2}{\sqrt{\pi }}{T}_g\right).\hfill \end{array}$$

(97)

For the bound-state calculation, to simplify the numerical computation, we make use of an approximate expression interpolating between the short-range and long-range behaviors.

$$\begin{array}{ccc}& & {(V_C)}_I\left(\mathbf{R}\right)=F(R){(V_C)}_S\left(\mathbf{R}\right)+\left(1-F(R)\right){(V_C)}_L\left(\mathbf{R}\right).\hfill \\ & & \mathrm{with}F(R)={\left[1+exp\left(\frac{R-(6/\sqrt{\pi })T_g}{0.1T_g}\right)\right]}^{-1}.\hfill \end{array}$$

(98)

This confining interaction is displayed in Figure 6 and Figure 7 for typical values of the parameters.

3.3.3. The Bound-State Equation

In this section, we give the bound-state equation for a constituent quark with mass $M=\overline{\mathcal{S}}(\rho )$ moving in a BCS-like state representing the broken (modified) QCD vacuum and submitted to the confining potential with the origin at the string junction point. All the details for the derivation of the results given below will be given in a long forthcoming paper (referred as [NJLFCM] [78]). The eigenvalue equation for the confined bound state described by the Dirac wave function, $\left({\psi }_n\right)(\mathbf{x})$, has the form:

$$\left[-i\mathit{\alpha }·{\overrightarrow{\nabla }}_{\mathbf{x}}+\beta M\right]\left({\psi }_n\right)(\mathbf{x})+\int d\mathbf{z}\beta M(\mathbf{x},\mathbf{z})\left({\psi }_n\right)(\mathbf{z})=E_n\left({\psi }_n\right)(\mathbf{x})$$

(99)

where the mass operator is [77,78]:

$$\begin{array}{ccc}\hfill M(\mathbf{x},\mathbf{y})& =& -C_F{\gamma }_0{\mathrm{\Lambda }}_{red}(\mathbf{x},\mathbf{y}){\gamma }_0\left(V_C(\mathbf{X}+V_C(\mathbf{Y})\right)\hfill \\ \hfill -{\mathrm{\Lambda }}_{red}(\mathbf{x},\mathbf{y})& =& \int \frac{d^{3}k}{{(2\pi )}^3}{e}^{-i\mathbf{k}·(\mathbf{x}-\mathbf{y})}{\mathrm{\Lambda }}_{red}(\mathbf{k})=\int \frac{d^{3}k}{{(2\pi )}^3}{e}^{-i\mathbf{k}·(\mathbf{x}-\mathbf{y})}\frac{s_k-c_k\mathit{\gamma }·\mathbf{k}}{2},\hfill \end{array}$$

(100)

and the chiral angle ${\phi }_k$ is defined by $s_k\equiv sin{\phi }_k=M/E_k$, $c_k\equiv cos{\phi }_k=k/E_k$ and $E_k=\sqrt{k^2+M^2}$. The reduced projector ${\mathrm{\Lambda }}_{red}$ is the difference between projectors in the positive- and negative-energy solutions of the Dirac equation. In ref. [77], a kernel with the same structure as Equation (84) was used (but in the limit of infinitely thin width, $T_g\to 0$). With such a kernel structure, the solution of this Dirac equation has the form of a quasi-plane wave state:

$$\begin{array}{ccc}\hfill \left({\psi }_n\right)(\mathbf{x})& =& \int \frac{d^{3}p}{{(2\pi )}^3}{e}^{-i\mathbf{p}·\mathbf{x}}({\tilde{\psi }}_n)(\mathbf{p})\hfill \\ \hfill ({\tilde{\psi }}_n)(\mathbf{p})& =& \sqrt{\frac{1}{2}}\left(\begin{array}{c}\sqrt{1+s_p}{\mathrm{\Phi }}_n(\mathbf{p})\\ \sqrt{1-s_p}\mathbf{\sigma }·\widehat{\mathbf{p}}{\mathrm{\Phi }}_n(\mathbf{p})\end{array}\right)=\sqrt{\frac{E_P+M}{2E_P}}\left(\begin{array}{c}{\mathrm{\Phi }}_n(\mathbf{p})\\ \frac{\mathbf{\sigma }·\mathbf{p}}{E_P+M}{\mathrm{\Phi }}_n(\mathbf{p})\end{array}\right).\hfill \end{array}$$

(101)

We are looking for normalized single-quark states, $|n;ljm\rangle$, with well-defined parity and total angular momentum, and the associated spinor wave function reads:

$$\begin{array}{ccc}\hfill {\mathrm{\Phi }}_{n;ljm}(\mathbf{p})& =& \sqrt{4}\pi R_{lj}(p){\mathrm{\Phi }}_{ljm}(\widehat{\mathbf{p}})\hfill \\ \hfill {\mathrm{\Phi }}_{ljm}(\widehat{\mathbf{p}})& =& \sum _{\mu s}{Y}_{l\mu }(\widehat{\mathbf{p}}){\chi }_s\langle l\mu ,\frac{1}{2}s|jm\rangle \mathrm{with}\mathit{\sigma }·\widehat{\mathbf{p}}{\mathrm{\Phi }}_{ljm}=-{\mathrm{\Phi }}_{l^{\prime }jm},l^{\prime }=2j-l.\hfill \end{array}$$

(102)

The bound-state equation finally reduces to a Schrödinger-like equation,

$$\begin{array}{ccc}& & E_p{\mathrm{\Phi }}_{n;ljm}(\mathbf{p})+\frac{1}{2}\int \frac{d^{3}q}{{(2\pi )}^3}{\tilde{W}}_S(\mathbf{p}-\mathbf{q})(\sqrt{1+s_p}\sqrt{1+s_q}\hfill \\ & & +\sqrt{1-s_p}\sqrt{1-s_q}\mathit{\sigma }·\mathbf{p}\mathit{\sigma }·\mathbf{q}){\mathrm{\Phi }}_{n;ljm}(\mathbf{p})=\left(E_{jl}+\frac{2\sigma T_g}{\sqrt{\pi }}\right){\mathrm{\Phi }}_{n;ljm}(\mathbf{p}),\hfill \end{array}$$

(103)

where ${\tilde{W}}_S(\mathbf{p})$ is the Fourier transform of an effective one-body confining potential $W_S(\mathbf{R})$ such that

$$W_S\left(\mathbf{R}\right)=\frac{\sigma }{2\sqrt{\pi }{T}_g}{\mathbf{R}}^{2}I\left(\mathbf{R},\mathbf{R};0\right)=\frac{\sigma }{2\sqrt{\pi }{T}_g}{\mathbf{R}}^2{\int }_0^{1}dv{\int }_0^{1}dw{e}^{-\left({\left(v-w\right)}^2\frac{{\mathbf{R}}^2}{4T_g^2}\right)},$$

(104)

with asymptotic behavior:

$$R<<T_g:W_S\left(\mathbf{R}\right)=\frac{\sigma }{2\sqrt{\pi }{T}_g}{R}^2,R>>T_g:W_S\left(\mathbf{R}\right)=\sigma \left(R-\frac{2}{\sqrt{\pi }}{T}_g\right).$$

(105)

Actually, this potential deviates from the above $V_C(\mathbf{R})$ by a factor of $1/C_F$ with the constant shift $2\sigma T_g/\sqrt{\pi }$ removed. The bound-state energy can also be written as

$$\begin{array}{ccc}\hfill E_{jl}& =& \int d\mathit{x}\left({\psi }_n\right)(\mathit{x})\left[-i\mathit{\alpha }·{\overrightarrow{\nabla }}_{\mathit{x}}+\beta M+W_S(\mathit{x})\right]\left({\psi }_n\right)(\mathit{x})-\frac{2\sigma T_g}{\sqrt{\pi }}\hfill \\ & \equiv & E_{jl,kin}+E_{jl,pot},\hfill \\ \hfill E_{jl,kin}& =& \int \frac{d\mathit{p}}{{(2\pi )}^3}{E}_p{R}_{jl}^2(p),\hfill \\ \hfill E_{jl,pot}& =& \int d\mathit{r}{W}_S(r)\left(|G_{1l}{(r)|}^2+{|G_{2l^{\prime }}(r)|}^2\right)-\frac{2\sigma T_g}{\sqrt{\pi }},\hfill \end{array}$$

(106)

with

$$G_{1l}(r)=\int \frac{d\mathit{q}}{{(2\pi )}^3}\sqrt{1+s_q}{R}_{jl}(q)j_l(qr),G_{2l^{\prime }}(r)=\int \frac{d\mathit{q}}{{(2\pi )}^3}\sqrt{1-s_q}{R}_{j{l}^{\prime }}(q)j_{l^{\prime }}(qr).$$

(107)

We look for a Gaussian solution for the lowest orbital with $j=1/2,l=0,l^{\prime }=1$,

$$R_0(p)={(2\pi )}^{\frac{3}{2}}{\left(\frac{b^2}{\pi }\right)}^{\frac{3}{4}}{e}^{-b^2{p}^2/2},$$

(108)

with b being a variational parameter. One can search for a solution for any value $M(s)$ of the in-medium constituent quark mass, hence producing a density-dependent nucleon mass as a function of the value of the “nuclear physics sigma meson” field s [26]. Here, we only need the solution for the vacuum case. With the QCD inputs given above, $\sigma =0.18{\mathrm{GeV}}^2$ and ${\mathcal{G}}_2=0.025{\mathrm{GeV}}^4$, yielding $T_g=0.286\mathrm{fm}$, $\eta =\sqrt{\sigma }{T}_g=0.615$ and $M_0=356.7$ MeV, the value of b minimizing the orbital energy is $b=0.978/\sqrt{\sigma }$. The resulting contribution to the quark orbital energy is $E_0=E_{0,kin}+E_{0,pot}=1.458\sqrt{\sigma }+0.636\sqrt{\sigma }$. The nucleon wave function is just the product of the three quark orbitals properly projected onto the color singlet state with $I=J=1/2$. After center-of-mass correction, the quark core contribution to the nucleon mass is:

$$M_N^{(core)}=\sqrt{9E_0^2-\langle P^2\rangle }=3\sqrt{E_0^2-\frac{1}{2b^2}}=3.62\sqrt{\sigma }.$$

(109)

As a side remark, it is possible to perturbatively incorporate within the model (one obtains $g_A=1.24$) the pion cloud and the chromomagnetic contribution to the nucleon mass to obtain a nucleon mass very close to the physical nucleon mass. We will discuss this point in a forthcoming paper [26], but here we consider only the quark core contribution to extract the response parameters, which can obtained either semi-analytically or numerically. One finds

$$g_S=\frac{M_0}{F_{\pi }}{\left(\frac{\partial M_N^{(core)}(\mathcal{S})}{\partial \mathcal{S}}\right)}_{\mathcal{S}=M_0}=6.52,C=\frac{M_0^2}{2M_N}{\left(\frac{{\partial }^2{M}_N^{(core)}(\mathcal{S})}{\partial {\mathcal{S}}^2}\right)}_{\mathcal{S}=M_0}=0.32.$$

(110)

The value of the C parameter looks a priori small, but since $g_S$ is significantly smaller than $M_N/F_{\pi }$, its effective value entering the parameters ${\tilde{C}}_3$ (64) and ${\tilde{C}}_L$ (65) is larger, e.g., $(M_N/g_S{F}_{\pi })C=0.50$. If we include the value of $C_{\chi }=0.488$, also derived from the microscopic FCM approach [26], one has:

$${\tilde{C}}_3=\frac{M_N}{g_S{F}_{\pi }}C+\frac{1}{2}{C}_{\chi }=0.75,{\tilde{C}}_L=\frac{M_N}{g_S{F}_{\pi }}C+\frac{3}{2}{C}_{\chi }=1.23.$$

(111)

It turns out that the s field whose value lies in the range $[0÷-F_{\pi }]$ does not exactly coincide with the canonical scalar field, $s_c$, of the bosonized NJL Lagrangian. According to [20], they are related by $s=z_S{s}_c$, where the dimensionless rescaling factor $z_S$ is expressible in terms of the NJL loop integrals. With the above NJL parameters, its numerical value is $z_S=1.284$. Hence, the values of the scalar coupling constant and of the scalar mass associated with the canonical scalar field are:

$$g_{\sigma }=z_S{g}_S=8.37,m_{\sigma }=z_S{M}_{\sigma }=919MeV.$$

(112)

This is the reason for which we quote those values of $g_{\sigma }$ and $m_{\sigma }$ in the parameters entering the NN potential in Equation (43). In addition, notice that the ratio $g_S/M_{\sigma }=g_{\sigma }/m_{\sigma }$ is not affected when passing from the field s to the canonical scalar field $s_c$, and this rescaling has no effect at the Hartree level but may have a small effect on the Fock terms. It is also possible to calculate the core rms and the vertex form factors in the various (scalar, axial and vector) Yukawa channels. Here, we give the cutoff values for the equivalent monopole form factors regularizing the corresponding Yukawa coupling to the nucleon:

$${\mathrm{\Lambda }}_{\sigma }=0.86GeV,{\mathrm{\Lambda }}_{\pi }=1.01GeV,{\mathrm{\Lambda }}_v=1.13GeV$$

(113)

This simple estimate just demonstrates that they are close to the $\mathrm{\Lambda }=1\mathrm{GeV}$ introduced above in Equation (43).

4. Equation of State and Correlations from the QCD-Connected Chiral Confining Model

To prepare the nuclear matter calculations, let us summarize the origin of the various input parameters distinguishing between those coming from the QCD-connected model and those coming from hadron phenomenology:

• Concerning the parameters entering the $NN$ interaction, the scalar and pionic sectors are entirely given or strongly constrained by the QCD-connected model: $g_{\sigma }=8.37,m_{\sigma }=919\mathrm{MeV},{\mathrm{\Lambda }}_S=1\mathrm{GeV},g_A=1.26,F_{\pi }=92.4\mathrm{MeV},m_{\pi }=140\mathrm{MeV},{\mathrm{\Lambda }}_{\pi }=1\mathrm{GeV}$. Notice that at this level, the cutoff parameters ${\mathrm{\Lambda }}_S,{\mathrm{\Lambda }}_{\pi }$, as well as the cutoff in the vector channel ${\mathrm{\Lambda }}_V$, are only approximately compatible with the QCD-connected model.
  The vector-(Lorentz) tensor $NN$ sector is constrained by well-established hadron phenomenology: $g_V=7.5,m_V=783\mathrm{MeV},g_{\rho }=g_V/3,m_{\rho }=770,\mathrm{MeV},{\kappa }_{\rho }=6$. Notice that in a version of the underlying NJL including pion-axial mixing (i.e., nonvanishing $G_2$), $g_V$ should be promoted to the rank of a QCD-connected parameter. We again stress that a large value of ${\kappa }_{\rho }$ (used in the Bonn potential [37]) is required to decrease the (Wigner) tensor force without obtaining too large a D-state probability in the deuteron.
• The two parameters entering the three-nucleon force, $C=0.32$ and $C_{\chi }=0.488$, are given by the QCD-connected model.
• The parameter $q_C=670\mathrm{MeV}$ entering the pair correlation function, $f(r)=1-j_0(q_{c}r)$, is obtained from the above $NN$ interaction with an adapted G-matrix calculation preserving the UV regularization of the loop integrals entering the correlation energy (see below). Finally, a rather large value of the cutoff, ${\mathrm{\Lambda }}_{\rho }=2\mathrm{GeV}$, for the tensor coupling of the rho meson is needed to obtain a sufficiently large value of the Landau–Migdal parameter $g^{\prime }$. This is the only constraint or requirement from nuclear matter phenomenology.
• At variance with two parallel works [25,79], the lattice parameters $a_2,a_4$ are not taken as inputs but are used to test the consistency of the model.

4.1. Binding Energy of Nuclear Matter

In a previous paper [20] devoted to a discussion of the chiral properties of nuclear matter, we calculated the equation of state with the inclusion of pion (and rho) loops on top of the Hartree mean-field approximation using the first version of the model, i.e., with the L$\sigma$M chiral potential. The energy density is written as:

$$\frac{E}{V}=\epsilon =\int \frac{4d^{3}p}{{(2\pi )}^3}\mathrm{\Theta }(p_F-p)E_p^{\ast }(\overline{s})+V_{\chi }(\overline{s})+\frac{g_V^2}{2m_V^2}{\rho }^2+{\epsilon }^{Loop},$$

(114)

where ${\epsilon }^{Loop}={\epsilon }^{Fock}+{\epsilon }^{Corr}$ summarizes the loop energy coming from the Fock terms and multi-loop correlations (with pion + rho + short-range spin–isospin effective interaction). In passing, we recall that the correlation piece (see Equation (116) below) calculated in the ring approximation actually includes the long-range RPA correlations but with a spin–isospin interaction modified by the short-range correlation. The calculation presented here is essentially the same as the one in [20] but with the following modifications:

• The L$\sigma$M potential (Equation (53)) is modified by incorporating the factor $1-C_{\chi }$ into the cubic $s^3$ term with the value $C_{\chi }=0.488$ obtained in the underlying QCD-connected NJL model, derived from the FCM approach. The parameters entering the in-medium nucleon mass, $g_S$ and C, are those given above by the same underlying microscopic model, whereas in [20], $g_S=M_N/F_{\pi }$ was taken from the original L$\sigma$M formulated at the nucleonic level, and C was a free parameter.
• The parameters governing the various interaction vertices are those listed in Equation (43). In particular, the parameters relative to the scalar sector ($g_{\sigma }$ and $m_{\sigma }$) and the cutoffs regularizing the various vertices are given by or compatible with the underlying FCM approach.
• The contribution of multi-pion (and multi-rho) exchange is explicitly included in the correlation energy, which thus incorporates the effect of the two-pion (and two-rho) exchange which is not reducible to iterative processes with NN intermediate states. The effect of these diagrams with at least one $\mathrm{\Delta }$ in the intermediate state (see Figure 10) was taken into account in the construction of the $NN$ potential (Section 2.3, Equation (44)) through the introduction of a second scalar field ${\sigma }^{\prime }$ by adding the Lagrangian ${\mathcal{L}}_{{\sigma }^{\prime }}=-g_{{\sigma }^{\prime }}\overline{N}{\sigma }^{\prime }N+{\partial }^{\mu }{\sigma }^{\prime }{\partial }_{\mu }{\sigma }^{\prime }/2-m_{{\sigma }^{\prime }}^2{\sigma }^{\prime }/2$ with $g_{{\sigma }^{\prime }}=4.8$ and $m_{{\sigma }^{\prime }}$ = 550 MeV. At the mean-field level, one has ${\overline{\sigma }}^{\prime }=-g_{{\sigma }^{\prime }}{\rho }_S/m_{{\sigma }^{\prime }}^2$. The net effect of this new scalar field is the modification of the Dirac scalar mass given by Equation (3) in [20], according to

  $$M_N^{\ast }(\overline{s})=M_N+g_S\overline{s}+\frac{1}{2}{\kappa }_{NS}{\overline{s}}^3+g_{{\sigma }^{\prime }}{\overline{\sigma }}^{\prime },$$

  where $\overline{s}$ is the solution of the in-medium gap equation. However, we do not include the associated energy per nucleon, $E_{{\sigma }^{\prime }}/A=g_{{\sigma }^{\prime }}{\overline{\sigma }}^{\prime }+m_{{\sigma }^{\prime }}^2{\overline{\sigma }}^{\prime 2}/2\rho \simeq -46(\rho /{\rho }_0)$ MeV, which is contained as a part of the correlation energy. Numerically, at the saturation density, the chemical potential is $\mu =\partial \epsilon /\partial \rho =E/A=-13.6$ MeV, which is very different from the Fermi energy when the ${\sigma }^{\prime }$ self-energy is ignored, $E_F^{no{\sigma }^{\prime }}=65$ MeV, in strong violation of the Hugenholtz–Van-Hove (HVH) theorem. The ${\sigma }^{\prime }$-field energy is actually not very far numerically from the dominant linear part of the correlation energy. It follows that the corresponding self-energy contributing to the chemical potential satisfies ${\Sigma }_{{\sigma }^{\prime }}\simeq \partial {\epsilon }_{corr}/\partial \rho$. Accordingly, the Fermi energy becomes $E_F^{with{\sigma }^{\prime }}=-25$ MeV, in much better agreement with the HVH theorem. Consequently, the unique effect of the ${\sigma }^{\prime }$ field is to decrease the in-medium Dirac nucleon mass to a value slightly above 750 MeV at normal nuclear matter density, as in [20]. (We recall in passing that this lack of thermodynamic consistency induces only a small non-relativistic correction to the binding energy at normal density.)
• For the calculation of the pion and rho Fock terms (Equations (23) and (24) of [20]) and of the correlation energy (Equations (15)–(20) in [20]), we replace the longitudinal and transverse spin–isospin interaction by the one given in Equations (47)–(50) based on the Jastrow ansatz (Equation (41)) for the correlation function with $q_c=670$ MeV. This procedure has the nice feature of satisfying the Beg–Agassi–Gal theorem, hence providing a natural UV regularization of the loop integrals. In the calculation of the correlation energy, we take the same effective interaction in the $N\mathrm{\Delta }$ and $\mathrm{\Delta }\mathrm{\Delta }$ channels, with the appropriate modification of the coupling constant with the $SU(3)$ flavor prescription for the ratio $R_{N\mathrm{\Delta }}=g_{\pi N\mathrm{\Delta }}/g_{\pi NN}=\sqrt{72/25}$. This generates a unique value for the Landau–Migdal parameters $g_{NN}^{\prime }=g_{N\mathrm{\Delta }}^{\prime }=g_{\mathrm{\Delta }\mathrm{\Delta }}^{\prime }=0.59$. The rather large value of the cutoff, ${\mathrm{\Lambda }}_{\rho }=2\mathrm{GeV}$, for the tensor coupling of the rho meson in a strong rho scenario, ${\kappa }_{\rho }=6$, has been chosen to obtain a sufficiently large value of this Landau–Migdal $g^{\prime }$ parameter compatible with phenomenology [40]. With the notations of [20], one has

  $$\begin{array}{ccc}\hfill {\epsilon }_L^{Corr}& =& \frac{3}{2}\int \frac{id\omega d\mathbf{q}}{{(2\pi )}^4}[-ln\left(1-V_{LNN}{\mathrm{\Pi }}_N^0-V_{L\mathrm{\Delta }\mathrm{\Delta }}{\mathrm{\Pi }}_{\mathrm{\Delta }}^0-(V_{LN\mathrm{\Delta }}^2-V_{LNN}{V}_{L\mathrm{\Delta }\mathrm{\Delta }}){\mathrm{\Pi }}_N^0{\mathrm{\Pi }}_{\mathrm{\Delta }}^0\right)\hfill \\ & & -V_{LNN}{\mathrm{\Pi }}_N^0-V_{L\mathrm{\Delta }\mathrm{\Delta }}{\mathrm{\Pi }}_{\mathrm{\Delta }}^0]\hfill \\ \hfill {\epsilon }_T^{Corr}& =& 3\int \frac{id\omega d\mathbf{q}}{{(2\pi )}^4}[-ln\left(1-V_{TNN}{\mathrm{\Pi }}_N^0-V_{T\mathrm{\Delta }\mathrm{\Delta }}{\mathrm{\Pi }}_{\mathrm{\Delta }}^0-(V_{TN\mathrm{\Delta }}^2-V_{TNN}{V}_{T\mathrm{\Delta }\mathrm{\Delta }}){\mathrm{\Pi }}_N^0{\mathrm{\Pi }}_{\mathrm{\Delta }}^0\right)\hfill \\ & & -V_{TNN}{\mathrm{\Pi }}_N^0-V_{T\mathrm{\Delta }\mathrm{\Delta }}{\mathrm{\Pi }}_{\mathrm{\Delta }}^0],\hfill \end{array}$$

  (115)

  where the second line of the two above expressions corresponds to the subtraction of the mean-field Fock terms.

Even if we are aware that each of the various approximations or prescriptions of this multi-step approach should or would deserve a more detailed study, we do not refrain from showing the results of the calculation without a further fine-tuning of the parameters. The saturation curve is displayed in Figure 8. The binding energy at the saturation point is too small, $E_{sat}/A=-13.18$ MeV, and the saturation density is slightly too large, $\rho$ = 1.05 ${\rho }_0=0.168{\mathrm{fm}}^{-3}$, whereas the incompressibility modulus, $K_{sat}$ = 215 MeV, is acceptable, although a bit too low. Even though the saturation point is not perfect, it is rather satisfactory to arrive at decent results for the properties of nuclear matter in an approach based on a microscopic model that mainly depends on QCD inputs, namely, the string tension and the gluon correlation length (or the gluon condensate) and parameters ($g_V/m_V,{\kappa }_{\rho }$) known from well-established hadronic phenomenology. A robust conclusion is that the pion and rho loops are necessary to obtain sufficient binding, a conclusion already reached in ref. [20]. Concerning the comparison with lattice data, our model calculation yields:

$$a_2=1.175{\mathrm{GeV}}^{-1},a_4=-0.615{\mathrm{GeV}}^{-3}.$$

(116)

The value of $a_2$ is compatible with lattice data and yields a nonpionic contribution to the sigma commutator ${\sigma }^{(s)}=a_2{M}_{\pi }^2=23.5$ MeV. The value of $a_4$, although a bit high, is essentially compatible with lattice data in the sense that it is much smaller than the value obtained in the simplest linear sigma model ignoring the nucleonic response (i.e., $C=0$) and the NJL correction (i.e., $C_{\chi }=0$), for which $a_4\simeq -3.5{\mathrm{GeV}}^{-3}$. Hence, the model generates the strong compensation required by lattice data from effects governing the three-body repulsive force needed for the saturation mechanism: compare Equations (65) and (64).

At variance with our original work [20], we should also introduce the Fock contribution to the binding energy for the scalar, the omega and and the (Lorentz vector piece of) rho meson exchanges, i.e., as in refs. [21,79]. For orientation, we can use a non-relativistic local approximation, replacing the exchange momentum in the meson propagators and vertex form factors by an averaged value ${\overline{q}}^2=6k_F^2/5$, which is valid if the Fermi momentum is small compared to the meson masses:

$$\begin{array}{ccc}\hfill \frac{E_F^{(s)}}{A}& =& \frac{g_{\sigma }^2}{16}\frac{{\mathrm{\Lambda }}_S^2({\overline{q}}^2)}{m_{\sigma }^2+{\overline{q}}^2}\frac{{\rho }^2+{\rho }_S^2}{\rho }\hfill \end{array}$$

(117)

$$\begin{array}{ccc}\hfill \frac{E_F^{(\omega )}}{A}& =& -\frac{g_V^2}{8}\frac{{\mathrm{\Lambda }}_V^2({\overline{q}}^2)}{m_V^2+{\overline{q}}^2}\frac{2{\rho }_S^2-{\rho }^2}{\rho },\frac{E_F^{(\rho )}}{A}=-\frac{g_V^2}{24}\frac{{\mathrm{\Lambda }}_V^2({\overline{q}}^2)}{m_{\rho }^2+{\overline{q}}^2}\frac{2{\rho }_S^2-{\rho }^2}{\rho }.\hfill \end{array}$$

(118)

In addition, this Fock term is included perturbatively in the sense that we remain in the Hartree basis, since we neglect the Fock contribution to the scalar self-energy of the nucleon. It turns out that there is a strong compensation between the scalar and omega Fock energy, but the rho Fock contribution provides a sizable attraction, and the net result at normal nuclear density is $E_F/A=-3.31$ MeV. We also introduce the effect of short-range correlations by subtracting from the Hartree and Fock contributions the same quantities but adding $q_C^2$ to the exchange momentum squared in the meson propagators and form factors. The net result at normal nuclear matter density is extremely small, $E_{SRC}/A=-0.03$ MeV. Of course, the effect of these Fock terms in the presence of SRCs is the significant modification of the saturation properties: the saturation density is now much too large, ${\rho }_{sat}=1.2{\rho }_0$, whereas the binding energy becomes $E_{sat}/A=-16.9$ MeV. To solve this problem, one certainly has to first incorporate the effect of the Fock terms in the nucleon self-energy, hence working in the fully self-consistent Hartree–Fock basis with a resulting modification of the Dirac effective mass. If we decide not to change the parameters that impact the bare NN interaction, one can also advocate a possible improvement in the FCM approach, which is flawed with various shortcomings, thus allowing variations in the C and $C_{\chi }$ parameters. Just to give a flavor of what might result, we modify the parameters governing the repulsive three-body force by taking a value of $C=0.4$, compatible with other confining models, and by decreasing the NJL parameter $C_{\chi }$ by $5\%$. The resulting saturation point (see Figure 9) is slightly improved, $E_{sat}/A=-14.85$ MeV, ${\rho }_{sat}=1.05{\rho }_0=0.168{\mathrm{fm}}^{-3}$, $K_{sat}=230$ MeV, whereas the lattice parameter $a_4=-0.41{\mathrm{GeV}}^{-3}$ is compatible with the lattice results.

It nevertheless remains true that further works are needed to adjust the parameters, with the QCD-connected model acting as a guideline, to address the question of the EOS at a higher density, particularly for the neutron star EOS.

4.2. The Interrelated Role of the Short-Range Correlations

Short-range correlations (SRCs) play a crucial role in the loop energy. The Fock contribution can be split into pion and rho contributions. At normal nuclear matter density, the presence of the SRC reduces the repulsive pion Fock term from $E_F^{\pi }/A=11.71$ MeV to $E_F^{\pi +SRC}/A=2.87$ MeV and transforms the repulsive rho Fock term, $E_F^{\rho }/$A = 8.86 MeV, into a strongly attractive contribution, $E_F^{\rho +SRC}/A=-20.5$ MeV. Consequently, the total Fock term, $E_F/A=-17.6$ MeV, appears as an important contribution to the binding energy of nuclear matter. The same conclusion has been reached in a parallel work, where all the Fock terms have been incorporated into a full HF approach [79]. The combined effects of the composite nature of the nucleon embedded in the form factors and of the SRC control the magnitude of the RPA correlation integrals. In particular, the UV behavior of the spin–isospin interaction linked to the Jastrow correlation function provides a natural regularization of the loop integral as a consequence of the Beg–Agassi–Gal theorem. It is interesting to make a comparison with the chiral perturbation calculation in ref. [80], where the effect of a short-range correlation is absent. In this work, the contribution of the iterated pion exchange is found to be $-68$ MeV, and the short-range physics is described by a unique cutoff regularizing the pion loop. This has to be compared with our result where $E_{corr}/A=-32.3$ MeV at normal nuclear matter density. This calculation of the correlation energy fully incorporates the recoil correction and the energy dependence of the spin–isospin interaction, which means that the $q^2$ appearing in the pion and rho propagators is systematically replaced by $q^2-{\omega }^2$, and the integral $\int id\omega d\mathbf{q}$ becomes $-\int dzd\mathbf{q}$ after a Wick rotation (see [20,22]), hence transforming $q^2-{\omega }^2$ into $q^2+z^2$.

In the following discussion, we consider the case of $\rho =0.8{\rho }_0$, i.e., $k_F=245$ MeV, appropriate for the carbon nucleus. The full correlation energy per nucleon is $E_{corr,full}=-27$ MeV. If we limit the calculation of the correlation energy to the two-loop order, which is essentially equivalent to second-order perturbation theory,

$$\begin{array}{ccc}\hfill {\epsilon }_L^{Corr}& \simeq & \frac{3}{2}\int \frac{id\omega d\mathbf{q}}{{(2\pi )}^4}\left[V_{LNN}^2{\mathrm{\Pi }}_N^{02}+2V_{LN\mathrm{\Delta }}^2{\mathrm{\Pi }}_N^0{\mathrm{\Pi }}_{\mathrm{\Delta }}^0+V_{L\mathrm{\Delta }\mathrm{\Delta }}^2{\mathrm{\Pi }}_{\mathrm{\Delta }}^{02}\right]\hfill \\ \hfill {\epsilon }_T^{Corr}& \simeq & 3\int \frac{id\omega d\mathbf{q}}{{(2\pi )}^4}\left[V_{TNN}^2{\mathrm{\Pi }}_N^{02}+2V_{TN\mathrm{\Delta }}^2{\mathrm{\Pi }}_N^0{\mathrm{\Pi }}_{\mathrm{\Delta }}^0+V_{T\mathrm{\Delta }\mathrm{\Delta }}^2{\mathrm{\Pi }}_{\mathrm{\Delta }}^{02}\right],\hfill \end{array}$$

(119)

its value is slightly increased to $E_{corr}=-30.7$ MeV. At this level, one can clearly distinguish the various contributions to the correlation energy per nucleon, corresponding to the presence of $NN,N\mathrm{\Delta },\mathrm{\Delta }\mathrm{\Delta }$ in the intermediate states (see Figure 10):

$$E_{corr}=E_{corr}^{NN}+E_{corr}^{N\mathrm{\Delta }}+E_{corr}^{\mathrm{\Delta }\mathrm{\Delta }}=-11.8-13-5.9=-30.7\mathrm{MeV}.$$

(120)

If we remove the energy dependence of the spin–isospin interaction, i.e., in a static approximation, one can show that this two-loop correlation energy takes the simple form

$$E_{corr,st}=\frac{1}{2A}\sum _{h_1,h_2,p_1,p_2}\frac{{\left|<p_1{p}_2|G_{\sigma \tau }|h_1{h}_2>\right|}^2}{{ϵ}_{h_1}+{ϵ}_{h_2}-{ϵ}_{p_1}-{ϵ}_{p_2}},$$

(121)

where $G_{\sigma \tau }$ is the static spin–isospin interaction given by Equations (47)–(50). Apart for the energy denominator, this expression has the same structure as the expression for the mean depopulation of the nucleon Fermi sea discussed above in Section 2.2:

$$\mathrm{\Delta }n=\frac{1}{2A}\sum _{h_1,h_2,p_1,p_2}{\left|\frac{<p_1{p}_2|G_{\sigma \tau }|h_1{h}_2>}{{ϵ}_{h_1}+{ϵ}_{h_2}-{ϵ}_{p_1}-{ϵ}_{p_2}}\right|}^2\equiv \kappa .$$

(122)

but omitting the Pauli exchange term, which is consistent with the ring summation. This allows the simultaneous discussion, at the same level of approximation, of the role of the SRC in the loop energy affecting the equation of state and in the depopulation of the Fermi sea, $\mathrm{\Delta }n$, which is an excellent indicator of the effect of the SRC on the nuclear response functions and which also coincides with the wound integral of the Brueckner theory. In the two above expressions, the particle p states can be either a nucleon above the Fermi sea or a $\mathrm{\Delta }$ state. The numerical calculation of the RPA-2p-2h correlation energy can be performed by using Equations (48)–(50) in our previous work [13] by just replacing the squared energy denominators by the (always negative) single-energy denominators. In this static approximation, the various contributions to the long-range correlation energy are:

$$E_{corr,st}=E_{corr,st}^{NN}+E_{corr,st}^{N\mathrm{\Delta }}+E_{corr,st}^{\mathrm{\Delta }\mathrm{\Delta }}=-11.1-14.8-8.2=-34.2\mathrm{MeV}$$

(123)

which is slightly above the nonstatic calculation. As a byproduct, we can also estimate the contribution to the correlation energy of the scalar ($s,{\sigma }^{\prime }$) and vector meson ($\omega$) exchanges that we have not considered explicitly before. This amounts to replacing, in the dotted lines of Figure 10, the spin–isospin interaction by the $\sigma +\omega$ interaction in the non-relativistic case; one finds numerically:

$$E_{corr,st}^{\sigma +\omega }=-3.6\mathrm{MeV},$$

(124)

which is probably overestimated since we did not incorporate the in-medium dropping of the scalar coupling constant.

We now come to the calculation of the Fermi sea depopulation and the kinetic energy per nucleon as a complement of the model results given in Section 5.2 of [13]. From baryon number conservation, the total depopulation per nucleon, $\mathrm{\Delta }n$, is equal to the fraction of states (nucleon or $\mathrm{\Delta }$) outside of the Fermi sea. One obtains:

$$\begin{array}{ccc}\hfill \mathrm{\Delta }n& =& \mathrm{\Delta }{n}^N+\mathrm{\Delta }{n}^{\mathrm{\Delta }}=0.113+0.019=0.132.\hfill \end{array}$$

(125)

The proportion of $\mathrm{\Delta }$ in the ground state induced by the short-range correlation is thus $2\%$. The contribution of the probably overestimated scalar+vector exchange is smaller than the pion+rho exchange:

$$\mathrm{\Delta }{n}^{\sigma +\omega }=0.059,$$

(126)

and the full depopulation of the Fermi sea (or the BHF wound integral) is $\mathrm{\Delta }n=0.19$, which is not far from the estimate, $\mathrm{\Delta }n\simeq 0.15÷0.2$, of references [6].

The calculation of the (non-relativistic) kinetic energy per nucleon, another test of the role of the SRC, can be performed using Equations (51)–(54) in ref. [13]. One finds:

$$<t>=\frac{3}{5}\frac{k_F^2}{2M}-<t{>}^{depop}+<t{>}_{tail}^N+<t{>}^{\mathrm{\Delta }}=19-2.72+10.3+5.9=32.5\mathrm{MeV}.$$

(127)

One can notice the sizable contribution of the deltas present in the ground state. The contribution of the scalar+vector exchange is small:

$$<t{>}^{\sigma +\omega }=-1.25+3.95=2.7\mathrm{MeV}.$$

(128)

The total kinetic energy is $<t>=35.2$ MeV, where a significant part $16\mathrm{MeV}$ comes from SRCs depopulating the Fermi sea. This value is also within the estimate, $<t>\simeq 35\mathrm{MeV}$, of [6,7,8], which is needed to explain the EMC effect.

5. Conclusions and Perspectives

The chiral confining model, with inputs linked to genuine QCD quantities (string tension, gluon condensate) and well-established hadronic phenomenology (${\kappa }_{\rho }\sim 6$), is able to generate, via an adapted G-matrix approach, a pair correlation function that ensures the natural UV regularization of the pionic correlation energy. It also gives, without further fine-tuning, reasonable results for the chiral properties of the nucleon ($a_2$ and $a_4$ parameters extracted from lattice data) and saturation properties of nuclear matter. In addition, it allows the simultaneous discussion of the role of SRCs in the loop energy affecting the equation of state and in the depopulation of the Fermi sea, $\mathrm{\Delta }n$, which is an excellent indicator of the effect of the SRC on the nuclear response. It is remarkable that one obtains good results for $\mathrm{\Delta }n$ and the kinetic energy per nucleon $<t>\sim 35$ MeV, which is needed to reproduce the celebrated EMC effect. Conversely, one can remark that these values of the Fermi sea depopulation and of the kinetic energy per nucleon induced by the SRC are compatible with a large value of the rho meson tensor coupling (${\kappa }_{\rho }\sim 6$ or $c_{\rho }\sim 2$) and of the associated cutoff (${\mathrm{\Lambda }}_{\rho }\sim 2GeV$), which, in turn, induces a large binding energy coming from pion and rho loops. This very interesting link can be captured by just looking at Equations (121) and (122).

Concerning the QCD-connected chiral confining model specifically, which is at the heart of the whole approach, the effect of the NJL-like enriched chiral effective potential significantly improves the compatibility of the lattice data with the model calculations of the dimensionless response parameter C. Another important result discussed in detail in [26] is the existence of a particular combination of C and $C_{\chi }$ (${\tilde{C}}_L$) constrained by LQCD, and a closely related combination (${\tilde{C}}_3=C+C_{\chi }/2$) governing the repulsive three-nucleon force ensuring the saturation mechanism. Hence, the fundamental QCD theory and nuclear matter modeling are linked by, on the one hand, the LQCD data $a_2$ and $a_4$ and, on the other hand, what we have called the “QCD connected parameters”, namely, the response parameters $g_s$ and C. Indeed, these results provide a link between the chiral properties of the nucleon and the saturation mechanism and/or a link between the fundamental theory of strong interaction and nuclear matter properties, although the results presented in this paper should be enriched in various aspects of this multi-step approach.

The above discussion opens perspectives that may impact the strategy to be adopted in the near future. First, as pointed out in Section 4.1, the Fock terms, including the rearrangement terms [21], have not been consistently incorporated at the level of nucleon self-energy. Although this problem is of minor importance for the binding energy around the nuclear matter density, the use of a Hartree–Fock basis in place of the Hartree basis for the nucleon Dirac wave function may play a very important role at a high density in limiting the maximum mass of a hyperonic neutron star, as pointed out in ref. [24]. The passage to this Hartree–Fock basis can be implemented without formal difficulty but at the expense of increasing numerical complexity [21,79]. However, the incorporation of this thermodynamic consistency (HVH theorem) for the correlation energy certainly necessitates additional formal developments. Second, there are still improvements to be made in several directions of the theoretical treatment of the FCM-QCD model, as discussed and suggested in Section 3.3. Given the remaining theoretical uncertainties, one possible strategy is to use a Bayesian method with QCD-connected parameters and lattice data ($a_2,a_4$) parameters, implemented with their associated uncertainties, as prior input variables. This may be implemented with the methodology used in two recent papers of the Lyon group [25,79].