Relativistic self-energy decomposition of nuclear symmetry energy and equation of state of neutron matter within QCD sum rules

TL;DR

This paper uses QCD sum rules to analyze the relativistic self-energy decomposition of nuclear symmetry energy and the equation of state of neutron matter, linking non-perturbative QCD vacuum properties to dense nuclear matter.

Contribution

It demonstrates that QCD sum rules can quantitatively explore asymmetric nuclear matter properties and connect them to the non-perturbative QCD vacuum, especially at lower densities.

Findings

QCDSR effectively describes properties of asymmetric nuclear matter at low densities.

The approach links the EOS of neutron matter to QCD vacuum characteristics.

Provides insights into the uncertainties of symmetry energy from a QCD perspective.

Abstract

Abstract (abridged edition): Properties of the relativistic nucleon self-energy decomposition of the symmetry energy as well as the equation of state (EOS) of pure neutron matter (PNM) are explored systematically within the QCD sum rules (QCDSR). Our results in the present work have demonstrated that the QCDSR approach can be used to explore the properties of asymmetric nuclear matter (ANM) in a quantitative manner, at least in lower density region. The QCDSR approach establishes a bridge connecting the EOS of ANM and the non-perturbative QCD vacuum, and thus provides a useful way to understand the properties of dense nucleonic matter from non-perturbative QCD vacuum. These studies are helpful to investigate the QCD origins about the uncertainties of nucleonic matter properties, e.g., the uncertainties of the symmetry energy. On the other hand, the exact knowledge on the EOS of ANM extracted from experiments, observations and model-independent calculations is also very useful for understanding the quark/gluon condensates in nuclear medium, which can provide important information on the chiral symmetry restoration phase transition in nuclear matter as well as the in-medium effects of hadron properties.

Equations

$$E(\rho ,\delta )\approx E_0(\rho )+E_{\text{sym}}(\rho ){\delta }^2+\mathcal{O}({\delta }^4),$$

$$E_{\mathrm{sym}}(\rho )$$

$$E_0(\rho )\approx E_0({\rho }_0)+\frac{1}{2}{K}_0{\chi }^2+\mathcal{O}({\chi }^3),$$

$$K_0={9{\rho }_0^2\frac{{\text{d}}^2{E}_0(\rho )}{\text{d}{\rho }^2}|}_{\rho ={\rho }_0},$$

$$E_{\text{sym}}(\rho )\approx E_{\text{sym}}({\rho }_{\text{r}})+L({\rho }_{\text{r}}){\chi }_{\text{r}}+\mathcal{O}({\chi }_{\text{r}}^2),$$

$$L({\rho }_{\text{r}})={3{\rho }_{\text{r}}\frac{\text{d}{E}_{\mathrm{sym}}(\rho )}{\text{d}\rho }|}_{\rho ={\rho }_{\text{r}}}.$$

$$\Sigma (|\text{k}|,k^0)=$$

$$=$$

$$+\vec{\gamma }\cdot {\text{k}}^0{\Sigma }_{\text{K}}(|\text{k}|,k^0),$$

$$G(k)=G^0(k)+G^0(k)\Sigma (k)G(k)$$

$$G^{-1}(k)={\gamma }_{\mu }[k^{\mu }+{\Sigma }^{\mu }(k)]-[M+{\Sigma }_{\text{S}}(k)].$$

$$M^{\ast }=$$

$${\text{k}}^{\ast }=$$

$$e^{\ast }=$$

$$k^{\ast \mu }=$$

$$G_{\text{F}}(k)=\frac{{\gamma }^{\mu }{k}_{\mu }^{\ast }+M^{\ast }}{k^{\ast ,2}-M^{\ast ,2}+i{0}^{+}},$$

$$G_{\text{D}}(k)=\frac{i\pi }{e^{\ast }}({\gamma }^{\mu }{k}_{\mu }^{\ast }+M^{\ast })\delta (k^0-e)\Theta (k_{\text{F}}-|\text{k}|).$$

$$e=$$

$$=$$

$$+{\Sigma }_{\text{V}}(|\text{k}|,e),$$

$$e(|\text{k}|)=$$

$$+{\Sigma }_{\text{V}}(|\text{k}|).$$

$$e_J(\rho ,\delta ,|\text{k}|)=$$

$$\displaystyle+\left[M+\Sigma_{\textrm{S}}^{J}(\rho,\delta,|\textbf{k}|)\right]^{2}\Big{]}^{1/2}+\Sigma_{\textrm{V}}^{J}(\rho,\delta,|\textbf{k}|),$$

$$e_J(\rho ,\delta ,|\text{k}|)=$$

$$+{\Sigma }_{\text{V}}^J(\rho ,\delta ,|\text{k}|).$$

$$E_0(\rho )=$$

$$E_{\text{n}}(\rho )=$$

$$E_{\text{sym}}(\rho )=$$

$$+\frac{1}{2}\left(\frac{M_0^{\ast }}{e_{\text{F}}^{\ast }}{\Sigma }_{\text{sym}}^{\text{S}}+{\Sigma }_{\text{sym}}^{\text{V}}\right),$$

$${\mathcal{A}}_{ab}^{\mu }=A^{A\mu }{t}_{ab}^A,$$

$$[t^A,t^B]=i{f}^{ABC}{t}^C,\text{tr}(t^A)=0,\text{tr}(t^A{t}^B)=\frac{1}{2}{\delta }^{AB},$$

$${\mathcal{G}}_{\mu \nu }=G_{\mu \nu }^A{t}^A\equiv D_{\mu }{\mathcal{A}}_{\mu }-D_{\nu }{\mathcal{A}}_{\mu },$$

$${\mathcal{G}}_{\mu \nu }=\frac{i}{g_{\text{s}}}[D_{\mu },D_{\nu }],$$

$$G_{\mu \nu }^A={\partial }_{\mu }{A}_{\mu }^A-{\partial }_{\nu }{A}_{\nu }^A+g_{\text{s}}{f}^{ABC}{A}_{\mu }^B{A}_{\nu }^C.$$

$${\Pi }_{\alpha \beta }(q)\equiv i\int {\mathrm{d}}^{4}x{e}^{iqx}⟨0|\text{T}{\eta }_{\alpha }(x){\overline{\eta }}_{\beta }(0)|0⟩,$$

$${\eta }_{\text{p}}={\epsilon }_{abc}({\text{u}}_a^{\text{T}}\text{C}{\gamma }_{\mu }{\text{u}}_b){\gamma }_5{\gamma }^{\mu }{\text{d}}_c,$$

Full text

I Introduction

The investigation of the equation of state (EOS) of isospin asymmetric nucleonic matter (ANM) from quantum chromodynamics (QCD) or QCD-based effective theories is one of longstanding issues in nuclear physics Bra14 ; ACCNP . The exact knowledge on the EOS of ANM provides important information on the in-medium nucleon-nucleon effective interactions, which play a central role, for instance, in understanding the structure and decay properties of finite nuclei Rin80 ; Bla80 ; Ben03 , the dynamical processes in nuclear reactions LiBA98 ; Dan02 ; Bar05 ; Ste05 ; Che07a ; LCK08 ; Tsa12 , the structure and evolution of neutron stars as well as the mechanism of core-collapse supernova explosion Gle00 ; Lat04 ; Lat12 ; SNRev ; Hem12 ; Oze16 ; Wat16 ; Oer17 ; Bal16 ; Cow19 , and the gravitational waves from binary neutron star merger or black-hole neutron star merger in the multimessenger era Abb16 ; Abb17 ; Bau12 ; Tak14 ; Ber15 ; Bos18 ; Fat18 ; Abb18 ; Mos18 ; De18 ; Abb19 ; Bru18 ; Bai17 ; Due19 ; NZhang18 ; NBZhang19 ; NBZhang19a ; ZhouY19 . Conventionally, the EOS of ANM defined as the binding energy per nucleon, is expanded around the symmetric nucleonic matter (SNM) as $E(\rho,\delta)\approx E_{0}(\rho)+E_{\textrm{sym}}(\rho)\delta^{2}+\mathcal{O}(\delta^{4})$ where $\rho=\rho_{\textrm{n}}+\rho_{\textrm{p}}$ and $\delta=(\rho_{\textrm{n}}-\rho_{\textrm{p}})/\rho$ are respectively the nucleon density and isospin asymmetry of the system in terms of the neutron and proton densities $\rho_{\textrm{n}}$ and $\rho_{\textrm{p}}$. In the above expansion, $E_{0}(\rho)\equiv E(\rho,0)$ is the EOS of SNM, while $E_{\textrm{sym}}(\rho)$ is the nuclear symmetry energy. The symmetry energy $E_{\textrm{sym}}(\rho)$ can be generally decomposed into kinetic and potential parts, i.e., $E_{\textrm{sym}}(\rho)=E_{\textrm{sym}}^{\textrm{kin}}(\rho)+E_{\textrm{sym}}^{\textrm{pot}}(\rho)$. Moreover, the potential part $E_{\textrm{sym}}^{\textrm{pot}}(\rho)$ can be further decomposed generally into several terms originated from the Lorentz structure of the nucleon self-energies Cai12 . A thorough understanding on the origin and properties of each part of the symmetry energy is useful in both nuclear physics and astrophysics EPJA . For instance, in simulating heavy-ion reactions using transport models one needs as an input the potential symmetry energy of quasi-nucleons LCK08 ; LiBA97 ; ChenLW05 ; Che14 ; Li15 ; Cai16 ; RWang19 , while the kinetic symmetry energy is found to strongly affect the critical formation density of the resonance state $\Delta(1232)$ in neutron stars Cai15a ; Dra14 . Furthermore, it has been widely discussed recently that the kinetic symmetry energy is closely related to the nucleon-nucleon short range correlations (SRC) Hen17 ; Due18 ; Sch19 ; Tor18 ; Wei18 ; Wei19 ; Fom17 ; CXu11x ; CXu13x ; Hen15b ; CaiLi15a ; CaiLi16a ; CaiLi16b ; CaiLi17 ; Yong17 ; Dut18 , indicating the very fundamental nature of the symmetry energy. After hard efforts made in the last few decades, the magnitude of the total symmetry energy at the nuclear saturation density $\rho_{0}\approx 0.16\,\textrm{fm}^{-3}$ is now best known to be around $32\pm 3\,\textrm{MeV}$ Che12a ; LiBA13 ; Oer17 ; BALi17xx ; Che17xx .

The determination on the EOS of ANM based on the phenomenological approaches such as the non-relativistic Skyrme–Hartree–Fock (SHF) models Sto06 ; Cha09 ; Zha16xx ; RWang18 and the relativistic mean field (RMF) models Ser86 ; Ser97 ; Rei89 ; Rin96 ; Men06 ; ChenLW07 has made great success in recent years. However, the origin of each individual term of the EOS in these phenomenological approaches is usually blurry. For example, the contribution from the $\rho$ meson to the symmetry energy largely depends on whether or not the $\delta$ meson and the cross interactions between the $\rho$ meson and isoscalar mesons are included in the model for the nonlinear RMF model Cai12 , indicating the contribution of the $\rho$ meson to the symmetry energy is only effective. On the other hand, from the viewpoint of many body theories Ser86 , the decomposition and analyses of a quantity via the Green’s functions or more precisely, according to the Lorentz structure of the quantity itself, is physical. In this sense, any effective method, especially that based on the QCD, encapsulating the proper Lorentz structure of the nucleon self-energies to investigate the EOS of ANM will be appealing.

The QCD sum rules (QCDSR) method Shi79a ; Shi79b ; Shi79c provides an important non-perturbative QCD approach to explore the properties of nucleonic matter (see, e.g., refs. Dru90 ). Intuitively, when the QCD coupling constant is small at high energies/small distances, the theory becomes asymptotically free, guaranteeing the applicability of perturbative calculations. As the energy scale decreases, the coupling constant of the theory becomes large, perturbative methods break down eventually and non-perturbative effects emerge. Among these effects, the most important is the appearance of the quark/gluon condensates. The QCDSR is actually based on some duality relations. More specifically, the basic idea of QCDSR for nucleonic matter calculations Dru90 ; Coh91 ; Coh92 ; Fur92 ; Jin93 ; Jin94 ; Dru04 ; Dru12 ; Dru17a ; Dru17b ; Dru19 ; Jeo13 ; Jeo16 ; Mar18 ; Son18 ; Coh95 ; Rei85 ; Nar89 ; Shi92 ; Dos94 ; Lei97 ; Nar04 ; Iof06 ; Iof11 ; Gub13 is to relate the condensates to the nucleon self-energies using the operator product expansion (OPE) technique Wil69 , where information on the nucleon self-energies is introduced via nucleon-nucleon correlation functions. On the other hand, the EOS of ANM can be obtained through the self-energy decomposition by analyzing the general Lorentz structure of the nucleon self-energies Cai12 . Within the QCDSR method, the exact information on the nucleon self-energies and nucleonic matter EOS can thus provide useful constraints on the in-medium quark/gluon condensates and vice verse. It should be noted that the in-medium quark condensates provide an order parameter of spontaneous chiral symmetry breaking in QCD. Practically, the QCDSR method is expected to work well at lower densities/momenta where effects of the complicated and largely unknown high mass-dimensional condensates as well as continuum effects are not important.

In this work, we use the QCDSR method to investigate the EOS of isospin asymmetric nucleonic matter through the nucleon self-energy decomposition formulae Cai12 based on the Hugenholtz–Van Hove (HVH) theorem Hug58 , and especially focus on the nuclear symmetry energy and the EOS of pure neutron matter (PNM). It is necessary to point out that the symmetry energy was also studied recently by the QCDSR through a specific approach mapping the nucleon self-energies to the symmetry energy Jeo13 . Moreover, some main results on the EOS of PNM via the QCDSR were already reported in ref. Cai17QCDSR , and more details of the calculations in ref. Cai17QCDSR will be given in the present work.

Before going to present the details on the QCDSR method and its application in EOS of ANM, we would like to give a brief overview on several celebrated or potentially useful QCD-based approaches to the study on the EOS of ANM. These approaches include: 1). Chiral perturbation theory (ChPT) Wei90 ; Wei91 ; Wei92 ; Gas84 ; Gas85 ; Gas88 ; Mei93 ; Ber95 ; Eck95 ; Pic95 ; Bed02 ; Epe06 ; Ber07 ; Bij07 ; Bur07 ; Ber08 ; Epe09 ; Bog10 ; Mac12 ; Sch12 ; Hol13 ; Ham13 ; Heb15 ; Hol16 . The basic idea of the ChPT is based on the general effective field theory Geo93 ; Wei79 ; Pol92 and the ChPT has become a well developed tool for systematically dealing with the EOS, e.g., of PNM Lut02 ; Kai02 ; Heb10 ; Lac11 ; Tew13 ; Kru13 ; Kru13a ; Dri14 ; Dri16 ; Hag14 ; Hol17 ; Eks18 ; Dri19 . However, due to the large uncertainties, e.g., on the many-nucleon forces Ham13 ; Heb15 and the nature of the method itself, the applicability of this approach is essentially limited to the relatively low density region. 2). Purturbative QCD (pQCD) Col11 ; Dok91 ; Ste95 . The pQCD approach is often used to explore the high-density behavior of the dense nucleonic/quark matter Gro81 ; Shu84 ; Sch98 ; Alf98 ; Rap98 ; Sho03 ; Hua03 ; Alf08 ; Fuk11 ; Fre77a ; Fre77b ; Fre77c ; Bal78 ; Kap79 ; Kur10 ; Jeo16b , where the QCD coupling constant $\alpha_{\textrm{s}}=g_{\textrm{s}}^{2}/4\pi$ is generally small, indicating the possibility of the perturbative schemes applied in dense matter at extremely high densities Kap06 ; Bel96 . 3). The large $N_{\textrm{c}}$ method Hoo74 ; Hoo74a ; Wit79 ; Wit79a ; Col85 ; Man98 . The QCD becomes solvable when the color number $N_{\textrm{c}}$ is large. This is an effective field theory based on the expansion of $N_{\textrm{c}}^{-1}$, and by using the techniques developed for large $N_{\textrm{c}}$ theories, information on the in-medium nucleon self-energies can be obtained Coh11 , which is potentially useful for the further study on nuclear matter EOS. 4). Method of Dyson-Schwinger equation (DSE) Dys49 ; Sch51 ; Sch51a ; Rob94 ; Rob00 ; Liu01 ; LiuY02 ; Yua06 ; Bas12 ; Zon05 ; Zon08 . The DSE method is one of non-perturbative approaches, in which the resummation techniques Bra90 ; Bra90a ; Fen90 are adopted to make effective approximations to the QCD Lagrangian. 5). Functional renormalization group (FRG) approach Wet93 ; Pol84 . The FRG method has been applied to investigate the EOS of ANM Dre13 ; Dre14 ; Dre15 ; Dre17 , where the quantum fluctuations are included non-perturbatively Bag01 ; Ber02 ; Pol03 ; Ros12 via the renormalization group techniques Pol84 ; Dre13 ; Dre14 ; Dre15 ; Dre17 . 6). The Skyrme model Sky61 . In the Skyrme model, the nucleons are described as a topological solitons (i.e., skyrmions) in a meson field theory and thus the nucleon properties in nuclear matter can be obtained Kal97 ; Rak98 , which can be further used to obtain nuclear matter EOS Lee11 ; Ma17 . 7). Last but not least, numerical approaches such as lattice QCD EPJA1 ; Lee08 ; Lee17 , quantum Monte Carlo (QMC) simulations Car15 ; Gan15 and several different types of ab-initia methods Dru10 ; Her16 ; Naz16 are important for the understanding of the root properties of finite nuclei and the dense nucleonic matter.

Early successes of QCDSR in nucleonic matter calculations can be traced back to the prediction on the large nucleon Lorentz covariant self-energies Coh91 . The present work is a natural generalization to the investigation on the EOS of ANM. Besides the prediction on the EOS of PNM by the QCDSR method reported earlier Cai17QCDSR , as we shall see, the present results on the nuclear symmetry energy obtained via the QCDSR method are also found to be consistent with predictions by other state-of-the-art microscopic many body theories, indicating that QCDSR can be applied to explore the EOS of ANM in a quantitative manner. The QCDSR method thus establishes a connection between the EOS of ANM and the non-perturbative QCD vacuum.

The paper is organized as follows. Section II gives a general description of the EOS of ANM and the Lorentz structure of nucleon self-energies in ANM. In Section III, the physical foundation of the QCDSR is introduced. In particular, the techniques on the calculation of the Wilson’s coefficients and the properties of the condensates used in this work are given. Section IV is devoted to the nucleon mass in vacuum within the QCDSR method, and the major task of this section is to determine the Ioffe parameter $t$ through the physical nucleon mass in vacuum, which is a starting point of the following investigation on the EOS of ANM. In Section V, the self-energy structure of the symmetry energy is analyzed via the simplified QCDSR, and the high order symmetry energy effects on the EOS of ANM are briefly discussed. Then Section VI studies these high order symmetry energy effects in some detail. Section IV, Section V, and Section VI are mainly qualitative, and are included to reveal some important features on the symmetry energy and nuclear matter EOS from QCDSR. Several important analytical expressions are given in these three sections. In Section VII, a full calculation on the EOS of SNM, the symmetry energy and the EOS of PNM is given, where the first parameter set QCDSR-1 (naive QCDSR) is constructed. In Section VIII, the effects of the higher order density terms in the chiral condensates on the symmetry energy and the PNM EOS are explored, and as a result, the second parameter set, i.e., QCDSR-2, is given. In Section IX, the contribution from the twist-four four-quark condensates to nucleonic matter EOS is studied with the third parameter set QCDSR-3 constructed. Section X gives the summaries and outlook of the work.

II EOS of ANM and Nucleon self-energies

In this section, we give a general description of the EOS of ANM and its relation to the in-medium nucleon self-energy structure in its relativistic form.

II.1 Definition of EOS of ANM

The EOS of ANM can be obtained through the total energy density $\varepsilon(\rho,\delta)$ by $E(\rho,\delta)={\varepsilon(\rho,\delta)}/{\rho}-M$ where $M=0.939\,$GeV is the nucleon rest mass. Moreover, the $E(\rho,\delta)$ can be expanded as a power series of even-order terms in $\delta$ as

[TABLE]

where $E_{0}(\rho)$ is the EOS of SNM, and the symmetry energy $E_{\mathrm{sym}}(\rho)$ is expressed as

[TABLE]

Around the saturation density $\rho_{0}$, the $E_{0}(\rho)$ can be expanded, e.g., up to 2nd-order in density, as,

[TABLE]

where $\chi=(\rho-\rho_{0})/3\rho_{0}$ is a dimensionless variable characterizing the deviations of the density from the saturation point $\rho_{0}$. The first term $E_{0}(\rho_{0})$ on the right-hand-side of Eq. (2.3) is the binding energy per nucleon in SNM at $\rho_{0}$ and the coefficient $K_{0}$ of the second term is

[TABLE]

which is the so-called incompressibility coefficient of SNM.

Similarly, one can expand the $E_{\mathrm{sym}}(\rho)$ around an arbitrary reference density $\rho_{\text{r}}$ as (see, e.g., ref. Zha13 )

[TABLE]

with $\chi_{\text{r}}=(\rho-\rho_{\text{{r}}})/3\rho_{\text{r}}$, and the slope parameter of the nuclear symmetry energy at $\rho_{\text{r}}$ is expressed as

[TABLE]

For $\rho_{\text{r}}=\rho_{0}$, the $L(\rho_{\text{r}})$ is reduced to the conventional slope parameter $L\equiv 3\rho_{0}{\text{d}E_{\mathrm{sym}}(\rho)}/{\text{d}\rho}|_{\rho=\rho_{0}}$.

II.2 Lorentz structure of nucleon self-energies

For the translational and rotational invariance, parity conservation, time-reversal invariance, and hermiticity in the rest frame of infinite nucleonic matter, the nucleon self-energy may be written generally in the relativistic case as Wal74 ; Mil74 ; Jam81 ; Hor83 ; Hor87 ; Ser86 ,

[TABLE]

where the isospin and density dependence of the nucleon self-energy are suppressed. The quantities $\Sigma_{\textrm{K}}(|\textbf{k}|,k^{0})$, $\Sigma_{\textrm{S}}(|\textbf{k}|,k^{0})$ and $\Sigma_{\textrm{V}}(|\textbf{k}|,k^{0})\equiv-\Sigma^{0}(|\textbf{k}|,k^{0})$ are Lorentz (rotational) scalar functions of $|\textbf{k}|$ and $k^{0}$ (the Minkovski metric is selected as $g_{\mu\nu}=(+,-,-,-)$ in the present work), $\textbf{k}^{0}=\textbf{k}/|\textbf{k}|$ is the unit vector along the direction of the momentum k. In the rest frame of infinite nucleonic matter, these invariants can be expressed in terms of $k^{0}$, $|\textbf{k}|$ (and $\rho$ as well as isospin $\delta$). The general proof of the decomposition Eq. (2.7) can be found in ref. Ser86 .

The effects of interactions between nucleons on the propagation of a nucleon in the medium can be included to all orders via Dyson’s equation Fet71 , i.e.,

[TABLE]

where $G^{0}(k)$ is the noninteracting nucleon Green’s function (propagator) and $\Sigma$ is the proper self-energy. Eq. (2.7) and Eq. (2.8) are completely general in the rest frame of infinite matter (in this case, the nucleon current-density four-vector has only a time-like non-vanishing component) and in principle could be used to determine $G$ exactly. Dyson’s equation Eq. (2.8) can be solved formally, yielding

[TABLE]

The location of the poles in $G(k)$ may be specified using the modified Feynman diagrams approach.

By defining the (Dirac) effective mass as well as the effective four-momentum of a nucleon, i.e.,

[TABLE]

one can rewrite the solution of Eq. (2.8) in a compact form, i.e., $G(k)=G_{\textrm{F}}(k)+G_{\textrm{D}}(k)$, where

[TABLE]

and

[TABLE]

$G_{\textrm{F}}$ and $G_{\textrm{D}}$ are two parts originated from the Pauli exclusion principle and the propagation of real nucleons in the Fermi sea in the interacting nucleonic matter, respectively Hor83 ; Hor87 .

The total single particle energy $e$ can be obtained from the dispersion relation, i.e., the solution of the following transcendental equation

[TABLE]

which evidently depends on $|\textbf{k}|$, the density $\rho$ and the energy $e$ itself. The above results are valid for any approximations to the self-energy in infinite matter. In order to arrive the Hartree–Fock approximation, for example, we include in $\Sigma$ only the contributions from tadpole and exchange diagrams Fet71 ; Hor83 ; Hor87 ; Ser86 . Moreover, if the self-energy has no explicit energy dependence, then one obtains

[TABLE]

When the above expression is generalized to ANM with any isospin asymmetry $\delta$, we then have

[TABLE]

where the isospin and density dependence of the quantity is recovered. Due to the general smallness of $\Sigma_{\textrm{K}}^{J}(\rho,\delta,|\textbf{k}|)$ Plo06 , we will neglect this term in the following study. Consequently, the single nucleon energy is given by

[TABLE]

In the present work, the EOS of ANM is obtained by the formulae based on the Hugenholtz–Van Hove (HVH) theorem Hug58 . More specifically Cai12 ,

[TABLE]

where “0” denotes SNM, i.e., $\Sigma_{\textrm{S(V)}}^{0/\textrm{n}}$ is the scalar (vector) self-energy in SNM/PNM, $e_{\textrm{F}}^{\ast}=({M_{0}^{\ast,2}+k_{\textrm{F}}^{2}})^{1/2}=[(M+\Sigma_{\textrm{S}}^{0})^{2}+k_{\textrm{F}}^{2}]^{1/2}$, $E_{\textrm{n}}$ is the EOS of PNM with $e_{\textrm{F,n}}^{\ast}=(M_{\textrm{n}}^{\ast,2}+k_{\textrm{F,n}}^{2})^{1/2}$, here $M_{\textrm{n}}^{\ast}$ is the neutron effective mass in PNM. Moreover, $k_{\textrm{F}}=\left({3\pi^{2}\rho}/{2}\right)^{1/3}$ ($k_{\textrm{F,n}}=2^{1/3}k_{\textrm{F}}$) is the Fermi momentum in SNM (PNM), $\Sigma_{\textrm{sym}}^{\textrm{S/V}}\equiv\Sigma_{\textrm{sym},1}^{\textrm{S/V}}$ is the first-order symmetry self-energy Cai12 . The main task of the present work is to explore the density/momentum dependence of the $\Sigma_{\textrm{S/V}}^{J}(\rho,\delta,|\textbf{k}|)$ via the QCDSR method, and obtain the EOS of SNM, the EOS of PNM and the symmetry energy through Eqs. (2.20), (2.21), and (2.22). Furthermore, the slope parameter of symmetry energy could be obtained by Eq. (2.6).

III Foundation of QCDSR

In this section, we briefly describe the physical foundation of QCDSR Coh91 ; Coh92 ; Fur92 ; Jin93 ; Jin94 ; Dru17a ; Dru17b ; Dru19 ; Dru04 ; Dru12 ; Jeo13 ; Coh95 . We first discuss the QCDSR in vacuum, which is relatively simple but contains all the important ingredients of the method. The generalization of the QCDSR in vacuum to finite densities is then followed. The quark/gluon condensates used in this work and the fitting scheme are finally given.

III.1 QCDSR in Vacuum

We start our discussions first by introducing the QCDSR in vacuum. In this work, $A_{\lambda}^{A}$ denotes the gluon field where $A=1\sim 8$ is the color index and $\lambda=0\sim 3$ is the space-time index Coh95 ; Pes95 ; Wei95 ; Zee11 ; Sre07 . The matrix form of the gluon field has the following structure,

[TABLE]

where $t^{A}=\lambda^{A}/2$ and $\lambda^{A}$’s are the Gell-Mann matrices which have the following basic properties

[TABLE]

with $f^{ABC}$ being the structure constant of the group SU(3). The strength tensor for gluon field is given by

[TABLE]

where $D_{\mu}=\partial_{\mu}-ig_{\textrm{s}}A_{\mu}$ is covariant derivative and $g_{\textrm{s}}$ is the coupling constant. Another form of the above equation is,

[TABLE]

and

[TABLE]

In order to discuss QCDSR in vacuum, one should introduce appropriate correlation functions of nucleons, here we adopt Coh95

[TABLE]

where $q$ is the momentum transfer between nucleons, $|0\rangle$ is the non-perturbative physical vacuum, $\eta_{\alpha}$ is the interpolation field of nucleons, and $\alpha,\beta$ are Dirac spinor index. The interpolation field for proton is given by Iof81 ,

[TABLE]

where u and d are quark fields and $a,b$ and $c$ are the color index. T represents transpose of the quark field in Dirac space, and C is the charge conjugation operator. The central quantity in QCDSR are the spectral functions (densities),

[TABLE]

Using the spectral functions $\rho_{\alpha\beta}(q)$ and $\widetilde{\rho}_{\alpha\beta}(q)$, we can rewrite the correlation function as

[TABLE]

with $q^{\prime}_{\mu}=({q}_{0}^{\prime},\textbf{q})$. In fact, the spectral functions can always be written in the following form after inserting a set of intermediate states,

[TABLE]

Using $\delta^{4}(q+P_{n})$ instead of $\delta^{4}(q-P_{n})$ gives a similar formula for $\widetilde{\rho}_{\alpha\beta}$ where $P_{n}^{\mu}$ is four momentum of state $n$.

Lorentz symmetry and parity invariance together mean that the general structure of $\rho_{\alpha\beta}$ is

[TABLE]

where $\rho_{\textrm{s}}$ and $\rho_{\textrm{q}}$ are scalar functions of $q$. Correspondingly, we have

[TABLE]

In the vacuum, we only need to contain integral for positive energy (which shall be modified in the finite density case), where the coefficients are Coh95

[TABLE]

with $s$ the threshold parameter ($\sim M^{2}$ for a nucleon). For example, the simplest phenomenological nucleon spectral densities take the form $\rho_{\textrm{s}}^{\textrm{phen}}(s)=M\delta(s-M^{2})$ and $\rho_{\textrm{q}}^{\textrm{phen}}(s)=\delta(s-M^{2})$, corresponding to $\Pi(q)=-(\not{q}+M)/(q^{2}-M^{2}+i0^{+})$, which is the standard nucleon propagator in vacuum, i.e., the two-point nucleon-nucleon correlation function.

The other important aspect of QCDSR is the OPE. For two local operators $A$ and $B$, we have Wil69 ; Pes95 ; Wei95

[TABLE]

where $C_{n}^{AB}$’s are the Wilson’s coefficients which can be calculated by standard perturbative methods, and $\mu$ is the renormalization energy scale. In the momentum space, we then have the correlation function from the OPE side as,

[TABLE]

where $Q^{2}=-q^{2}$ and $\langle\mathcal{O}_{n}\rangle$’s are the different types of quark/gluon condensates. Notice that OPE is only applicable in the large $Q^{2}$ region, i.e., in the deep space-like region.

Physically, it is no prior that the correlation functions from OPE side should be same as these from the phenomenological side, and they could also be very different from each other. The basic assumption of QCDSR is that in some range of $q^{2}$ these different correlation functions are physically equivalent. This range of $q^{2}$, or equivalently, the range of applicability, is called QCDSR window. At this point, it should be pointed out that the QCDSR approach is usually expected work well at lower densities/momenta. The nucleon spectral functions in nuclear medium are very complicated, and at low densities/momenta there exists a very narrow resonance state corresponding to the nucleon which can be described as a delta function. As density/momentum increases, continuum excitations emerge and these high density/momentum states will have increasing importance at high densities/momenta. However, in QCDSR, the contributions from high order states (e.g., high density/momentum momentum states) are significantly suppressed by the Borel transformation of correlation functions, leading to that QCDSR shall be mainly applicable in the low density/momentum region.

According to the above general analysis of spectral functions, one can write out the general structure of the nucleon spectral densities as

[TABLE]

Delta function indicates that it is a resonance (the nucleon) and the ellipsis denotes high order states. In the above expressions, $\lambda$ is the constant related to two physical states, $|0\rangle$ and $|q\rangle$ connected through $\langle 0|\eta(0)|q\rangle=\lambda u(q)$ with $q^{2}=M^{2}$ and $u(q)$ the Dirac spinor, $s$ is a threshold parameter. Correspondingly, the correlation functions are given by

[TABLE]

which shall be rewritten in a unified form,

[TABLE]

When the Borel transformation is made on the correlation functions from both the phenomenological side and the OPE side, one obtains the QCDSR equations, which connect the nucleon self-energies appearing on the phenomenological side and the quark/gluon condensates on the OPE side. Before giving the Borel transformation of the correlation functions, we discuss the essential procedures on QCDSR calculations:

1. Firstly, we determine the interpolation field to be studied, for example, the interpolation field for proton, Eq. (3.7), or the more general expression,

[TABLE]

with two independent terms,

[TABLE]

In Eq. (3.22), $t$ is a parameter whose natural value is around $-1$ Iof81 . The interpolation field for proton with $t=-1$ is called the Ioffe interpolation field. In Section IV, we determine the $t$ in a self-consistent manner instead of using $t=-1$. In order to obtain the interpolation field for the neutron, one just needs to exchange u and d in Eq. (3.22).

2. The second step is to determine the tensor structure of the spectral functions. For instance, there are only “s” and “q” parts for the nucleon correlation functions in vacuum while a new term will emerge at finite densities.

3. Then one writes down the dispersion relations for the nucleon correlation functions on the phenomenological side, e.g., Eq. (3.14), which is a fundamental step in QCDSR, and in the next subsection we give the general correlation functions for nucleon at finite densities.

4. At the same time, one writes down the OPE for the interpolation fields in terms of quark/gluon condensates, where the central quantities in this step are the Wilson’s coefficients. Using perturbative method from standard quantum field theories will furnish this calculation.

5. Finally, one makes the Borel transformation on two types of the correlation functions, i.e., one from the phenomenological side and the other from the OPE side, and then obtains the QCDSR equations. Solving the QCDSR equations and analyzing the results are the next procedure in the whole program.

At last of this subsection, we discuss the Borel transformation. For any function of the momentum transfer, $f(Q^{2})$, the Borel transformation

[TABLE]

is defined through

[TABLE]

where $\mathscr{M}$ is the Borel mass. For instance, Borel transformation of some typical functions are given by

[TABLE]

where $m$ equals to $0,1,2,\cdots$, $k$ equals to $1,2,3,\cdots$, and $\gamma_{\textrm{E}}\approx 0.577$ is the Euler constant.

Under Borel transformation, the correlation function (3.14) becomes

[TABLE]

where the polynomials in Eq. (3.14) disappear. The disappearance of the polynomials is the most important approximation made by the Borel transformation. As discussed above, the high order states, and/or, the continuum excitations including polynomials in Eq. (3.14), become more and more important at high densities/momenta, the disappearance of the polynomials under Borel transformation would make the physical predictions eventually incredible at high densities/momenta. This is the main reason why QCDSR should mainly be used in low density/momentum region. Moreover, the high-$s$ states become unimportant due to the suppression factor $e^{-s/\mathscr{M}^{2}}$, and they can be even removed (as the polynomials in Eq. (3.14)). As a rough example on the density region above which the QCDSR should be broken down, we consider the formation of the $\Delta$ resonance as an excited state in dense nucleonic matter. As shown in ref. Cai15a , the formation density of the first charged state of $\Delta(1232)$ could be smaller than 2$\rho_{0}$, even to be around the saturation density. Thus it is conservative to expect that the QCDSR should not be applied at densities around or larger than 2$\rho_{0}$. However, a comprehensive analysis of the applicable region of the QCDSR deserve more further work.

III.2 QCDSR at Finite Densities

In this subsection, we generalize the QCDSR in vacuum to finite densities. The nucleon propagator in medium at finite densities is given by Ser86 ; Hor83 ; Hor87 ; Rin96 ; Ser97 ,

[TABLE]

where $\Phi_{0}$ is the physical ground state for the infinite nucleonic matter and $\psi$ is the corresponding nucleon field. The nucleon self-energy $\Sigma(q)$ is defined through Dyson’s equation in the following form (with isospin index suppressed),

[TABLE]

which can be decomposed into

[TABLE]

by symmetry principles, where $u_{\mu}$ is the nucleon four-velocity and $qu=q_{\mu}u^{\mu}$. The “u” term, i.e., $G_{\textrm{u}}(q^{2},qu)\not{u}$, is new at finite densities. Similarly, we decompose the self-energy into the corresponding terms,

[TABLE]

with

[TABLE]

Defining the scalar self-energy in nuclear medium as $\Sigma_{\textrm{S}}=M^{\ast}-M$ with

[TABLE]

and the vector self-energy as

[TABLE]

we can then rewrite the propagator of a nucleon as

[TABLE]

where $\lambda^{\ast,2}$ is the residual factor Coh95 , and ${e}$ and $\overline{{e}}$ are the poles of positive energy branch and negative energy branch, i.e.,

[TABLE]

with

[TABLE]

The discontinuity passing through the real axis of $q$ represents the spectral function of the correlation function, i.e.,

[TABLE]

Based on the nucleon propagator given above, the nucleon correlation functions can be obtained correspondingly,

[TABLE]

where $\eta$ is proton’s interpolation field. Very similarly, we decompose $\Pi(q)$ into three parts,

[TABLE]

with

[TABLE]

Furthermore, from the discussions on the propagator above, we have

[TABLE]

and their Borel transformations are $\lambda^{\ast,2}M^{\ast}e^{-({e}^{2}-\textbf{q}^{2})/\mathscr{M}^{2}}$, $\lambda^{\ast,2}e^{-({e}^{2}-\textbf{q}^{2})/\mathscr{M}^{2}}$, and $\lambda^{\ast,2}\Sigma_{\textrm{V}}e^{-({e}^{2}-\textbf{q}^{2})/\mathscr{M}^{2}}$, respectively.

On the other hand, the correlation functions constructed from quark/gluon condensates are

[TABLE]

where $\langle\mathcal{O}_{n}\rangle_{\rho}=\langle\Psi_{0}|\mathcal{O}_{n}|\Psi_{0}\rangle$ are quark/gluon condensates at finite densities. In this work, the quark/gluon condensates at finite densities up to mass dimension-6 are included in the QCDSR equations, i.e., $\langle\overline{q}q\rangle$, $\left\langle({\alpha_{\textrm{s}}}/{\pi})G^{2}\right\rangle$, $\langle g_{\textrm{s}}\overline{q}\sigma\mathcal{G}q\rangle$, $\langle g_{\textrm{s}}{q}^{{\dagger}}\sigma\mathcal{G}q\rangle$, $\langle\overline{q}\Gamma_{1}q\overline{q}\Gamma_{2}q\rangle$ and $\langle\overline{q}\Gamma_{1}\lambda^{A}q\overline{q}\Gamma_{2}\lambda^{A}q\rangle$. Properties of them will be given in the following subsections.

III.3 OPE Coefficients

Wilson’s coefficients (OPE coefficients) $C_{n}^{i}(q^{2},qu)$ could be calculated by standard perturbative method Coh91 ; Coh92 ; Fur92 ; Jin93 ; Jin94 ; Jeo13 ; Coh95 . To make the discussions simpler, we calculate the OPE coefficients in SNM Coh95 based on Fock–Schwinger gauge in background field method Nov84 . The generalizations to the ANM are straightforward without difficulty, which will be given in the last of this subsection. The Fock–Schwinger gauge is

[TABLE]

with $\mathcal{A}^{\mu}\equiv A^{A\mu}t^{A}$. Eq. (3.55) indicates

[TABLE]

Moreover, with

[TABLE]

one then has

[TABLE]

Using $y^{\nu}=\alpha x^{\nu}$, we then have

[TABLE]

so

[TABLE]

then

[TABLE]

In the background field method, the non-perturbative effects of quarks are represented by Grassmann background fields, $\chi_{a\alpha}^{q},\overline{\chi}_{a\alpha}^{q}$, while the effects of gluons are represented by their classical fields, $F_{\mu\nu}^{A}$. The propagator of two quarks in coordinate space reads Nov84

[TABLE]

The products of Grassmann fields and classical fields can be written as products of matrix elements of the ground states of quarks and gluons, i.e.,

[TABLE]

etc., where all the values of fields are calculated at $x=0$. Then we can write the propagator of quarks as

[TABLE]

In Fig. 1, we show the quark propagator in nuclear medium graphically, where the last three terms represent non-perturbative effects.

The matrix element $\langle q_{a\alpha}(x)\overline{q}_{b\beta}(0)\rangle_{\rho}$ can be projected as,

[TABLE]

where parity symmetry and color neutrality of the ground state of nucleonic matter are taken into account when writing down the above expression. At short distance, we expand the quark field as

[TABLE]

Then using the Fock–Schwinger gauge, one obtains an identity for gluons,

[TABLE]

Each term above is zero, and we thus have

[TABLE]

and similarly, we have

[TABLE]

etc. After these simplifications on gluon fields, the quark fields then can be re-expanded as

[TABLE]

Finally, we obtain

[TABLE]

where all the derivatives are calculated at $x=0$.

The calculations of Wilson’s coefficients of condensates $\langle\overline{q}D_{\mu_{1}}D_{\mu_{2}}\cdots D_{\mu_{n}}q\rangle_{\rho}[\langle\overline{q}\gamma_{\mu}D_{\mu_{1}}D_{\mu_{2}}\cdots D_{\mu_{n}}q\rangle_{\rho}]$ are very similar to that of $\langle\overline{q}q\rangle_{\rho}[\langle\overline{q}\gamma_{\mu}q\rangle_{\rho}]$. For example, in the coordinate space, we have

[TABLE]

In the momentum space, we have

[TABLE]

The condensate terms in Eq. (III.3) can be furnished by writing out all the possible terms required by Lorentz symmetry, while the coefficients of these terms can be obtained through their traces of the decomposition, i.e.,

[TABLE]

where $g_{\mu\nu}=(+,-,-,-)$ and $u\cdot D=uD=u_{\mu}D^{\mu}$, $\sigma\cdot\mathcal{G}=\sigma\mathcal{G}$.

Equations of motion are useful when deriving these expressions, e.g., for the second identity in the second line, we have used the relation $D_{\mu}\equiv 2^{-1}(\gamma_{\mu}\not{D}+\not{D}\gamma_{\mu})$, combined with the translation invariance gives $\langle\overline{q}i\not{D}\not{u}q\rangle_{\rho}=-\langle\overline{q}i\overleftarrow{\not{D}}\not{u}q\rangle_{\rho}$. Very similarly, we have

[TABLE]

where the relation $\langle\overline{q}\not{u}D^{2}q\rangle_{\rho}=2^{-1}\langle g_{\textrm{s}}\overline{q}\not{u}\sigma\mathcal{G}q\rangle_{\rho}$ is used. Expressions for other condensates with different dimensions are similarly obtained, i.e.,

[TABLE]

where we have $(u\cdot G)^{2}=u^{\lambda}G_{\lambda\nu}^{A}u_{\mu}G^{A\mu\nu},\widetilde{G}^{A\mu\nu}=2^{-1}\varepsilon^{\mu\nu\kappa\lambda}G_{\kappa\lambda}^{A}$. Furthermore, the decomposition of the four-quark condensates has the following form,

[TABLE]

In order to translate the above expressions from coordinate space into momentum space, we use the following formulae,

[TABLE]

where $P_{m}(q^{2})$ is the polynomial of $q^{2}$ of order $m$. We decompose the correlation functions into their odd and even parts, i.e.,

[TABLE]

with

[TABLE]

where we have $q_{0}^{2}=(q\cdot u)^{2},u=(1,\textbf{0})$, and $c_{1}=7t^{2}-2t-5,c_{2}=1-t^{2},c_{3}=2t^{2}-2t-1,c_{4}=5t^{2}+2t+5$ and $c_{5}=7t^{2}+10t+7$ Jin94 . In this work, the following condensates are included for the QCDSR equations,

[TABLE]

III.4 QCDSR Equations in ANM

Now we generalize the above results of SNM to the case of ANM, and the QCDSR equations for the proton are given by Jeo13 ,

[TABLE]

where $\mathscr{M}$ is the Borel Mass and $\lambda_{\textrm{p}}^{\ast}$ is the residual for quasi-proton, and

[TABLE]

together with

[TABLE]

are the total effective single particle energy for proton (proton-hole). In order to write down the QCDSR equations for the neutron, one can exchange the quark fields of d and u. In the above equations, $\langle\cdots\rangle_{\rho,\delta}$ represents the condensates at finite density $\rho$ and isospin asymmetry $\delta$. $B_{\textrm{tw4}}^{\textrm{II}}$ and $B_{\textrm{tw4}}^{\textrm{III}}$ are contributions from twist-four four-quark condensates Jeo13 , which will be omitted until Section IX where we discuss their physical effects in some detail. The wave line in the four-quark condensates means Coh95

[TABLE]

with $f$ an effective parameter introduced Cai17QCDSR , and

[TABLE]

characterizes the anomalous dimension of the interpolation fields, $\Lambda_{\textrm{QCD}}\approx 0.17\,\textrm{GeV}$ is the QCD energy scale, and $\mu\approx 0.5\,\textrm{GeV}$ is the corresponding renormalization scale. The energy dependence of the QCD coupling constant, i.e., $\alpha_{\textrm{s}}\equiv g_{\textrm{s}}^{2}/4\pi$ is given by $\alpha_{\textrm{s}}(\textbf{q}^{2})={4\pi}/{9\ln(\textbf{q}^{2}/\mu^{2})}$. The factor $L$ contains radiation effects effectively Jeo13 . Furthermore, the contributions from continuum excitations whose physical origin will be discussed in Section IV, are included by the following functions,

[TABLE]

where we have $s_{0}^{\ast}=\omega_{0}^{2}-\textbf{q}^{2}$ with $\omega_{0}$ representing the effect of continuum excitations.

In Eq. (3.100), Eq. (3.101) and Eq. (3.102), the five parameters, i.e., $\mathscr{C}_{4},\mathscr{C}_{5},\mathscr{C}_{\textrm{H}},\mathscr{C}_{\textrm{B}},\mathscr{C}_{\textrm{D}}$ are introduced. If they take the value “+1” then the corresponding contributions are included, otherwise if they take the value “0” then the corresponding contributions are absent. More specifically, $\mathscr{C}_{4}$ characterizes the four-dimensional condensates, $\mathscr{C}_{5}$ for the five-dimensional condensates, $\mathscr{C}_{\textrm{H}}$ for the contributions from quasi-hole effects, while $\mathscr{C}_{\textrm{B}}$ and $\mathscr{C}_{\textrm{D}}$ characterize the four-quark condensates of the types of $\langle\overline{q}q\rangle_{\rho,\delta}^{2}$ and $\langle q^{{\dagger}}q\rangle_{\rho,\delta}^{2}$, respectively. In the following analysis, effects from the anomalous dimension and the continuum excitations denoted by the two parameters $\mathscr{C}_{\textrm{A}}$ and $\mathscr{C}_{\textrm{C}}$, will be also studied. If $\mathscr{C}_{\textrm{A}}=0$, then $L=1$ and if $\mathscr{C}_{\textrm{C}}=0$, the functions $E_{0},E_{1}$ and $E_{2}$ take the value 1.

III.5 Quark and Gluon Condensates Used in This Work

In this subsection, we discuss the properties of quark/gluon condensates with their density dependence, which are used as input for QCDSR equations. For the three-dimensional chiral condensate $\langle\overline{q}q\rangle_{\rho,\delta}$, one could introduce Jeo13 ,

[TABLE]

where $|\textrm{p}\rangle$ is the proton state. The mass of a nucleon in vacuum can be represented through the trace of energy-momentum tensor, i.e., $M\langle\textrm{N}|\overline{\psi}\psi|\textrm{N}\rangle=\langle\textrm{N}|\mathcal{T}|\textrm{N}\rangle$, where

[TABLE]

with $\mathcal{T}^{\mu\nu}$ the energy momentum tensor, $h$ denotes the heavy quark field, $\overline{\beta}=-9\alpha_{\textrm{s}}^{2}/2\pi$ is the reduced Gell-Mann–Low functions. For the baryon octet, we have

[TABLE]

with $A=\langle(\overline{\beta}/4\alpha_{\textrm{s}})G^{2}\rangle_{\textrm{p}},B_{\textrm{u}}=\langle\overline{\textrm{u}}\textrm{u}\rangle_{\textrm{p}}$, etc. After straightforward calculations, we obtain

[TABLE]

with $m_{\textrm{q}}=2^{-1}({m_{\textrm{u}}+m_{\textrm{d}}})$. On the other hand, we have

[TABLE]

in which the nucleon sigma term is $\sigma_{\textrm{N}}\equiv m_{\textrm{q}}{\mathrm{d}M}/{\mathrm{d}m_{\textrm{q}}}\approx 45\,\textrm{MeV}$ Gas91 (see also ref. Ala12 ). Introducing $\langle\textrm{p}|\overline{\textrm{u}}\textrm{u}|\textrm{p}\rangle\pm\langle\textrm{p}|\overline{\textrm{d}}\textrm{d}|\textrm{p}\rangle=a_{\pm}\langle\textrm{p}|\overline{\textrm{u}}\textrm{u}|\textrm{p}\rangle$, then we can rewrite the asymmetric part in terms of the symmetric part as

[TABLE]

where Jeo13

[TABLE]

The masses of the quarks and the baryons are given by PDG ,

[TABLE]

according to the particle data group (PDG), we then have

[TABLE]

In the following studies we denote $\alpha\equiv a_{-}$ and $\beta\equiv a_{+}$ and $\xi=\alpha/\beta\approx 0.1$. Collecting all the elements discussed above we finally obtain the chiral condensate at finite densities,

[TABLE]

where “$-$” is for the u quark and “+” for the d quark. The corresponding condensate in the vacuum takes the following value Jeo13 ; Coh95

[TABLE]

In Eq. (3.126), only the linear term in density $\rho$ is considered. In Section VIII, we will consider possible higher order terms in density and study how these higher order terms affect the nucleonic matter EOS.

Finally, we list other condensates of quarks and gluons. Another three-dimensional condensate is $\langle q^{{\dagger}}q\rangle_{\rho,\delta}^{\textrm{u,d}}$ Coh95 ,

[TABLE]

which is the actually the quark density Coh92 ; Fur92 ; Jin93 ; Jin94 ; Coh95 ; Jeo13 . Other quark/gluon condensates include Coh92 ; Fur92 ; Jin93 ; Jin94 ; Coh95 ; Jeo13 :

1. Four-dimensional condensates,

[TABLE]

where $\vartheta_{1}\approx 0.35,\varphi_{1}\approx 0.55$ Jeo13 . The lowest order of gluon condensates have mass dimension four, i.e.,

[TABLE]

where $G^{2}=2(\textbf{B}^{2}-\textbf{E}^{2})$ with B and E the magnetic field and electrical fields of QCD, respectively.

2. Five-dimensional condensates. In dimension five, we have several condensates constructed from quarks and gluons,

[TABLE]

where $\varphi_{2}\approx 0.34$, and,

[TABLE]

with $(620\,\textrm{MeV})^{2}\leq\langle g_{\textrm{s}}\overline{q}\sigma\mathcal{G}q\rangle_{\textrm{sym}}^{\textrm{p}}\leq(3\,\textrm{GeV})^{2}$, and

[TABLE]

with $\vartheta_{3}\approx 0.51,\varphi_{3}\approx 0.145$. Similarly, we have

[TABLE]

where $-(330\,\textrm{MeV})^{2}\leq\langle g_{\textrm{s}}{q}^{{\dagger}}\sigma\mathcal{G}q\rangle_{\textrm{sym}}^{\textrm{p}}\leq(660\,\textrm{MeV})^{2}$.

3. For the six-dimensional condensates, we consider the effective four-quark condensates defined in Eq. (3.105). And in Section IX, we study the twist-four four-quark condensates effects on the quantities we are interested in.

For more detailed physical discussions on the condensates, see, e.g., refs. Coh92 ; Fur92 ; Jin93 ; Jin94 and ref. Coh95 .

III.6 Fitting Scheme

In this work, the quark/gluon condensates at finite densities up to mass dimension-six are included in the QCDSR equations, see Eqs. (3.96), (3.97), (3.98), and (3.99), and the default central values are listed in the last subsection, see Eqs. (3.126) to (3.136). Moreover, the quark masses are taken to be $m_{\textrm{q}}=3.5\,\textrm{MeV}$, $m_{\textrm{s}}=95\,\textrm{MeV}$, the $\sigma$-N term $\sigma_{\textrm{N}}=45\,\textrm{MeV}$, the Borel mass $\mathscr{M}^{2}=1.05\,\textrm{GeV}^{2}$ Iof84 , and the threshold parameter defined in Eqs. (3.107), (3.108), and (3.109), as $\omega_{0}=1.5\,\textrm{GeV}$ Coh95 except a slight different $\omega_{0}=1.4\,\textrm{GeV}$ in Section V and Section VI.

Starting from Section VIII, we consider the following quark chiral condensates,

[TABLE]

The motivation for including the last term “$\Phi(1\mp g\delta)\rho^{2}$” in Eq. (3.137) is as follows Cai17QCDSR : As the density increases, the linear approximation for the chiral condensates becomes worse eventually, and higher order terms in density should be included in the $\langle\overline{q}q\rangle_{\rho,\delta}$. However, the density dependence of the chiral condensates is extremely complicated, and there is no general power counting scheme to incorporate these higher order density terms. Besides the $\rho^{2}$ term we adopted here, for instance, based on the chiral effective theories Kai09 ; Kru13a , a term proportional to $\rho^{5/3}$ was found in the perturbative expansion of $\langle\overline{q}q\rangle_{\rho,\delta}$ in $\rho$. On the other hand, using the chiral Ward identity God13 , a $\rho^{4/3}$ term was found in the density expansion in the chiral condensates. In our work, including the higher-order $\rho^{2}$ term is mainly for the improvement of describing the empirical EOS of PNM around and above saturation density, for which we use the celebrated Akmal–Pandharipande–Ravenhall (APR) EOS APR . In this sense, the $\Phi$-term adopted here is an effective correction to the chiral condensates beyond the linear leading-order. Two aspects related to the $\Phi$-term should be pointed out: 1). Without the higher-order $\rho^{2}$ term, the EOS of PNM around and above saturation density can not be adjusted to be consistent with that APR EOS, i.e., there exists systematic discrepancy between the QCDSR EOS and the APR EOS; 2). Using an effective correction with a different power in density, e.g., a $\rho^{5/3}$ term, the conclusion does not change, i.e., the EOS of PNM around and above saturation density can still be adjusted to fit the APR EOS, and the sign of the coefficients $\Phi$ and $g$ will not change although their absolute values change, see the results in Fig. 31. In the following, we abbreviate the QCDSR using the chiral condensate without the last term in Eq. (3.137) as the “naive QCDSR” Cai17QCDSR , which will be explored in detail in Subsection VII.4.

In carrying out the QCDSR calculations, we fix the central value of the $E_{\textrm{n}}(\rho)$ at a very low density $\rho_{\textrm{vl}}=0.02\,\textrm{fm}^{-3}$ to be consistent with the prediction by the ChPT Tew13 ; Kru13 , i.e., $E_{\textrm{n}}(\rho_{\textrm{vl}})\approx 4.2\,\textrm{MeV}$, the central value of the symmetry energy $E_{\textrm{sym}}(\rho)$ at a critical density $\rho_{\textrm{c}}=0.11\,\textrm{fm}^{-3}$ to be $E_{\textrm{sym}}(\rho_{\textrm{c}})\approx 26.65$ MeV Zha13 , and fit the EOS of PNM be close to the APR EOS as much as possible, via varying $\Phi$, $g$ and $f$. We note that the parameter $f$ defined in Eq. (3.105) is essentially determined by $E_{\textrm{n}}(\rho_{\textrm{vl}})$, and the overall fitting of the EOS of PNM to the APR EOS and the symmetry energy at $\rho_{\textrm{c}}$ determines the other two parameters $\Phi$ and $g$ Cai17QCDSR .

IV Nucleon Mass in Vacuum

In this section, we use the QCDSR to study the nucleon mass in vacuum. The motivation is twofold: firstly the connection between the nucleon mass and the quark/gluon condensates in QCDSR is explored, and most importantly the scheme for determining the Ioffe parameter $t$ Iof81 is given here, i.e., via the relation $M_{\textrm{QCDSR}}^{\textrm{vac,static}}(t)=939\,\textrm{MeV}$.

Conventionally, the high energy states or the continuum states are not given by the QCDSR method itself, and in order to model these states in the nucleon spectral densities (3.17) and (3.18), one usually adopts the corresponding results from the OPE calculations Coh95 . More specifically, the nucleon correlation functions from the OPE side including only the lowest order terms are given by

[TABLE]

where the last term in $\Pi_{\textrm{q}}^{\textrm{OPE}}(q^{2})$ is the four-quark condensates expressed in terms of the square of the chiral condensate. Moreover, $t=-1$ is adopted here for simplicity. The Borel transformations of them are given by

[TABLE]

Consequently, the high energy/continuum states are approximated by the equivalent OPE terms, starting at the sharp threshold $s_{0}=\omega_{0}^{2}$, i.e. Coh95 ,

[TABLE]

see Fig. 2 for the sketch of the effects of $s_{0}$.

After putting these spectral functions into Eq. (3.14) and finishing the Borel transformation, we obtain the correlation functions at the phenomenological side in QCDSR,

[TABLE]

Then according to the QCDSR, i.e.,

[TABLE]

where $\widehat{\Pi}_{\textrm{s}}^{\textrm{OPE}}(\mathscr{M}^{2})$ and $\widehat{\Pi}_{\textrm{q}}^{\textrm{OPE}}(\mathscr{M}^{2})$ are given by Eq. (4.3) and Eq. (4.4), respectively, one obtains the expression for nucleon mass in vacuum as Footnotexx ,

[TABLE]

where $E_{0},E_{1}$ and $E_{2}$ are defined in Eq. (3.107), Eq. (3.108) and Eq. (3.109), and the following abbreviations are introduced,

[TABLE]

Furthermore, if we neglect all the high-dimensional except the three-dimensional chiral condensates, a very simple formula for the nucleon mass in vacuum is obtained, i.e.,

[TABLE]

demonstrating that the nucleon mass in vacuum is roughly determined by the the chiral condensates in vacuum.

Moreover, the above discussion clearly shows that the nucleon mass in vacuum obtained in QCDSR is not necessarily about its realistic value 939 MeV. For instance, taking $\mathscr{M}^{2}=1.05\,\textrm{GeV}^{2}$, $\langle\overline{q}q\rangle_{\textrm{vac}}=-(0.252\,\textrm{GeV})^{3}$, $\langle(\alpha_{\textrm{s}}/\pi)G^{2}\rangle_{\textrm{vac}}=(0.33\,\textrm{GeV})^{4}$, and $s_{0}=\omega_{0}^{2}=2.25\,\textrm{GeV}^{2}$, one obtains the mass $M\approx 0.847\,\textrm{GeV}$, which has about a 92 MeV deviation from 939 MeV (roughly an effect about 10%), and in Fig. 3, the nucleon mass in vacuum as a function of the Borel mass (squared) obtained by Eq. (4.11) is shown. In this work, we put an extra constraint on the QCDSR method, i.e., the nucleon mass in vacuum is set to be 939 MeV. Consequently the Ioffe parameter $t$ should be determined uniquely once the other parameters are fixed. Specifically, when we rewrite the QCDSR equations in vacuum under the static condition, i.e., taking $\rho=0$ and $\textbf{q}=\textbf{0}$, one obtains,

[TABLE]

where $c_{1}=7t^{2}-2t-5$ and $c_{4}=5t^{2}+2t+5$. Then the Ioffe parameter $t$ could be expressed as

[TABLE]

with

[TABLE]

In Fig. 4, the Ioffe parameter $t$ as a function of the Borel mass squared $\mathscr{M}^{2}$ with different threshold parameter $s_{0}$ is shown. The Borel mass squared was extensively studied in the literature and the value of it was found to fall within the range $0.8\,\textrm{GeV}^{2}\lesssim\mathscr{M}^{2}\lesssim 1.4\,\textrm{GeV}^{2}$ Iof84 ; Coh95 . In this work, we choose $\mathscr{M}^{2}\approx 1.05\,\textrm{GeV}^{2}$ as our default value for the Borel mass squared, and the Ioffe parameter obtained via Eq. (IV) is found to be around $-1.2$ considering $2\,\textrm{GeV}^{2}\lesssim s_{0}\lesssim 3\,\textrm{GeV}^{2}$, which is close to the Ioffe value ($t_{\textrm{Ioffe}}=-1$). For example, taking $s_{0}=3\,\textrm{GeV}^{2}$, the Ioffe parameter at $\mathscr{M}^{2}=1.05\,\textrm{GeV}^{2}$ is found to be about $t\approx-1.04$. Moreover, the Ioffe parameter will be set to its default value ($t_{\textrm{Ioffe}}=-1$) in Section V and Section VI mainly for qualitative discussions, and starting from Section VII, $t$ would be self-consistently determined by the scheme given here, i.e., Eq. (IV). The treatment on the Ioffe parameter given here is natural, thus one can obtain the physical nucleon mass in vacuum when the QCDSR equations are applied at zero density.

V Main Structure of The Symmetry Energy

In this section, we discuss the structure of the symmetry energy through the simplified QCDSR equations, which will be given shortly. The main purpose of this section is to qualitatively analyze the nucleon self-energy structure of the symmetry energy through Eq. (2.22), and leave the detailed/quantitative calculation on the $E_{\textrm{sym}}(\rho)$ to the following sections.

In the simplified version, the QCDSR equations for proton have the following form,

[TABLE]

i.e., they only include the three-dimensional chiral condensates, and all the high-dimensional condensates as well as the contributions from continuum excitations, quasi-hole effects are neglected. In the above expressions, $M_{\textrm{p}}^{\ast}$ is the proton’s Dirac effective mass at finite densities introduced through the following dispersion relation,

[TABLE]

The discussion in this section on the symmetry energy from QCDSR is more qualitative and thus the Ioffe parameter $t$ is set at its natural value, i.e., $t=-1$, as mentioned in the last section. Via exchanging the d quark and the u quark, one can obtain the corresponding sum rule equations for neutron, i.e.,

[TABLE]

and the total single particle energy for neutron is defined similarly as Eq. (5.4). The QCDSR with Eqs. (5.1), (5.2), (5.3), (5.5), (5.6), and (5.7) and the dispersion relation for proton/neutron, is called the mostly simplified QCDSR (“msQCDSR”). It is necessary to point out that in “msQCDSR”, the nucleon self-energies have no momentum dependence, thus the corresponding momentum-dependent terms in Eq. (2.22) are absent naturally. The momentum-dependent terms and their consequences will be explored in some detail in Section VII, and the corresponding contributions to the symmetry energy are found to be small (see the blue and green lines of Fig. 17).

The nucleon Dirac effective mass $M_{J}^{\ast}$ and the vector self-energy $\Sigma_{\textrm{V}}^{J}$ could be obtained immediately in msQCDSR as,

[TABLE]

The nucleon Dirac effective mass in SNM is then given by

[TABLE]

and the vector-self energy by

[TABLE]

with $\langle q^{{\dagger}}q\rangle_{\rho}=3\rho/2$.

Similarly, the nucleon Dirac effective mass and the vector self-energy in ANM are given by

[TABLE]

and

[TABLE]

respectively. The definitions of $\alpha$ and $\beta$ are given in Eq. (3.125). Consequently, the first-order symmetry scalar/vector self-energy is obtained as

[TABLE]

where the symmetry self-energy is generally defined as Cai12

[TABLE]

In the above definition, we use $\tau_{3}^{\textrm{n}}=+1$ and $\tau_{3}^{\textrm{p}}=-1$.

Through Eq. (2.22), one can obtain the symmetry energy in msQCDSR as

[TABLE]

with its decomposition given by

[TABLE]

where $e_{\textrm{F}}^{\ast}=(k_{\textrm{F}}^{2}+M_{0}^{\ast,2})^{1/2}$, and the sum of the last two terms is the potential part of the symmetry energy in msQCDSR $E_{\textrm{sym}}^{\textrm{pot}}(\rho)=E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)+E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)$. Several features of these expressions should be pointed out:

1. At low densities, the leading order term in the symmetry energy is the kinetic part, which roughly scales as $\rho^{2/3}$ since at low densities $e_{\textrm{F}}^{\ast}\approx M$. The potential part of the symmetry energy contributes starting from the linear terms in density $\rho$. Moreover, both the scalar and vector self-energy contributions to the symmetry energy are positive in msQCDSR, i.e., $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)>0$ and $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)>0$, and this finding will be verified in the full QCDSR.

2. Perhaps much more important is that, through Eq. (5.21) and Eq. (5.22), one establishes the connection between the symmetry energy and some other fundamental quantities, such as the quark mass $m_{\textrm{q}}$ and the nucleon-sigma term $\sigma_{\textrm{N}}$. For instance, a larger $m_{\textrm{q}}$ corresponds to a smaller $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ (dependence of $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ on $\sigma_{\textrm{N}}$ is non-trivial since $\sigma_{\textrm{N}}$ also affects the ratio $\alpha/\beta$, and will be discussed in detail in Section VII, see the fourth panel of Fig. 20).

3. Tracing back to the physical origin of the isospin factor $\alpha/\beta$, see Eq. (3.126), one can simply obtain the $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ if other approaches could give the density/isospin dependence of the chiral condensates. For instance, ref. Lac10 gave a different isospin effect of the chiral condensate, i.e., $\alpha/\beta\leftrightarrow c/f_{\pi}^{2}$ with $c$ a low energy coefficient and $f_{\pi}$ the pion decay constant, then

[TABLE]

where the last relation is obtained via the Gell-Mann–Oakes–Renner relation $m_{\pi}^{2}f_{\pi}^{2}=-m_{\textrm{q}}\langle\overline{q}q\rangle_{\textrm{vac}}$. The above relation shows the $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ scales as the square of the ratio $m_{\pi}/m_{\textrm{q}}$. It is fair to say that the $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ is the most non-trivial term in the decomposition of the symmetry energy and its density dependence largely characterizes the change of the density behavior of the total symmetry energy. More discussions on this issue will be given in the following sections.

4. On the other hand, the density dependence of $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)$ is much simpler and even when other contributions are included in the QCDSR equations, it only changes slightly.

An important problem related to the structure of the symmetry energy is the parabolic approximation of the EOS of ANM. Let us briefly discuss it in msQCDSR and let the more detailed investigations to the next section. Parabolic approximation of EOS of ANM could be conventionally characterized as

[TABLE]

or equivalently the parabolic approximation of the symmetry energy,

[TABLE]

i.e., the $E_{\textrm{sym}}^{\textrm{para}}(\rho)$ is the difference between the EOS of PNM and that of SNM. Consequently, the high order effect of the EOS of ANM is given by,

[TABLE]

or,

[TABLE]

since $E_{\textrm{n}}(\rho)=E_{0}(\rho)+E_{\textrm{sym}}(\rho)+E_{\textrm{sym,4}}(\rho)+E_{\textrm{sym,6}}(\rho)+\cdots$, where $E_{\textrm{sym,4}}(\rho)$, $E_{\textrm{sym,6}}(\rho)$, $\cdots$ are the fourth-order symmetry energy, sixth-order symmetry energy, $\cdots$  Cai12xx . In Fig. 5, the high order EOS $E_{\textrm{HO}}(\rho)$ as a function of density in msQCDSR is shown. At the saturation density, the high order effect is found to be about $E_{\textrm{HO}}(\rho_{0})\approx-2.4\,\textrm{MeV}$. Although many investigations on the high order symmetry energy both from microscopic calculations and phenomenological models LCK08 ; Cai12xx ; Ste06 ; Kai15 ; Nan16 ; Agr17 ; PuJ17 ; ZhZ17 ; RWang17 ; ZWZhang19 indicates that the magnitude of the high order term, especially the fourth-order symmetry energy, is small, e.g., roughly $|E_{\textrm{HO}}(\rho\lesssim\rho_{0})|\lesssim 1\,\textrm{MeV}$ at $\rho_{0}\approx 0.16\,\textrm{fm}^{-3}$, there are still no fundamental symmetries and/or principles guaranteeing its smallness. Thus it is an interesting issue to explore the high order EOS of the ANM in the framework of QCDSR, and this is the main subject of the next section.

Finally, let us qualitatively discuss the possible relation between the symmetry energy and the chiral condensates. In order to make the discussion general, we assume that the chiral condensates has the following density structure,

[TABLE]

with $\theta>1$, and we also assume that

[TABLE]

then we have

[TABLE]

Neglecting the high order term, i.e., the $\theta$ term temporarily, we have

[TABLE]

and consequently the first order symmetry scalar self-energy is given by

[TABLE]

and the corresponding contribution to the symmetry energy is

[TABLE]

Due to the empirical knowledge on the nucleon Dirac effective mass in SNM, the scalar self-energy part of the symmetry energy strongly depends on the structure of the u/d quark condensates, especially the isospin part, i.e., the properties of function $\mathscr{F}_{\textrm{u/d}}^{1}$. The $\mathscr{D}^{1}$ is given by $\mathscr{D}^{1}={\sigma_{\textrm{N}}}/{2m_{\textrm{q}}}$, thus

[TABLE]

If the u quark in PNM restores its chiral symmetry first, i.e., $\langle\overline{\textrm{u}}\textrm{u}\rangle_{\rho}$ approaches to zero earlier than $\langle\overline{\textrm{d}}\textrm{d}\rangle_{\rho}$ as the density increases, then $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ is positive; on the opposite side, $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ is negative. In the msQCDSR, one has $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)>0$, i.e., the u quark in PNM restores its chiral symmetry earlier than d quark. When we consider higher order terms in density (characterized by $\mathscr{F}^{\theta}_{\textrm{u/d}}$) in the chiral condensates, different patterns may emerge. It indicates that one can constrain the density dependence of the chiral condensates via the empirical knowledge of the symmetry energy within a reasonable density region. This is one of the motivations to include the higher order density terms in the chiral condensates as in Eq. (3.137), which will be further explored in detail in Section VIII.

VI More Discussions on $E_{\textmd{HO}}(\rho)$

In this section, we study the effects of the high order effects of the EOS of ANM. In order to study the parabolic approximation of the EOS of ANM, we add terms into the QCDSR equations eventually. For example, based on the msQCDSR, we add the four-quark condensates of the form $\langle\overline{q}q\rangle^{2}_{\rho,\delta}$, i.e.,

[TABLE]

into the QCDSR equation, and call the corresponding QCDSR the simple QCDSR (abbreviated as sQCDSR), which is denoted by $\mathscr{C}_{\textrm{B}}=1$. Specifically, the relevant QCDSR equations are modified as (where $t=-1$ is adopted in Eq. (6.2) and Eq. (6.3))

[TABLE]

compared with Eq. (5.2) and Eq. (5.6) for the msQCDSR, respectively. Similarly, $\mathscr{C}_{4}=1$ means four-dimensional condensates are included based on $\mathscr{C}_{\textrm{B}}=\mathscr{C}_{\textrm{A}}=1$, etc. The whole order of adding different types of effects is as follows: “msQCDSR” (only three-dimensional condensates) $\to$ $\mathscr{C}_{\textrm{B}}=1$ (sQCDSR with $B_{4}^{q}(f)$ included) $\to$ $\mathscr{C}_{\textrm{A}}=1$ (anomalous effects $L^{-2\Gamma_{\eta}+\Gamma_{\mathcal{O}_{n}}}$) $\to$ $\mathscr{C}_{\textrm{4}}=1$ (four-dimensional condensates) $\to$ $\mathscr{C}_{\textrm{5}}=1$ (five-dimensional condensates) $\to$ $\mathscr{C}_{\textrm{H}}=1$ (quasi-hole effects) $\to$ $\mathscr{C}_{\textrm{D}}=1$ (four-quark condensates of type $\langle q^{{\dagger}}q\rangle_{\rho,\delta}^{2}$, which is absent in vacuum) $\to$ “full QCDSR (fQCDSR)” ($\mathscr{C}_{\textrm{C}}=1$).

In these calculations, the Ioffe parameter $t=-1$ is adopted. In Fig. 6, the symmetry energy and its parabolic approximation are shown as functions of density in the above order. It is obvious from the figure that the parabolic approximation well behaves until the continuum excitations ($\mathscr{C}_{\textrm{C}}=1$), characterized by the functions $E_{0}$, $E_{1}$ and $E_{2}$, i.e., Eqs. (3.107), (3.108), and (3.109), are included. For example, in the full QCDSR (abbreviated as fQCDSR), the $E_{\textrm{HO}}$ at $\rho_{0}$ is found to be about $E_{\textrm{HO}}(\rho_{0})\approx 5.7\,\textrm{MeV}$. However, the $E_{\textrm{HO}}(\rho_{0})$ in the QCDSR to order $\mathscr{C}_{\textrm{D}}=1$ is about 0.7 MeV.

The investigations on the the dependence of the high order effect $E_{\textrm{HO}}(\rho)$ on the intrinsic parameters of the QCDSR will be useful for understanding the behavior of the parabolic approximation of the ANM EOS. In the following, we study the dependence of $E_{\textrm{HO}}(\rho)$ on the continuum excitation parameter $\omega_{0}$, the Borel mass $\mathscr{M}$, and the effective four-quark condensates parameter $f$. In Fig. 7, we show the $E_{\textrm{HO}}(\rho)$ and $E_{\textrm{sym}}(\rho)$ at (the cross density) $\rho_{\textrm{c}}=0.11\,\textrm{fm}^{-3}$ Zha13 ; Che12a ; Li12 ; Dan14 as functions of the continuum excitations factor $\omega_{0}$ (upper) and Borel mass squared $\mathscr{M}^{2}$ (lower) using the fQCDSR equations. Since the symmetry energy $E_{\textrm{sym}}(\rho)$ at the cross density $\rho_{\textrm{c}}$ is better constrained Zha13 ; Dan14 than that at the saturation density, we will also focus on the symmetry energy at $\rho_{\textrm{c}}$ besides its value at $\rho_{0}$ in the following sections. Two main features of Fig. 7 are necessary to be pointed out: 1). With the value of $\omega_{\textrm{0}}$ increasing, i.e., as the effects of the continuum excitations reduces, the parabolic approximation of the EOS of ANM becomes better, and this finding is consistent with the results shown in panels (g) and (h) of Fig. 6 (if $\omega_{0}\to\infty$, there will be no continuum excitation effects). 2). As the Borel mass (squared) increases, the parabolic approximation becomes better, however, the symmetry energy also becomes larger. The symmetry energy and the $E_{\textrm{HO}}$ display almost opposite variation tendency with the $\omega_{0}$ and $\mathscr{M}^{2}$. For example, the symmetry energy at $\rho_{\textrm{c}}$ with $\mathscr{M}^{2}=1.05\,\textrm{GeV}^{2}$ is found to be about 43.1 MeV, which is much larger than the empirical value about 26.6 MeV Zha13 , however, the corresponding $E_{\textrm{HO}}(\rho_{\textrm{c}})$ is found to be about 5.1 MeV. These results show that the parabolic approximation of the EOS of ANM in QCDSR is heavily broken, indicating the fourth-order symmetry energy even the higher order terms maybe large.

Now let us study the four-quark effective parameter $f$ on the EOS qualitatively. The effects of $f$ can be demonstrated semi-analytically in the sQCDSR (including the contributions from the four-quark condensates of the type $\langle\overline{q}q\rangle_{\rho,\delta}^{2}$ based on the msQCDSR). and then according to Eq. (2.20) and Eq. (2.21), one can obtain the $E_{0}(\rho)$ and $E_{\textrm{n}}(\rho)$, respectively. The $\langle\overline{q}q\rangle_{\rho,\delta}^{2}$ type of the four-quark condensates in SNM can be approximated as

[TABLE]

where $B_{4}$ (without superscript) denotes the corresponding term in SNM, $a=\sigma_{\textrm{N}}/2m_{\textrm{q}},A=\langle\overline{q}q\rangle_{\textrm{vac}}^{2},B=2a\langle\overline{q}q\rangle_{\textrm{vac}}$. In the above approximation, the term proportional to $\rho^{2}$ is omitted, since we are only interested in low density behavior of the EOS. The nucleon effective mass is given by (Ioffe parameter $t=-1$ in this estimate)

[TABLE]

where $\langle\overline{q}q\rangle_{\rho}$ is the chiral condensate in SNM, and

[TABLE]

Thus, $M_{0}^{\ast}(\rho)\approx{I}_{1}+{I}_{2}\rho+\mathcal{O}(\rho^{2})$, with

[TABLE]

In a very similar manner, one can obtain the approximation for the vector self-energy

[TABLE]

using the expression for $\langle q^{{\dagger}}q\rangle_{\rho}$, then $\Sigma_{\textrm{V}}^{0}(\rho)\approx{S}_{2}\rho+\mathcal{O}(\rho^{2})$ with ${S}_{2}={\mathscr{M}^{4}}/{\pi^{2}D}$. Then according to Eq. (2.20), one obtains

[TABLE]

to order $\rho$. The first term on the right hand side should be zero since $E_{0}(0)=0$, leading to $M=-{\mathscr{M}^{4}\langle\overline{q}q\rangle_{\textrm{vac}}}/{4\pi^{2}D}$ Footnote3 , and one obtains the following approximation for the EOS of SNM,

[TABLE]

where $E_{0}^{\textrm{FFG}}(\rho)=3k_{\textrm{F}}^{2}/10M$ is the free Fermi gas (FFG) prediction on the EOS of SNM, using the expressions for $I_{2}$ and $S_{2}$, the EOS of SNM can finally be written as

[TABLE]

Similarly, the nucleon effective mass in PNM could be approximated as $M_{\textrm{n}}^{\ast}(\rho)\approx{I}_{1}+W_{2}\rho$, with ${W}_{2}=\left(1-{\alpha}/{\beta}\right){I}_{2}$, and $\Sigma_{\textrm{V}}^{\textrm{n}}(\rho)\approx 5S_{2}\rho/4+\mathcal{O}(\rho^{2})$, then the EOS of PNM is given approximately by

[TABLE]

where $E_{\textrm{n}}^{\textrm{FFG}}(\rho)=3k_{\textrm{F,n}}^{2}/10M\sim\rho^{2/3}$ is the FFG prediction on the EOS of PNM. This expression is already very interesting. For example, since the EOS of PNM at very low densities (say densities smaller than $0.01\,\textrm{fm}^{-3}$) could be determined very accurate by simulations or microscopic calculations, there exists a relation between several fundamental quantities and the model parameters in QCDSR, such as the four-quark effective parameter $f$, the nucleon-sigma term $\sigma_{\textrm{N}}$ (the factor $\alpha/\beta$ also depends on $\sigma_{\textrm{N}}$), the chiral condensate in vacuum, the light quark mass, and the Borel mass. Specifically, a positive $f$ parameter leads to a reduction on the $E_{\textrm{n}}(\rho)$, which will be verified numerically in the following sections (e.g., see Fig. 28). Moreover, as discussed in Subsection III.6, the parameter $f$ would be essentially determined by the EOS of PNM at a very low density $\rho_{\textrm{vl}}\approx 0.02\,\textrm{fm}^{-3}$, and it is also indicated in Eq. (6.13). Based on the approximations for $E_{0}(\rho)$ and $E_{\textrm{n}}(\rho)$, the symmetry energy obtained in the parabolic approximation is roughly given by

[TABLE]

where the kinetic symmetry energy in the parabolic approximation is given by

[TABLE]

On the other hand, when using the exact nucleon self-energy decomposition formula Eq. (2.22) for calculating the symmetry energy, one obtains $E_{\textrm{sym}}^{\textrm{kin}}(\rho)=(3/5)E_{0}^{\textrm{FFG}}(\rho)$, $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)\approx-2^{-1}(\alpha/\beta)\rho I_{2}$, and $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)\approx 8^{-1}\rho S_{2}$, in the same approximation level (i.e., in sQCDSR). Interestingly, the potential part of the symmetry energy in the two approaches is found to be same at this order (i.e., order $\rho$), and only the kinetic part introduces the corresponding high order effects, which is about $-0.01E_{0}^{\textrm{FFG}}(\rho)$. Going beyond the linear order in density, the $f$ parameter will come into play in the high order effects of the EOS $E_{\textrm{HO}}(\rho)$. Moreover, the above qualitative analysis (approximation) on the $E_{\textrm{n}}(\rho)$, $E_{0}(\rho)$, and $E_{\textrm{sym}}(\rho)$ is useful for further investigations on, e.g., the correlation between the role played by the four-quark condensates and the EOS of ANM, as mentioned just above. Furthermore, if one takes $f=0$ in Eq. (6.13), then the $E_{\textrm{n}}(\rho)$ could be written as Cai17QCDSR

[TABLE]

which depends only on several fundamental quantities, such as $m_{\textrm{q}}$, $\langle\overline{q}q\rangle_{\textrm{vac}}$, and $\sigma_{\textrm{N}}$, and not on the effective parameters $f$ and $\mathscr{M}^{2}$. Despite its simplicity, Eq. (6.16) already has the power of quantitative predictions at very low densities. We will discuss more on this point in the following sections.

In Fig. 8, the effects of the parameter $f$ on the EOS of SNM and that of PNM in the fQCDSR are shown. Although the predictions on $E_{0}(\rho)$ and $E_{\textrm{n}}(\rho)$ are quantitatively incorrect compared to their empirical constraints, the effects of $f$ are obvious, which are consistent with the estimate given in Eq. (VI) and Eq. (6.13), i.e., as $f$ increases, the $E_{0}(\rho)$ and $E_{\textrm{n}}(\rho)$ are both enhanced correspondingly. However, the enhancement on $E_{\textrm{n}}(\rho)$ due to $f$ is less than that on $E_{0}(\rho)$ as shown from Eq. (VI) and Eq. (6.13) by a factor $\alpha/\beta$, leading to the result that as $f$ increases, the the symmetry energy is reduced, see Eq. (VI). The numerical results on the $E_{\textrm{sym}}(\rho)$ and $E_{\textrm{HO}}(\rho)$ using fQCDSR equations are shown in Fig. 9. For example, the high order effects in the EOS at $\rho_{\textrm{c}}=0.11\,\textrm{fm}^{-3}$ with $f=0$, $0.25$ and $0.5$ are found to be about 7.5 MeV, 6.6 MeV, and 6.3 MeV, respectively, and those values at $\rho_{0}=0.16\,\textrm{fm}^{-3}$ are found to be about 8.8 MeV, 7.4 MeV, and 6.4 MeV, respectively. These results again, when combining the features shown in Fig. 7 and that the Borel mass squared was constrained to fall within about $0.8\,\textrm{GeV}^{2}\lesssim\mathscr{M}^{2}\lesssim 1.4\,\textrm{GeV}^{2}$ Iof84 together with the continuum excitation factor about $2\,\textrm{GeV}^{2}\lesssim s_{0}^{\ast}\sim\omega_{0}^{2}\lesssim 3\,\textrm{GeV}^{2}$ Coh95 , indicate that the high order term $E_{\textrm{HO}}(\rho)$ at the cross density $\rho_{\textrm{c}}$ and at the saturation density $\rho_{0}$ are generally not small (e.g., $\lesssim 1\,\textrm{MeV}$).

Finally, it is useful to generally analyze the origin of the possible breakdown of the parabolic approximation of the EOS of ANM in QCDSR framework. In the msQCDSR, the density dependence of the chiral condensates has the following structure,

[TABLE]

see Eq. (3.126), where $a=\sigma_{\textrm{N}}/2m_{\textrm{q}}$ and $\xi=\alpha/\beta$. Moreover, the nucleon Dirac effective mass and the vector self-energy are given as $M_{J}^{\ast}\sim\langle\overline{q}q\rangle_{\rho,\delta}^{\textrm{u,d}}$ and $\Sigma_{\textrm{V}}^{J}\sim\langle q^{{\dagger}}q\rangle_{\rho,\delta}^{\textrm{u,d}}$, see Eq. (5.8) and Eq. (5.9), respectively. Then, the corresponding contribution originated from the scalar self-energy to the EOS of SNM via Eq. (2.20) can be obtained from the following estimate,

[TABLE]

or, roughly $\sim a\rho/2$, where “$\cdots$” in the above expression denotes relativistic corrections. Moreover, the relevant contribution to EOS of PNM is roughly given by $a(1-\xi)\rho/2$ via Eq. (2.21). Consequently, the symmetry energy obtained using the parabolic approximation, i.e., $E_{\textrm{sym}}^{\textrm{para}}(\rho)$, is roughly $-a\xi\rho/2$. It can be easily shown that the contribution of the vector self-energy to the symmetry energy exactly cancels in the parabolic approximation and in the self-energy decomposition for the $E_{\textrm{sym}}(\rho)$ in msQCDSR, both are $4\pi^{2}\rho/\mathscr{M}^{2}$, see Eq. (5.22). Consequently, the scalar self-energy contribution to the symmetry energy via the decomposition is roughly given by $-({M_{0}^{\ast}}/{{e}_{\textrm{F}}^{\ast}})a\xi\rho/2$ via simply subtracting $a(1-\xi)\rho/2$ and $a\rho/2$ and multiplying a factor $M_{0}^{\ast}/e_{\textrm{F}}^{\ast}$. Thus $E_{\textrm{HO}}(\rho)=E_{\textrm{sym}}^{\textrm{para}}(\rho)-E_{\textrm{sym}}(\rho)$ is found to be

[TABLE]

where the overall factor is unimportant for the qualitative demonstration. At low densities, $e_{\textrm{F}}^{\ast}=(M_{0}^{\ast,2}+k_{\textrm{F}}^{2})^{1/2}\approx k_{\textrm{F}}^{2}/2M_{0}^{\ast}+M_{0}^{\ast}\approx M_{0}^{\ast}$, indicating that the expression in the parentheses in Eq. (6.18) is small. As density increases, the $E_{\textrm{HO}}(\rho)$ by Eq. (6.18) increases eventually, leading to the increasing of the high order effects eventually. For example, keeping only the leading order term in $(k_{\textrm{F}}/M_{0}^{\ast})^{2}$, one obtains the high order effect as

[TABLE]

in msQCDSR. It shows that the high order effect $E_{\textrm{HO}}(\rho)$ is negative (see the left panel of Fig. 5).

Furthermore, when considering higher order terms in density in the chiral condensates, e.g.,

[TABLE]

and assuming other quark/gluon condensates remain unchanged, then the integration of EOS from Eq. (2.20) and Eq. (2.21) gives roughly $-({\theta+1})^{-1}b\zeta\rho^{\theta}$ to the symmetry energy from the $\rho^{\theta}$ term. And the corresponding term using the self-energy decomposition of the symmetry energy is given by $-(M_{0}^{\ast}/e_{\textrm{F}}^{\ast})b\zeta\rho^{\theta}/2$, consequently the difference becomes now

[TABLE]

the last line is valid at low densities where $M_{0}^{\ast}\approx e_{\textrm{F}}^{\ast}$ is a good approximation. It is obvious from this expression the $E_{\textrm{HO}}(\rho)$ will not be small even at low densities, and the high order effects become more and more important as density increases since when keeping term in $(k_{\textrm{F}}/M_{0}^{\ast})^{2}$, one has

[TABLE]

The first term is absent if $\theta=1$ and in this case the second term is relatively small owing to the small factor $(k_{\textrm{F}}/M_{0}^{\ast})^{2}$. On the other hand, the higher order terms in density in the chiral condensates are naturally essential as density increases, e.g., the linear approximation (3.126) breaks down eventually. Combining these analyses, it is intuitive to conclude that one of the main reasons for the breakdown of the parabolic approximation for the EOS of ANM in QCDSR maybe the higher order density terms in the chiral condensates, and this finding will be justified in the numerical studies in the following sections. The other effects, such as the continuum excitations, may also lead to the breakdown of the parabolic approximation.

VII Full Calculations on EOS of ANM

After three sections on the nucleon mass in vacuum, the self-energy structure of the symmetry energy and the high order effects of the EOS of ANM (mainly qualitatively) given above, we now systematically investigate the EOS of ANM through the full QCDSR calculations. It is necessary to mention again that the default parameters used are the same as these used in Fig. 5 except $\omega_{0}$, i.e., $m_{\textrm{q}}=3.5\,\textrm{MeV},m_{\textrm{s}}=95\,\textrm{MeV},f=0,\mathscr{M}^{2}=1.05\,\textrm{GeV}^{2}$, $\sigma_{\textrm{N}}=45\,\textrm{MeV}$, and $\omega_{0}=1.5\,\textrm{GeV}$. In this section, the dependence of the physical quantities on $\omega_{0},\mathscr{M}^{2},f$ and $\sigma_{\textrm{N}}$ will be studied carefully. Moreover, the Ioffe parameter $t\approx-1.22$ is determined in this section self-consistently via the scheme for the nucleon mass in vacuum given in Section IV. This section is organized as follows: The nucleon self-energies in SNM will be explored first in Subsection VII.1; then the symmetry energy with its self-energy decomposition will be studied in Subsection VII.2; Subsection VII.3 studies the correlation between the symmetry energy and the quark/gluon condensates, and finally the first QCDSR parameter set in this work, i.e., QCDSR-1 (or the naive QCDSR in the sense only the linear density terms in the chiral condensates (i.e., Eq. 3.137) are included), will finally be given in Subsection VII.4. Shortcoming of the QCDSR-1 together with the possible improvements will also be given.

VII.1 Nucleon Self-energies in SNM

In Fig. 10, the nucleon self-energies in SNM as functions of density with different values of $\omega_{0},\mathscr{M}^{2}$ or $f$, respectively, are shown. It is obvious that these three parameters all have obvious effects on the vector self-energy $\Sigma_{\textrm{V}}^{0}$, while $\omega_{0}$ and $\mathscr{M}^{2}$ have little impacts on the scalar self-energy $\Sigma_{\textrm{S}}^{0}$, and the effects of the parameter $f$ on both the $\Sigma_{\textrm{S}}^{0}$ and $\Sigma_{\textrm{V}}^{0}$ are shown to be sizable (similar phenomena are also displayed in Fig. 8).

Fig. 11 shows the nucleon Dirac effective mass in SNM $M_{\textrm{D}}^{\ast}\equiv M_{0}^{\ast}=M+\Sigma_{\textrm{S}}^{0}(\rho)$, as a function of density with different $f$, where the empirical constraint on the nucleon Dirac effective mass in SNM about $M_{\textrm{D}}^{\ast}(\rho_{0})/M\approx 0.6\pm 0.1$, e.g., by analyzing the energy level splitting data in several typical finite nuclei (e.g., see ref. Men06 ), is also shown. QCDSR gives a reasonable result for $M_{\textrm{D}}^{\ast}$, for instance, the $M_{\textrm{D}}^{\ast}(\rho_{0})/M$ is found to be about 0.45 for $f=0$, 0.53 for $f=0.25$, and 0.64 for $f=0.5$. In Fig. 12, the ratio of $\Sigma_{\textrm{V}}^{0}$ and $\Sigma_{\textrm{S}}^{0}$ as a function of density with different $f$ or $\omega_{0}$ is shown. Specifically, ${\Sigma_{\textrm{V}}^{0}(\rho)}/{\Sigma_{\textrm{S}}^{0}(\rho)}\approx-0.7\sim-0.3$, or ${\Sigma_{\textrm{V}}^{0}(\rho)}/{\Sigma_{\textrm{S}}^{0}(\rho)}\approx-0.5\pm 0.2$, within a wide range of densities, is obtained. The nearly constant ratio ${\Sigma_{\textrm{V}}^{0}(\rho)}/{\Sigma_{\textrm{S}}^{0}(\rho)}$ indicates the density dependence of the self-energies is almost linear. In Fig. 13 and Fig. 14, the $\Sigma_{\textrm{V}}^{0}$ and $\Sigma_{\textrm{S}}^{0}$ as functions of the Borel mass squared or the continuum excitation factor $\omega_{0}$ with different $f$ at the cross density $\rho_{\textrm{c}}=0.11\,\textrm{fm}^{-3}$ are shown. The $\mathscr{M}^{2}$–dependence of $\Sigma_{\textrm{V}}^{0}$ and $\Sigma_{\textrm{S}}^{0}$ at other densities is also studied, which shows very similar behaivor. Interestingly, the QCDSR window for the Borel mass squared $\mathscr{M}^{2}$, i.e., the region where the dependence of the quantities on $\mathscr{M}^{2}$ is weak even to be insensitive, from Fig. 13 is found to be roughly about $1\,\textrm{GeV}^{2}\lesssim\mathscr{M}^{2}\lesssim 1.5\,\textrm{GeV}^{2}$, which is consistent with the early studies on it, e.g., ref. Iof84 gave about $0.8\,\textrm{GeV}^{2}\lesssim\mathscr{M}^{2}\lesssim 1.4\,\textrm{GeV}^{2}$. Furthermore, the effects of $f$ on both the scalar and the vector self-energies are found to be obvious from Fig. 13, see also Fig. 10. Figs. 10, 11, 12, 13, and 14 together demonstrate that the four-quark condensates have sizable effects on the self-energies. Investigations on the four-quark condensates from more fundamental approaches (instead of only using the effective parameter $f$) thus will be extremely important for making further progress in QCDSR method. In Fig. 15, the density dependence of $\Sigma_{\textrm{V}}^{0}$ and $\Sigma_{\textrm{S}}^{0}$ with different $\sigma_{\textrm{N}}$ is shown. It is clearly shown that the sigma term $\sigma_{\textrm{N}}$ strongly affects the scalar self-energy. The reason is that the density dependence of $\Sigma_{\textrm{S}}^{0}$ is directly determined by the $\sigma_{\textrm{N}}$. For instance, in the msQCDSR, one has,

[TABLE]

i.e., $\Sigma_{\textrm{S}}^{0}(\rho)$ linearly decreases as the density $\rho$ increases. On the other hand, the effects of $\sigma_{\textrm{N}}$ on the vector self-energy $\Sigma_{\textrm{V}}^{0}$ are found to be essentially small, through the full QCDSR calculations. Actually, the $\Sigma_{\textrm{V}}^{0}$ is independent of the $\sigma_{\textrm{N}}$ in the msQCDSR.

In Fig. 16, the density dependence of $\Sigma_{\textrm{V}}^{0}$ and $\Sigma_{\textrm{S}}^{0}$ obtained order by order according to the same scheme used in Fig. 6 is shown. Since the density dependence of the quark/gluon condensates adopted in this section is linear (without the higher order density terms), the predictions on the dependence of the self-energies are also found to be roughly linear. The change on the scalar self-energy $\Sigma_{\textrm{S}}^{0}$ from msQCDSR to fQCDSR is relatively smaller than that on the vector self-energy $\Sigma_{\textrm{V}}^{0}$. For example, at $\rho_{0}=0.16\,\textrm{fm}^{-3}$, $\Sigma_{\textrm{S}}(\rho_{0})$ changes from $-464$ MeV in msQCDSR to $-506$ MeV in fQCDSR. However, for the $\Sigma_{\textrm{V}}(\rho_{0})$, it changes from 370 MeV in msQCDSR to 170 MeV in fQCDSR, inducing a $-54\%$ relative change.

Another feature in Fig. 16 is that the predictions on the scalar self-energy in the msQCDSR (black solid line) and in the sQCDSR (red dash line) are exactly the same. This can be proved as follows. In the msQCDSR, the expression for the effective mass is given by (in order to make a general proof here the Ioffe parameter can take any value instead of being set to $-1$, $\langle\overline{q}q\rangle_{\rho}$ is the chiral condensate at finite densities in SNM)

[TABLE]

or equivalently ${c_{1}}/{c_{4}}=-{\mathscr{M}^{2}{M}_{0}^{\ast}}/{16\pi^{2}\langle\overline{q}q\rangle_{\rho}}$, taking it at zero density leads to $c_{1}/c_{4}=-{\mathscr{M}^{2}{M}}/{16\pi^{2}\langle\overline{q}q\rangle_{\textrm{vac}}}$, thus one obtains the following formula for the nucleon scalar self-energy as a function of density,

[TABLE]

which is independent of the parameters $c_{1}$ and $c_{4}$. Similarly, in the sQCDSR, one has an extra term, i.e., $c_{1}\langle\overline{q}q\rangle_{\textrm{vac}}^{2}/6\equiv c_{1}\mathcal{I}$ in the QCDSR equations, and the corresponding effective mass is given by

[TABLE]

solving it gives

[TABLE]

or equivalently,

[TABLE]

thus

[TABLE]

i.e., the scalar self-energy $\Sigma_{\textrm{S}}^{0}(\rho)$ obtained in sQCDSR is still given by Eq. (VII.1), furnishing the proof. When other condensates are included in the QCDSR equations, however, there will be no close expression for the scalar self-energy.

VII.2 The Symmetry Energy

In this subsection, we generally study the self-energy decomposition of the symmetry energy through the full QCDSR equations without any fitting scheme. The five terms in Eq. (2.22) are the kinetic symmetry energy $E_{\textrm{sym}}^{\textrm{kin}}(\rho)$, contributions from momentum dependence of the scalar self-energy $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho)$ and that of the vector self-energy $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho)$, first order symmetry scalar self-energy $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ and the first order symmetry vector self-energy $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)$, respectively. In Fig. 17, each term as a function of density is shown. For example, at the cross density $\rho_{\textrm{c}}=0.11\,\textrm{fm}^{-3}$, one has $E_{\textrm{sym}}^{\textrm{kin}}(\rho_{\textrm{c}})\approx 14.9\,\textrm{MeV}$, $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho_{\textrm{c}})\approx-4.3\,\textrm{MeV}$, $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho_{\textrm{c}})\approx 2.5\,\textrm{MeV}$, $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho_{\textrm{c}})\approx 12.3\,\textrm{MeV}$, and similarly $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho_{\textrm{c}})\approx 18.6\,\textrm{MeV}$, leading to the total symmetry energy about $E_{\textrm{sym}}(\rho_{\textrm{c}})\approx 44.1\,\textrm{MeV}$. Obviously, the value of the total symmetry energy at $\rho_{\textrm{c}}$ is much larger than the empirical constraints on it, e.g., the one obtained by analyzing the binding energy difference between a heavy isotope pair gives $E_{\textrm{sym}}(\rho_{\textrm{c}})\approx 26.65\pm 0.2\,\textrm{MeV}$ Zha13 . As a reference, we also give these symmetry energy terms at the saturation density $\rho_{0}$, they are about 23.9 MeV, $-$5.0 MeV, 5.2 MeV, 14.7 MeV, and 25.4 MeV, leading to $E_{\textrm{sym}}(\rho_{0})\approx 64.2\,\textrm{MeV}$. Despite the discrepancy between the symmetry energy from the fQCDSR and the empirical constraints is large, QCDSR predicts the sign of the momentum-dependence of the nucleon self-energies, i.e.,

[TABLE]

and these relations will be further studied in Subsection VII.4 and the following sections when the fitting scheme (Subsection III.6) is adopted.

The slope parameter of the symmetry energy, i.e., $L(\rho)=3\rho{\text{d}E_{\textrm{sym}}(\rho)}/{\text{d}\rho}$ (it is easy to recognize meaning of $L$ from the context since we use the same letter for the slope parameter of the symmetry energy and the anomalous dimension of the interpolation field defined in Eq. (3.106)), can similarly be decomposed as different nucleon self-energy terms, see ref. Cai12 for more details. However, the most straightforward manner is to calculate each symmetry energy decomposition term, i.e., $L^{j}(\rho)=3\rho\text{d}E_{\textrm{sym}}^{j}(\rho)/\text{d}\rho$ with $j=$“kin”, “mom,0,S”, “mom,0,V”, “1st,S”, and “1st,V”. In Fig. 18, the density-dependence of these terms is shown, and for instance, at the cross density $\rho_{\textrm{c}}$, one obtains $L^{\textrm{kin}}(\rho_{\textrm{c}})\approx 51.4\,\textrm{MeV}$, $L^{\textrm{mom,0,S}}(\rho_{\textrm{c}})\approx-8.7\,\textrm{MeV}$, $L^{\textrm{mom,0,V}}(\rho_{\textrm{c}})\approx 13.7\,\textrm{MeV}$, $L^{\textrm{1st,S}}(\rho_{\textrm{c}})\approx 23.0\,\textrm{MeV}$, and $L^{\textrm{1st,V}}(\rho_{\textrm{c}})\approx 45.8\,\textrm{MeV}$, leading to $L(\rho_{\textrm{c}})\approx 125.1\,\textrm{MeV}$. The value of $L(\rho_{\textrm{c}})$ is similarly found to be much larger than the constraint on it, e.g., $L(\rho_{\textrm{c}})\approx 46\pm 4.5\,\textrm{MeV}$ from ref. Zha13 by analysing the correlation between the neutron skin thickness of the neutron-rich heavy nuclei and the $L$ parameter. As a reference, we also list these terms at the saturation density, they are 100.6 MeV, $-$4.9 MeV, 32.6 MeV, 14.9 MeV, and 65.9 MeV, respectively, and $L\equiv L(\rho_{0})\approx 209.2\,\textrm{MeV}$. Thus, the prediction on the $L$ parameter from the QCDSR without any fitting scheme is much larger than the empirical one both at the cross density and at the saturation density, for example, the $L$ parameter at $\rho_{0}$ was nowadays better constrained to be around about $60\pm 30$ MeV (see, e.g., refs. Che12a ; LiBA13 ; Oer17 ; Che17xx ).

In Fig. 19, the symmetry energy as a function of density with different $\omega_{0},\mathscr{M}^{2}$ or $f$ is shown. It is clearly shown from the figure that the effects of $\omega_{0}$ and $\mathscr{M}^{2}$ are comparatively larger than those of $f$. This could be easily understood, e.g., according to Eq. (VI), that the effects of the $f$–term on the symmetry energy is mainly characterized by the factor $\xi\equiv(\alpha/\beta)f\approx 0.1f$. Thus there is no surprise that the effects of the parameter $f$ on the symmetry energy is smaller than these on the EOS of SNM or that of PNM since the relevant term in $E_{0}(\rho)$ or $E_{\textrm{n}}(\rho)$ is directly proportional to $f$ or $(1-\alpha/\beta)f$, see the relevant analyses in Section VI, i.e., Eq. (VI) and Eq. (6.13). It means that the effects of $f$ on the symmetry energy is roughly cancelled, see the left panel of Fig. 9 for a similar calculation.

In Fig. 20, the density dependence of the self-energy decomposition of the symmetry energy is shown with different $\sigma_{\textrm{N}}$. The effects of $\sigma_{\textrm{N}}$ on the kinetic symmetry energy $E_{\textrm{sym}}^{\textrm{kin}}(\rho)$ is found to be large and it could be understood as follows: as $\sigma_{\textrm{N}}$ increases, then according to, e.g., Eq. (5.10) from the msQCDSR, one finds that the effective mass $M_{0}^{\ast}(\rho)$ is reduced. Consequently, $e_{\textrm{F}}^{\ast}$ is also reduced, leading to the enhancement on the kinetic symmetry energy $k_{\textrm{F}}^{2}/6e_{\textrm{F}}^{\ast}$. For instance, the $E_{\textrm{sym}}^{\textrm{kin}}(\rho_{0})$ changes from about 19.2 MeV at $\sigma_{\textrm{N}}=30\,\textrm{MeV}$ to 30.9 MeV at $\sigma_{\textrm{N}}=60\,\textrm{MeV}$, roughly a 61% relative increase. However, on the other hand, the factor $y\equiv M_{0}^{\ast}/e_{\textrm{F}}^{\ast}\approx 1-x^{2}/2$ decreases as $x\equiv k_{\textrm{F}}/M_{0}^{\ast}$ increases. Then according to Eq. (5.21), whether $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)\sim y\sigma_{\textrm{N}}$ is enhanced or reduced depends on the competition between the enhancement factor $\sigma_{\textrm{N}}$ and the reduction factor $y$, and the final effect is found to reduce the $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ at a fixed density, as shown in the panel (d) of Fig. 20, i.e., the factor $y$ wins the competition over $\sigma_{\textrm{N}}$. Besides $E_{\textrm{sym}}^{\textrm{kin}}(\rho)$ and $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$, the effects of $\sigma_{\textrm{N}}$ on the other three terms are very non-trivial and could not be analyzed in a semi-analytical manner. The $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho)$ and $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho)$, for instance, are found to be enhanced from about $-$4.8 MeV to $-$3.6 MeV and 2.2 MeV to 3.0 MeV as $\sigma_{\textrm{N}}$ changes from 30 MeV to 60 MeV, respectively, and the $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)$ is almost unaffected. Moreover, as shown in Fig. 20, the enhancement due to the increasing of $\sigma_{\textrm{N}}$ on the kinetic symmetry energy is much larger than the reduction on $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$, leading to the enhancement on the total symmetry energy, as shown in the left panel in Fig. 21. For example, if one takes $\sigma_{\textrm{N}}=45\pm 15\,\textrm{MeV}$ as shown in the left panel in Fig. 21, then the uncertainty on the symmetry energy is found to be about 3.5 MeV (12.3 MeV) at $\rho_{\textrm{c}}$ ($\rho_{0}$). Furthermore, the $\sigma_{\textrm{N}}$ affects the density dependence of the symmetry energy in the sense that a smaller $\sigma_{\textrm{N}}$ induces a softer symmetry energy, i.e., a smaller slope parameter $L$. Specifically, the relative change generated by varying $\sigma_{\textrm{N}}$ from 30 MeV to 60 MeV on the $L$ parameter shown in the left panel of Fig. 21 is about 73%. The connection between the nucleon-sigma term $\sigma_{\textrm{N}}$ and the parameter $L$ itself provides a possible mechanism to investigate the density dependence of the symmetry energy.

Interestingly, as shown in the right panel in Fig. 21, a larger mass of the s quark will induce a smaller symmetry energy. In fact, one can find that the effects of s quark mass $m_{\textrm{s}}$ come into play in the symmetry energy through the parameters $\alpha$ and $\beta$ defined in Eq. (3.121). For instance, taking a smaller (larger) value of the s quark mass, e.g., $m_{\textrm{s}}=60\,\textrm{MeV}$ ($m_{\textrm{s}}=130\,\textrm{MeV}$), one then obtains $\alpha\approx 0.29$ ($0.14$) and $\beta\approx 1.71$ ($1.86$), leading to $\alpha/\beta\approx 0.17$ ($0.08$). Then according to the expression for $\Sigma_{\textrm{sym}}^{\textrm{S}}(\rho)$ from the msQCDSR (see Eq. (5.16)), one can easily find that a smaller (larger) $m_{\textrm{s}}$ induces a larger (smaller) $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$, shown in the left panel in Fig. 22 (e.g., the $\Sigma_{\textrm{sym}}^{\textrm{1st,S}}(\rho_{\textrm{c}})$ changes from about 8.2 MeV to 21.5 MeV as $m_{\textrm{s}}$ changes from 130 MeV to 60 MeV, inducing an effect about 13.3 MeV).

On other hand, the s quark mass only affects the $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)$ in a more minor manner through the total QCDSR equations, shown in the right panel in Fig. 22 (the corresponding effect is found to be about 2.4 MeV at $\rho_{\textrm{c}}$), i.e., the change on the $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)$ due to the $m_{\textrm{s}}$ is relatively small compared to its effects on $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$. From the decomposition (2.22), $m_{\textrm{s}}$ obviously has no effects on the other three terms, i.e., $E_{\textrm{sym}}^{\textrm{kin}}(\rho)$, $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho)$ and $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho)$ (since these terms only involve the symmetric quantities). These analyses finally lead to the conclusion that considering the $m_{\textrm{s}}$ uncertainties will induce sizable effects on the total symmetry energy. To our best knowledge, since there exists no similar analyses on the relation between the symmetry energy and the s quark mass, our studies on it may provide new insights into the physical origin of as well as the uncertainties on the nuclear symmetry energy with respect to the s quark mass.

In Fig. 23, we show the density dependence of the symmetry energy obtained order by order. One can find that the four-quark condensates of the type (3.105) essentially have large effects on the symmetry energy (from the msQCDSR with black solid line to the sQCDSR with red dash line denoted by “$\mathscr{C}_{\textrm{B}}=1$”). For instance, the symmetry energy at $\rho_{\textrm{c}}$ ($\rho_{0}$) in the msQCDSR is found be about 61.2 MeV (88.9 MeV). On the other hand, the $E_{\textrm{sym}}$ in the sQCDSR at $\rho_{\textrm{c}}$ ($\rho_{0}$) is found to be about 49.8 MeV (72.3 MeV), generating a reduction about 11.4 MeV (16.6 MeV), respectively. Once Eq. (3.105) is included in the QCDSR equations, the following effects are found to be much smaller compared to those from Eq. (3.105). These estimates demonstrate again the importance of the four-quark condensates (3.105).

Finally, the self-energy decomposition of the symmetry energy with different $\omega_{0}$/$\mathscr{M}^{2}$/$f$ is shown in Fig. 24. The order by order results on these terms are also shown in the fourth row. Several features are necessary to be pointed from these figures: 1). At a fixed density, the effects of the continuum excitation factor $\omega_{0}$ on each decomposition term are positive correlated, i.e., a smaller (larger) $\omega_{0}$ corresponds to a smaller (larger) $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho)$ or $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho)$ or $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ or $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)$, and $\omega_{0}$ almost does not affect $E_{\textrm{sym}}^{\textrm{kin}}(\rho)$. 2). The effects of the Borel mass squared $\mathscr{M}^{2}$ on the $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho)$ or $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho)$ or $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ are similar like $\omega_{0}$. However, as $\mathscr{M}^{2}$ increases, $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ is shown to be reduced, and the overall effects of $\mathscr{M}^{2}$ are to enhance the symmetry energy (see also Fig. 19). 3). As the four-quark condensates parameter $f$ increase, the effective mass $M_{0}^{\ast}$ (and consequently $e_{\textrm{F}}^{\ast}$) increases according to Eq. (6.5), leading to the reduction on the kinetic symmetry energy at a fixed density, as shown in the left panel of row 3 in Fig. 24. 4). For the order by order calculations, the continuum excitation effects are found to be large on $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho)$, $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho)$, and $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$, see the black dot lines shown in the second, third and fourth panels (from left) in the fourth row of Fig. 24. Features shown in these figures clearly establish the close relationships between the self-energy decomposition terms of the symmetry energy and the effective parameters appearing in the QCDSR, providing important guidelines to understand the physical origin of the symmetry energy as well as the corresponding uncertainties.

VII.3 Correlation between $E_{\textmd{sym}}(\rho)$

and Condensates

In this subsection, we study the correlation of the nucleon self-energies and the symmetry energy with the quark/gluon condensates. These explorations are useful for generally establishing the connection between the EOS of ANM in nuclear physics, and the condensates encapsulating the degrees of freedom of quarks and gluons in hadronic physics (and related issues in QCD). The condensates and their uncertainties are listed as follows (see Subsection III.5),

[TABLE]

where the last term in each line represents the central value of the quantity, and

[TABLE]

Other parameters introduced in QCDSR are $\alpha_{\textrm{s}}=0.5,f=0$, as well as $\mu=0.5\,\textrm{GeV},\mathscr{M}^{2}=1.05\,\textrm{GeV}^{2},\omega_{0}=1.5\,\textrm{GeV},\sigma_{\textrm{N}}=45\,\textrm{MeV},m_{\textrm{q}}=3.5\,\textrm{MeV}$ and $m_{\textrm{s}}=95\,\textrm{MeV}$.

In general, one can in the framework of QCDSR study different correlations between the condensates and other physical quantities, through the QCDSR equations. For example, in Fig. 25, the correlation between the nucleon self-energies in SNM $\Sigma_{\textrm{S/V}}^{0}$(upper panel)/nuclear symmetry energy (low panel) at the cross density $\rho_{\textrm{c}}$ and the condensate properties is shown. Similarly, a very similar correlation pattern will be found if one studies the correlation at other densities, e.g., at $\rho_{0}$. For the scalar self-energy $\Sigma_{\textrm{S}}^{0}$, since, e.g., $\Sigma_{\textrm{S}}^{0}\propto-\langle\overline{q}q\rangle_{\textrm{vac}}$ in the msQCDSR (see Eq. (7.1)), a larger chiral condensate in vacuum naturally leads to a smaller $\Sigma_{\textrm{S}}^{0}$. Besides the strong dependence on $\langle\overline{q}q\rangle_{\textrm{vac}}$, both $\Sigma_{\textrm{S}}^{0}$ and $\Sigma_{\textrm{V}}^{0}$ are independent of the other condensate properties except that $\Sigma_{\textrm{S}}^{0}$ weakly depends on the five-dimensional parameter $\langle g_{\textrm{s}}\overline{q}\sigma\mathcal{G}q\rangle_{\textrm{sym}}^{\textrm{p}}$.

On the other hand, the dependence of the symmetry energy on the quark/gluon condensates show very fruitful patterns (see the lower panel of Fig. 25). For instance, the positive correlation between $E_{\textrm{sym}}(\rho)$ and $\langle\overline{q}q\rangle_{\textrm{vac}}$. This could be understood as follows: as the $\langle\overline{q}q\rangle_{\textrm{vac}}$ increases (e.g., moving from left to right on the horizon axes), the nucleon Dirac effective mass $M_{0}^{\ast}=M+\Sigma_{\textrm{S}}^{0}$ is reduced correspondingly (since $\Sigma_{\textrm{S}}^{0}$ increases), leading to an enhancement on the factor $e_{\textrm{F}}^{\ast,-1}=[k_{\textrm{F}}^{2}+(M+\Sigma_{\textrm{S}}^{0})^{2}]^{-1/2}$ at a fix density, thus consequently the kinetic symmetry energy $E_{\textrm{sym}}^{\textrm{kin}}(\rho)=k_{\textrm{F}}^{2}/6e_{\textrm{F}}^{\ast}$ is enhanced. While on the other hand, other symmetry energy decomposition terms (defined in Eq. (2.22)) depend weakly on the chiral condensate in vacuum. More specifically, the symmetry energy at $\rho_{\textrm{c}}$ is changed from 35.1 MeV at $\langle\overline{q}q\rangle_{\textrm{vac}}=-(255\,\textrm{MeV})^{3}$ to 61.0 MeV at $\langle\overline{q}q\rangle_{\textrm{vac}}=-(220\,\textrm{MeV})^{3}$, generating an enhancement about 25.9 MeV. More interestingly, the symmetry energy is found to depend on several other quantities characterizing the high mass dimensional condensates. And among which the symmetry energy displays the strongest correlation with the $\vartheta_{3}$ and the $\langle g_{\textrm{s}}q^{{\dagger}}\sigma\mathcal{G}q\rangle_{\textrm{sym}}^{\textrm{p}}$ parameters. It is not intuitive to understand this strong correlation, however, these two parameters together determine the density dependence of the mixing condensates of quarks and gluons, i.e., $\langle g_{\textrm{s}}q^{{\dagger}}\sigma\mathcal{G}q\rangle_{\rho,\delta}\approx(1\mp\vartheta_{3}\delta)\langle g_{\textrm{s}}q^{{\dagger}}\sigma\mathcal{G}q\rangle_{\textrm{sym}}^{\textrm{p}}$, see Eq. (3.136), and to our best knowledge, it is the first time to relate the nuclear symmetry energy to the quark and gluon mixing condensates, with the latter the very fundamental quantities in hadronic physics. Thus the strong connection between $E_{\textrm{sym}}(\rho)$ and $\langle g_{\textrm{s}}q^{{\dagger}}\sigma\mathcal{G}q\rangle_{\rho,\delta}$ provides a useful bridge to explore the properties of EOS of ANM using the knowledge from other physics branches (e.g., hadronic physics here).

VII.4 QCDSR-1 (naive QCDSR)

Based on the qualitative analyses given in the above sections and the full QCDSR calculations given in this section, we now give the first set of the QCDSR parameter to study the EOS of ANM. In this subsection and the following two sections, the main attentions will be given to the density dependence of the EOS of SNM, the symmetry energy and the EOS of PNM, while the other physical quantities such as dependence of the self-energies on the Borel mass squared, etc., will not be given (which actually could be obtained very similarly as the analyses in the above subsections). According to the fitting scheme given in Subsection III.6, the four-quark condensates parameter $f$ is found to be about $f\approx 0.50$. Interestingly, with only the $f$ parameter, the PNM EOS $E_{\textrm{n}}(\rho)$ at the density $\rho_{\textrm{vl}}\approx 0.02\,\textrm{fm}^{-3}$ and the symmetry energy $E_{\textrm{sym}}(\rho)$ at $\rho_{\textrm{c}}\approx 0.11\,\textrm{fm}^{-3}$ could be fitted reasonably within their empirical ranges (Subsection III.6), i.e., $E_{\textrm{n}}(\rho_{\textrm{vl}})\approx 4.2\,\textrm{MeV}$ Tew13 ; Kru13 and $E_{\textrm{sym}}(\rho_{\textrm{c}})\approx 26.65\pm 0.2\,\textrm{MeV}$ Zha13 . In the following, we call this QCDSR parameter set the QCDSR-1, or the naive QCDSR, in the sense that only the linear approximation of the chiral condensates is used (i.e., without the $\Phi$-term in Eq. (3.137)). Moreover, the twist-four four-quark condensates will not be included in this section and in the next section, and we explore this type of condensates on the EOS of ANM in some detail in Section IX.

In Fig. 26, the symmetry energy as well as its self-energy decomposition (see Eq. (2.22)) are shown. Specifically, the symmetry energy at cross density $\rho_{\textrm{c}}=0.11\,\textrm{fm}^{-3}$ is found to be about $E_{\textrm{sym}}(\rho_{\textrm{c}})\approx 26.9\,\textrm{MeV}$ (slightly larger than the one adopted in the fitting scheme), while the corresponding decomposition terms are found to be about $E_{\textrm{sym}}^{\textrm{kin}}(\rho_{\textrm{c}})\approx 11.6\,\textrm{MeV}$, $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho_{\textrm{c}})\approx-6.6\,\textrm{MeV}$, $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho_{\textrm{c}})\approx 0.6\,\textrm{MeV}$, $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho_{\textrm{c}})\approx 12.9\,\textrm{MeV}$, and $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho_{\textrm{c}})\approx 8.4\,\textrm{MeV}$, respectively. As a reference, the symmetry energy as well as the corresponding decomposition terms at $\rho_{0}$ are found to be about 45.1 MeV and 17.4 MeV, $-$9.1 MeV, 1.2 MeV, 21.4 MeV, and 14.2 MeV, respectively. It is obvious that the $E_{\textrm{sym}}(\rho_{0})$ from the naive QCDSR is still some larger than its empirical constraints (e.g., about $32\pm 3\,\textrm{MeV}$). As is shown in Fig. 26, the momentum dependence of the nucleon scalar (vector) self-energy is negative (positive) within the density region, which is consistent with the results obtained in the last section without adopting the fitting scheme. However, as mentioned in the above sections, the magnitude of the momentum dependence of the nucleon self-energies is essentially weak (both the blue and green lines), leading to a relative smaller contribution to the symmetry energy. For comparison, we also include in Fig. 26 the corresponding predictions from the phenomenological nonlinear RMF model with the celebrated parameter set FSUGold  Tod05 . Since there is no scalar and isovector channel (characterized by the so-called $\delta$ meson) in the FSUGold parameter set, the $E_{\textrm{sym}}^{\textrm{1st,S}}$ is identically zero, and since the conventional RMF approaches lack momentum-dependent interactions, the corresponding momentum-dependent terms $E_{\textrm{sym}}^{\textrm{mom,0,S/V}}$ are zero too. More specifically, the total symmetry energy in the FSUGold parameter set at the density $\rho_{\textrm{c}}=0.11\,\textrm{fm}^{-3}$/the density $\rho_{0}=0.16\,\textrm{fm}^{-3}$ is about 27.1 MeV /34.2 MeV. In addition, it is seen from Fig. 26 that the $E_{\textrm{sym}}^{\textrm{1st,V}}$ from FSUGold is significantly different from the prediction of QCDSR-1. Furthermore, the symmetry energy (2.22) could also be rewritten as Cai12 ; LCCX18

[TABLE]

where $M_{\textrm{L}}^{\ast}$ is the nucleon Landau effective mass in SNM, defined as LCCX18

[TABLE]

with $e_{\textrm{F}}(\rho,|\textbf{k}|)$ the total single nucleon energy in SNM. Consequently, from the sum of the kinetic symmetry energy and the two terms related to the momentum dependence of the self-energies, one can easily obtain the Landau mass. In QCDSR-1, the nucleon Landau mass in SNM is given by $M_{\textrm{L}}^{\ast}/M\approx 1.29$ at the saturation density (the Dirac mass will be given shortly). It is also interesting to notice that the overall density dependence of the symmetry energy is very different for the QCDSR-1 and FSUGold, reflecting that the slope parameter $L$ are different in these predictions. Specifically, the slope parameter of the symmetry energy could be obtained directly through its definition, and at $\rho_{\textrm{c}}$ the $L$ parameter $L(\rho_{\textrm{c}})$ in QCDSR-1 is found to be $105.9$ MeV, which is much larger than the empirical constraints, see, e.g., refs. Zha13 ; Dan14 . On the other hand, the $L(\rho_{\textrm{c}})$ in FSUGold is about $50.0$ MeV. While at $\rho_{0}=0.16\,\textrm{fm}^{-3}$, the $L$ parameter in QCDSR-1 (FSUGold) is found to be about $196.8$ MeV/ ($63.9$ MeV), once again showing the one obtained from QCDSR-1 is much larger than the phenomenological prediction. The large symmetry energy at the saturation density as well as the $L$ parameter both at the cross density and the saturation density indicate that the linear approximation of the condensates, especially the chiral condensates (3.137), already breaks down at densities significantly less than the saturation density $\rho_{0}$. It is necessary to consider the effective higher order terms in density in the chiral condensate, and this is the main task of the next section.

In Fig. 27, the total symmetry energy obtained in QCDSR-1 as a function of density is shown. Although the symmetry energy roughly passes through the constraint at $\rho_{\textrm{c}}\approx 0.11\,\textrm{fm}^{-3}$ Zha13 ; Dan14 ; Zha15 ; Bro14 , its density behavior both at low densities $\lesssim 0.05\,\textrm{fm}^{-3}$ and at densities $\gtrsim 0.12\,\textrm{fm}^{-3}$ is shown to be inconsistent with the predictions from other approaches in the sense that the symmetry energy obtained in QCDSR-1 is too stiff. Several other typical constraints on the symmetry energy are also shown in Fig. 27 include the results from the analysis of isobaric analog states (IAS) Dan14 while the shaded area enclosed by the solid cyan line and labeled by “IAS+NSkin” results when the IAS analysis is supplemented with additional constraints from neutron skin data, the constraint on the symmetry energy from heavy ion collisions (HIC) which is labeled by the blue band Tsa12 . A recent study on the symmetry energy at low densities (around $\rho_{0}/3$) using the electric dipole polarizability in 208Pb is also shown for comparison Zha15 by the magenta band. Two studies on the symmetry energy at the cross density Zha13 are shown, they are the constrain on the symmetry energy at $\rho_{\textrm{c}}\approx 0.11\,\textrm{fm}^{-3}$ to be $E_{\textrm{sym}}(\rho_{\textrm{c}})\approx 26.65\pm 0.20\,\textrm{MeV}$ using isotope binding energy difference Zha13 which is labeled by a magenta star and that from fit of ground state properties of double magic nucleus using Skyrme CSkp functionals which found a value of $E_{\text{sym}}(\rho_{\text{c}})\approx 25.4\pm 0.8\,\text{MeV}$ at $\rho_{\text{c}}\approx 0.10\,\text{fm}^{-3}$ Bro14 (blue circle). Moreover, since we are mainly interested in the symmetry energy from QCDSR method at densities $\lesssim\rho_{0}$, some other constraints on $E_{\textrm{sym}}(\rho)$ at densities $\gtrsim\rho_{0}$ will not be compared here, e.g., ref. Rus16 constraints the symmetry energy between $\rho_{0}$ to about $2\rho_{0}$.

In Fig. 28, we show the EOS of PNM as a function of density. Also included in Fig. 28 are the results from ChPT Tew13 ; Kru13 (green band), QMC simulations combined with chiral force to next-to-next-to-leading order (N2LO) with Tew16 (blue band) and without Gez13 (magenta band) leading-order chiral three-nucleon interactions forces, next-to-leading order (NLO) lattice calculation Epe09a (magenta circle), and QMC simulations for PNM at very low densities Gez10 (green diamond). The result from analyzing experimental data on the electric dipole polarizability $\alpha_{\text{D}}$ in 208Pb Zha15 is also shown for comparison. The inset in Fig. 28 shows the EOS of PNM at very low densities. Interestingly, by artificially neglecting the contributions from dimension-four and higher order terms (and only keeping the four-quark condensates of the type (3.105)), one obtains an effective approximation for EOS of PNM, see Eq. (6.13) and Eq. (6.16), and the latter (Eq. (6.16)) clearly demonstrates how the chiral condensate goes into play in the EOS of PNM, i.e., the second term characterized by several constants ($\xi,\sigma_{\textrm{N}},m_{\textrm{q}}$ and $\langle\overline{q}q\rangle_{\textrm{vac}}$) is negative, leading to a reduction on the $E_{\textrm{n}}(\rho)$ compared to the FFG prediction. In Fig. 28, we also plot the results from Eq. (6.16) at densities $\lesssim 0.02\,\textrm{fm}^{-3}$ (violet solid square). One can see that the approximation Eq. (6.16) can already produce reasonably the $E_{\textrm{n}}(\rho)$ at low densities. The FFG prediction on the $E_{\textrm{n}}^{\textrm{FFG}}(\rho)$, however, already becomes larger at, e.g., $\rho_{\textrm{vl}}=0.02\,\textrm{fm}^{-3}$, and even a relativistic correction $-k_{\textrm{F,n}}^{4}/56M^{3}$ to the $E_{\textrm{n}}^{\textrm{FFG}}(\rho)$ will not make the situation much better, strongly indicating that the non-interacting Fermi model lacks the fundamental information between nucleons and nucleons to produce the correct EOS. Furthermore, it is seen from Fig. 28 that the prediction on the $E_{\textrm{n}}(\rho)$ from QCDSR-1 is consistent with several QMC simulations and lattice computation at densities $\lesssim 0.02\,\textrm{fm}^{-3}$, showing that QCDSR is a reliable approach in the study of EOS of PNM, especially at lower densities, where the naive QCDSR is good enough.

Although the only adjustable parameter $f$ in the naive QCDSR is fixed by the $E_{\textrm{n}}(\rho)$ at the very low density of $0.02\,\textrm{fm}^{-3}$ (thus the $E_{\textrm{n}}(\rho)$ at densities greater than $0.02\,\textrm{fm}^{-3}$ has no fitting requirements), the prediction on the EOS of PNM at densities $\lesssim 0.1\,\textrm{fm}^{-3}$ in QCDSR-1 is found to be well-behaved compared with the APR EOS, demonstrating that the QCDSR with the linear density approximation for the chiral condensates can be quantitatively applied to study the EOS of PNM within these density region. However, as density even increases, the systematic deviation between the $E_{\textrm{n}}(\rho)$ obtained by QCDSR-1 and that predicted by APR EOS becomes large and this can not be improved by simply adjusting the parameter $f$, indicating on the other hand that the leading-order linear density approximation for the chiral condensates does not work well enough and the higher order density terms in the chiral condensates are needed for the PNM calculations at these densities. For example, at $\rho_{0}\approx 0.16\,\textrm{fm}^{-3}$, the difference between the ARP EOS and the $E_{\textrm{n}}(\rho)$ from QCDSR-1 is found to be about 5.4 MeV. Once one considers the term $\Phi g\rho^{2}$ in Eq. (3.137) for PNM, and recalculate the $E_{\textrm{n}}(\rho)$ under the fitting scheme, we find that compared with the case of the naive QCDSR, the obtained prediction can be largely improved to fit the APR EOS. This feature suggests that the QCDSR with effective higher order density terms in quark condensates can be used to study the EOS of dense nucleonic matter at higher densities. The relevant investigations will be given in the next section.

Since in our fitting scheme (see Subsection III.6), the constraints on the $E_{\textrm{n}}(\rho)$ at a very low density $\rho_{\textrm{vl}}$ and on the symmetry energy at $\rho_{\textrm{c}}$ are used without using any empirical constraints on the EOS of SNM (even the symmetry energy constraint is not used in QCDSR-1), it is interesting to explore the consequent prediction on the $E_{0}(\rho)$ from the QCDSR-1. In Fig. 29, we show the density dependence of the EOS of SNM obtained in QCDSR-1. The empirical constraints on the saturation density about $\rho_{0}\approx 0.16\pm 0.02\,\textrm{fm}^{-3}$ and the corresponding binding energy $E_{0}(\rho_{0})\approx-16\pm 2\,\textrm{MeV}$ are also shown in Fig. 29. For instance, the $E_{0}(\rho)$ at $0.16\,\textrm{fm}^{-3}$ ($0.11\,\textrm{fm}^{-3}$) in QCDSR-1 is found to be about $-$34.0 MeV ($-$19.7 MeV). And in the meanwhile, the saturation density of QCDSR-1 itself is found to be about $\rho_{0}^{\textrm{QCDSR-1}}\approx 0.6\,\textrm{fm}^{-3}$, with the corresponding binding energy about $-$100.0 MeV, showing that the symmetric matter in QCDSR-1 is very deep bounded. It actually is another cue that the effective $\Phi$-term in Eq. (3.137) is important, indicating the breakdown of the chiral condensates at linear order at densities even smaller than the saturation density. Moreover, it is really a very difficult problem on how to obtain the correct (even reasonable) saturation properties of the SNM in the microscopic theories (see, e.g., ref. Dri19 in the framework of ChPT). Based on the symmetry energy and the EOS of SNM and the EOS PNM, one can estimate the high order effects in the EOS $E_{\textrm{HO}}(\rho)$, and the detailed results on $E_{\textrm{HO}}(\rho)$ will not be given here (see Fig. 36). However as a reference, for example, the $E_{\textrm{HO}}(\rho)$ is found to be about 3.6 MeV (1.1 MeV) at $\rho_{\textrm{c}}$ ($\rho_{0}$), demonstrating again that the $E_{\textrm{HO}}(\rho)$ is generally non-negligible in QCDSR.

In Fig. 30, the nucleon self-energies in SNM and the nucleon Dirac effective mass as functions of density are shown. The predictions on the self-energies by the Dirac–Brueckner–Hartree–Fock (DBHF) calculations Bro90 are also shown for comparison. It is obvious from the figure that although the overall tendency of the density dependence of the self-energies is consistent compared to the DBHF predictions, the density behavior roughly characterized, e.g., by the index $\sigma$ in $\rho^{\sigma}$, is very different. It demonstrates onces again that the higher density terms in the chiral condensates are important, since $\Sigma_{\textrm{S}}^{0}(\rho)$ is directly related to the chiral condensates $\langle\overline{q}q\rangle_{\rho}$. Moreover, the nucleon Dirac effective mass at $\rho_{0}=0.16\,\textrm{fm}^{-3}$ is found to be $M_{\textrm{D}}^{\ast}/M\approx 0.65$, which is consistent with the empirical constraints LCK08 (since from Fig. 30 it is clearly shown that around $0.16\,\textrm{fm}^{-3}$ the prediction on $\Sigma_{\textrm{S}}^{0}$ is consistent with the DBHF prediction).

As a short summary of this section, we determine the effective parameter $f$ in the QCDSR to fix the $E_{\textrm{n}}(\rho)$ at $\rho_{\textrm{vl}}\approx 0.02\,\textrm{fm}^{-3}$ with the prediction from ChPT, and the symmetry energy at the cross density $\rho_{\textrm{c}}$ is found to be consistent with the empirical constraints. However, the symmetry energy at the saturation density and the slope parameter $L$ at both the saturation and the cross density are found to be too large, clearly indicating the breakdown of the linear approximation for the chiral condensates. Besides, the EOS of PNM at densities larger around $0.1\,\textrm{fm}^{-3}$ also shows systematic deviation from the APR EOS. Improving the EOS of PNM and the density behavior of the symmetry energy is one of the main motivations to include higher order terms in density in the chiral condensates, which is the main issue in the next section.

VIII Higher Order Density Terms in Chiral Condensates, Symmetry Energy and QCDSR-2

As studied in the last section, the symmetry energy and its slope parameter in QCDSR-1 are found to be too large at the saturation density, indicating the linear approximation for the chiral condensates (i.e., Eq. (3.126)) breaks down already at densities less than the saturation point since the density dependence of the chiral condensates is directly related to the symmetry energy, e.g., see the expression for the symmetry energy obtained in the msQCDSR, i.e., Eq. (5.19) and Eq. (5.21). In this section, we study the possible higher order terms in density in the chiral condensates introduced by the $\Phi$-term in Eq. (3.137). Once the $\Phi$-term is included in the chiral condensates, two other effective parameters, i.e., $\Phi$ and $g$ are introduced, and $\Phi$, $g$ and $f$ will be determined by the $E_{\textrm{n}}(\rho)$ at a very low density $\rho_{\textrm{vl}}=0.02\,\textrm{fm}^{-3}$, the symmetry energy at $\rho_{\textrm{c}}=0.11\,\textrm{fm}^{-3}$, and adjusting the PNM EOS to fit the APR EOS as much as possible, see the fitting scheme given in Subsection III.6. Consequently, these parameters will be re-adjusted. Besides, other parameters in the QCDSR are still taken as $\mathscr{M}^{2}=1.05\,\textrm{GeV}^{2},\omega_{0}=1.5\,\textrm{GeV},\sigma_{\textrm{N}}=45\,\textrm{MeV},m_{\textrm{q}}=3.5\,\textrm{MeV}$ and $m_{\textrm{s}}=95\,\textrm{MeV}$. One then obtains $f\approx 0.43$, $\Phi^{\prime}\equiv\Phi\times\langle\overline{q}q\rangle_{\textrm{vac}}\approx 3.45$ Cai17QCDSR , and $g=-0.64$. Moreover, the Ioffe parameter $t\approx-1.22$ is independent of the $\Phi$-term introduced, since when using the nucleon mass in vacuum to fix the parameter $t$ it only depends on the vacuum properties of the condensates.

We abbreviate the corresponding QCDSR parameter set the QCDSR-2. Based on the obtained $\Phi$ and $g$, we can estimate the density below which the $\Phi$-term has minor contribution to the quark condensates. This density can be estimated as $|\Phi(1-g)\rho^{2}|\ll|(\sigma_{\textrm{N}}/2m_{\textrm{q}})(1-\xi)\rho|$, i.e., the last term in Eq. (3.137) is significantly less than the second term in Eq. (3.137), and we thus obtain $\rho\ll\rho_{\textrm{es}}\approx 2.13\,\textrm{fm}^{-3}$. Therefore, the effects of $\Phi$ and $g$ on the $E_{\textrm{n}}(\rho)$ are trivial at substantial densities, e.g., when one artificially takes $\Phi^{\prime}=0$ and keeping $f$ fixed, the $E_{\textrm{n}}(\rho_{\textrm{vl}})$ $(E_{\textrm{n}}(0.1\,\textrm{fm}^{-3})$ changes from 4.20 MeV to 4.22 MeV (from 11.15 MeV to 10.01 MeV). Thus it is reasonable to expect that effects of $\Phi$ and $g$ on the $E_{\textrm{n}}(\rho)$ at low densities $\lesssim 0.1\,\textrm{fm}^{-3}$ are small. However, as the densities even increase, there is no guarantee that the $\Phi$-term still has small effects on the EOS of PNM since the $E_{\textrm{n}}(\rho)$ is obtained by integrating over the density, see Eq. (2.21) Cai17QCDSR . In Fig. 31, the $E_{\textrm{n}}(\rho)$ obtained in QCDSR-2 is shown, and the prediction on the EOS of PNM from QCDSR-1 is also shown for comparison. Other constraints on $E_{\textrm{n}}(\rho)$ shown in the figure are the same as those in Fig. 28. Compared with the prediction on the $E_{\textrm{n}}(\rho)$ by the naive QCDSR, once one considers the term $\Phi g\rho^{2}$ in Eq. (3.137) for PNM, and recalculates the corresponding EOS, one finds that the obtained prediction is largely improved to fit the APR EOS. For instance, the EOS of PNM at $0.12\,\textrm{fm}^{-3}$ is now found to be 12.9 MeV, which is very close to the APR prediction 13.3 MeV. Moreover, the overall fitting between the EOS of PNM from QCDSR-2 and the APR EOS is much better compared with the prediction by QCDSR-1, e.g., at densities $\lesssim 0.16\,\textrm{fm}^{-3}$. These features suggest that the QCDSR with effective higher order density terms in quark condensates can be used to study the EOS of dense nucleonic matter at higher densities. It is also necessary to point out that using a different higher order density term in Eq. (3.137) and re-fix the parameters $f$, $\Phi$ and $g$ by the same fitting scheme, the density behavior of the $E_{\textrm{n}}(\rho)$ is almost unchanged. For example, when adopting a $\rho^{5/3}$ term, i.e., $\Phi(1\mp g\delta)\rho^{5/3}$, then $f\approx 0.46$, $\Phi^{\prime}\equiv\Phi\times\langle\overline{q}q\rangle_{\textrm{vac}}^{2/3}\approx 1.61$ and $g\approx-0.34$ could be obtained, and the corresponding $E_{\textrm{n}}(\rho)$ is shown in Fig. 31 by the magenta dot line. It is obvious seen that using a different higher order density term in the chiral condensates will not change our conclusions on the EOS of PNM Cai17QCDSR . In the following, we will not consider other higher order density terms in Eq. (3.137) except the one has the form $\Phi(1\mp g\delta)\rho^{2}$.

In Fig. 32, the density dependence of the symmetry energy obtained in QCDSR-2 is shown, and the $E_{\textrm{sym}}(\rho)$ from QCDSR-1 is also shown here for comparison. Compared to the prediction on the symmetry energy from QCDSR-1, the $E_{\textrm{sym}}(\rho)$ from QCDSR-2 is now largely improved, e.g., from about $0.04\,\textrm{fm}^{-3}$ to about the saturation density $\rho_{0}$. For example, the $E_{\textrm{sym}}(\rho)$ now can roughly pass through the constraints obtained from the IAS studies (green band). More specifically, when evaluating at $\rho\approx 0.04\,\textrm{fm}^{-3}$, the symmetry energy changes from 8.0 MeV in the naive QCDSR to about 10.8 MeV in QCDSR-2, introducing a relative 35% change, and more interestingly now $E_{\textrm{sym}}(0.04\,\textrm{fm}^{-3})$ is very close to the lower limit predicted by the IAS studies about 11.1 MeV. Similarly, the symmetry energy at $\rho_{0}=0.16\,\textrm{fm}^{-3}$ changes from 45.1 MeV in QCDSR-1 to about 35.3 MeV in QCDSR-2, and the relative change is about $-$22%. All these features demonstrate that the higher order density terms in the chiral condensates (3.137) are essentially needed to describe a reasonable density behavior of the symmetry energy.

Similarly, the nucleon self-energy decomposition (Eq. (2.22)) of the symmetry energy obtained in QCDSR-2 is shown in Fig. 33. The corresponding results from the RMF model with FSUGold  Tod05 are also included for comparison. At the cross density $\rho_{\textrm{c}}=0.11\,\textrm{fm}^{-3}$, for instance, one now obtains $E_{\textrm{sym}}^{\textrm{kin}}(\rho_{\textrm{c}})\approx 11.7\,\textrm{MeV}$, $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho_{\textrm{c}})\approx-5.3\,\textrm{MeV}$, $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho_{\textrm{c}})\approx 0.8\,\textrm{MeV}$, $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho_{\textrm{c}})\approx 7.4\,\textrm{MeV}$, and $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho_{\textrm{c}})\approx 12.1\,\textrm{MeV}$, and consequently $E_{\textrm{sym}}(\rho_{\textrm{c}})\approx 26.7\,\textrm{MeV}$ (the fitting scheme). Moreover, the symmetry energy at the saturation density is found to be about 35.2 MeV with its self-energy decomposition terms about 16.6 MeV, $-$7.8 MeV, 1.5 MeV, 3.0 MeV, and 21.8 MeV, respectively. Besides the weak momentum dependence of the self-energies in SNM revealed several times in the above sections, one can find interestingly that the $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ starts to decrease at a critical density about $0.09\,\textrm{fm}^{-3}$ and to even be negative around the saturation density. This change could be traced back to the effective $\Phi$-term introduced in the chiral condensates (3.137). These features demonstrate again the importance of the higher order terms in density in the chiral condensates. It is necessary to point out that using a different higher order terms in density (such as a term proportional to $\rho^{4/3}$ or $\rho^{5/3}$ as discussed in Subsection III.6) and refits the parameters of the model, the changes in the density behavior of the symmetry energy will be almost the same as the one shown in Fig. 33, i.e., the $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ starts to decreases at a critical density and even to become negative. Moreover, the $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)$ becomes dominant at densities $\gtrsim 0.1\,\textrm{fm}^{-3}$ compared with the kinetic symmetry energy, e.g., $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho_{0})\approx 22.0\,\textrm{MeV}$ and $E_{\textrm{sym}}^{\textrm{kin}}(\rho_{0})\approx 16.7\,\textrm{MeV}$, leading to $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho_{0})/E_{\textrm{sym}}^{\textrm{kin}}(\rho_{0})\approx 1.32$. Furthermore, as the effective $\Phi$-term is included, the density dependence of the symmetry energy obtained from QCDSR-2 is found to be closer to the one from the phenomenological FSUGold parameter set. However, as in the case of QCDSR-1, the $E_{\textrm{sym}}^{\textrm{1st,V}}$ from FSUGold is again significantly different from the prediction of QCDSR-2. Compared with the QCDSR approach, therefore, the phenomenological nonlinear RMF model with FSUGold exhibits very different self-energy decomposition of the symmetry energy.

The $L$ parameter could be obtained through the density behavior of the symmetry energy. Specifically, one finds that at $\rho_{\textrm{c}}$ ($\rho_{0}$) the $L$ parameter is about 64.7 MeV (67.5 MeV). Compared with the predictions on the $L$ parameter by the QCDSR-1, i.e., 105.9 MeV (at $\rho_{\textrm{c}}$) and 196.8 MeV (at $\rho_{0}$), the $\Phi$-term introduces a relative amount about $-$39% ($-$66%) on the $L$ at $\rho_{\textrm{c}}$ ($\rho_{0}$). Interestingly, although the $L(\rho_{\textrm{c}})$ is slightly larger than nowadays best empirical constraints (e.g., $L(\rho_{\textrm{c}})$ about $46\pm 4.5\,$MeV from ref. Zha13 ), the $L$ parameter at the saturation density is found to be consistent with its empirical value about $60\pm 30\,\textrm{MeV}$ (e.g., see refs. Che12a ; LiBA13 ; Oer17 ). Moreover, a relevant quantity for the discussion on the density behavior of the symmetry energy is given by ref. Dan14

[TABLE]

The $\gamma$ parameter introduced in the parameterized form of the symmetry energy, e.g., $E_{\textrm{sym}}(\rho)\sim(\rho/\rho_{0})^{\gamma}$, is often used in heavy-ion collisions simulations Tsa12 . For instance, the $\gamma$ parameter of the potential part was constrained to be about $\gamma\approx 0.72\pm 0.19$ in ref. Rus16 from the comparison of the elliptic flow ratio of neutrons with respect to charged particles based on UrQMD predictions. In our calculation in QCDSR-2, we find that the $\gamma$ parameter is about 0.81 (0.64) at $\rho_{\textrm{c}}$ ($\rho_{0}$). Ref. Dan14 gives the constraint on $\gamma$ obtained by the IAS analyses, shown in Fig. 34. It is interesting to see that the $\gamma$ parameter obtained in QCDSR-2 at $\rho_{\textrm{c}}$ is shown to be slightly larger than the blue band, however, the $\gamma$ parameter at $\rho_{0}$ is well consistent with the one given by ref. Dan14 .

In Fig. 35, we show the EOS of SNM obtained in QCDSR-2. Due to the effective $\Phi$-term introduced in the chiral condensate, the saturation density of QCDSR-2 is improved to be about $0.2\,\textrm{fm}^{-3}$, with the corresponding binding energy about $-$26.7 MeV.

The improvements on the EOS of PNM, the EOS of SNM and the symmetry energy together verify again that the higher order terms in density in the chiral condensates are important. Similarly in Fig. 36, the $E_{\textrm{HO}}(\rho)\approx E_{\textrm{sym,4}}(\rho)+E_{\textrm{sym,6}}(\rho)+\cdot\cdot\cdot$ as a function of density is shown in both the QCDSR-1 and QCDSR-2. As discussed in Section VI, the higher order terms in density in the chiral condensates may eventually break the parabolic approximation of the EOS of ANM, see the discussion given around Eq. (6.21), we find in QCDSR-2 the $E_{\textrm{HO}}(\rho)$ is about 7.1 MeV (4.4 MeV) at $\rho_{0}$ ($\rho_{\textrm{c}}$). Although the high order term $E_{\textrm{HO}}(\rho)$ is generally believed to have very little effects on, e.g., the nuclear structure problems, it may induce sizable influence on the quantities in neutron stars Cai12xx ; PuJ17 , such as the core-crust transition density and the transition pressure. Detailed investigation on the $E_{\textrm{HO}}(\rho)$ in QCDSR will be extremely useful for further studies on the relevant issues, which however is beyond the main scope of the present work.

In Fig. 37, the density dependence of the nucleon self-energies in SNM is shown. Compared with the density behavior of the self-energies obtained in QCDSR-1 shown in the left panel of Fig. 30, the overall density behavior is largely improved in QCDSR-2. For instance, the $\Sigma_{\textrm{S}}^{0}$ ($\Sigma_{\textrm{V}}^{0}$) at $\rho_{0}\approx 0.16\,\textrm{fm}^{-3}$ is now shown to be about $-$312.3 MeV (258.5 MeV), which is close to the prediction by the DBFH theories about $-$339.7 MeV (267.7 MeV). Moreover, the nucleon Dirac effective mass in SNM is found to be about $M_{\textrm{D}}^{\ast}(\rho_{0})/M\approx 0.68$, via the $\Sigma_{\textrm{S}}^{0}(\rho)$ since $M_{\textrm{D}}^{\ast}=M+\Sigma_{\textrm{S}}^{0}$. Furthermore, the nucleon Landau effective mass could be obtained by the sum of $E_{\textrm{sym}}^{\textrm{kin}}(\rho)$, $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho)$, and $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho)$, see Eq. (7.23), consequently $M_{\textrm{L}}^{\ast}(\rho_{0})/M\approx 1.19$ through the decomposition of the symmetry energy (see Fig. 33).

Since in QCDSR-2, the parameters $\Phi$ and $g$ are fixed, one can inversely study their effects on the chiral condensates (3.137). In Fig. 38, the density dependence of the chiral condensates in QCDSR-2 as well as the corresponding predictions from ChPT Kru13a ; Lac10 ; Kai09 and the FRG method Dre14 is shown Cai17QCDSR . As demonstrated in Eq. (3.137) and shown in the figure, the chiral condensate at low densities is dominated by the linear density term. More specifically, one has $(\langle\overline{\textrm{u}}\textrm{u}\rangle_{\rho}-\langle\overline{\textrm{d}}\textrm{d}\rangle_{\rho})/\langle\overline{q}q\rangle_{\textrm{vac}}\approx-\rho\sigma_{\textrm{N}}\xi/m_{\textrm{q}}\langle\overline{q}q\rangle_{\textrm{vac}}>0$ at low densities, since $\langle\overline{q}q\rangle_{\textrm{vac}}$ is negative. As density increases, the $\Phi$-term in Eq. (3.137) starts to dominate and even to flip the relative relation of the magnitude between $\langle\overline{\textrm{u}}\textrm{u}\rangle_{\rho}$ and $\langle\overline{\textrm{d}}\textrm{d}\rangle_{\rho}$, leading to $\langle\overline{\textrm{u}}\textrm{u}\rangle_{\rho}/\langle\overline{q}q\rangle_{\textrm{vac}}<\langle\overline{\textrm{d}}\textrm{d}\rangle_{\rho}/\langle\overline{q}q\rangle_{\textrm{vac}}$ when the density is larger than about $0.15\,\textrm{fm}^{-3}$. For instance, the $\langle\overline{\textrm{d}}\textrm{d}\rangle_{\rho_{0}}/\langle\overline{q}q\rangle_{\textrm{vac}}\,(\langle\overline{\textrm{u}}\textrm{u}\rangle_{\rho_{0}}/\langle\overline{q}q\rangle_{\textrm{vac}})$ in PNM is found to change from 0.45 (0.56) in the linear density approximation to 0.60 (0.59) with the inclusion of the $\Phi$-term, leading to an enhancement of about 33% (5%). It is necessary to point out that this flip is a direct consequence of the inclusion of the higher order $\Phi$-term in Eq. (3.137).

More interestingly, one can find that the higher order $\Phi$-term in Eq. (3.137) stabilizes the chiral condensate both for u and d quarks at larger densities, while the leading-order linear density approximation Eq. (3.137) leads to chiral symmetry restoration at a density about $2\rho_{0}$ Cai17QCDSR . The hindrance of the chiral symmetry restoration due to the higher order density terms in quark condensates has important implications on the physical degrees of freedom in the core of compact stars such as the neutron stars/quark stars where the matter is very close to the PNM. These features are consistent with a recent analysis on the same issue based on the FRG method Dre14 .

Furthermore, if one uses a different $\sigma_{\textrm{N}}$, then by readjusting the values of the parameters $f$, $\Phi$ and $g$ based on the fitting scheme (see Subsection III.6), it can be demonstrated that the $\sigma_{\textrm{N}}$ has very little influence on the EOS of ANM. Different $\sigma_{\textrm{N}}$ leads to different $\Phi$ and $g$, but the density dependence of the chiral condensates will change only quantitatively, instead of qualitatively since the $\sigma_{\textrm{N}}$ term (linear order) is a perturbation to the vacuum chiral condensates, and similarly, the $\Phi$-term is a perturbation to the linear term. It should be pointed out that the exploration on the $\sigma_{\textrm{N}}$ itself is an important issue in nuclear physics, and the exact knowledge on the $\sigma_{\textrm{N}}$ will certainly help improving our understanding on the relevant aspects of the in-medium strong interaction.

Now, let us vary the values of the parameters $\Phi$ and $g$, while at the same time keep all the other parameters fixed, and meanwhile the EOS of PNM at $0.02\,\textrm{fm}^{-3}$ and the symmetry energy at $0.11\,\textrm{fm}^{-3}$ are fixed by the fitting scheme (Subsection III.6). Since $\Phi$ and $g$ essentially have no effect on $E_{\textrm{n}}(\rho_{\textrm{vl}})$ and the parameter $f$ is determined by $E_{\textrm{n}}(\rho_{\textrm{vl}})$, thus $f\approx 0.43$ remains unchanged. Moreover, fixing of the symmetry energy at $\rho_{\textrm{c}}$ indicates that $\Phi g\approx\textrm{const.}$, since the effects of $\Phi$ and $g$ on the EOS of SNM and that of PNM are roughly given by

[TABLE]

and thus $\Phi$ and $g$ contribute to the symmetry energy roughly as $\int\Phi g\rho^{2}\mathrm{d}\rho$. In Fig. 39, we show the $E_{\textrm{n}}(\rho)$ with the upper (lower) black dash line corresponding to $\Phi^{\prime}\approx 4.00,g\approx-0.552$ ($\Phi^{\prime}\approx 2.70,g\approx-0.817$), by maximally expanding the greed band predicted by the ChPT. In this sense we can study the extra constraints on the $E_{\textrm{n}}(\rho)$ from QCDSR. It is at this time necessary to point out that although the error on $E_{\textrm{n}}(\rho_{\textrm{vl}})$ is relatively small by the state-of-the-art microscopic many-body calculations and simulations, it still has some uncertainties (e.g., the uncertainties generated by the nucleon-sigma term $\sigma_{\textrm{N}}$). However in our scheme, the $E_{\textrm{n}}(\rho_{\textrm{vl}})$ is fixed at a certain value, thus the constraints on the $E_{\textrm{n}}(\rho)$ given in the following paragraphs should be thought only as a roughly estimate. On the other hand, it is useful to study the saturation properties of the EOS of SNM, since fixing the $E_{\textrm{n}}(\rho_{\textrm{vl}})$ actually gives a relevant estimate on the uncertainties on the saturation density of $E_{0}(\rho)$ due to the higher order density terms in Eq. (3.137). It is clearly seen from Fig. 39 that in such a way, a much stronger constraint on the EOS of PNM in the density region from about $0.04\,\textrm{fm}^{-3}$ to $0.12\,\textrm{fm}^{-3}$ is obtained (the inset of Fig. 39). For example, one obtains $6.8\,\textrm{MeV}\,(10.0\,\textrm{MeV})\lesssim E_{\textrm{n}}(0.05\,\textrm{fm}^{-3})\,(E_{\textrm{n}}(0.1\,\textrm{fm}^{-3}))\lesssim 7.3\,\textrm{MeV}\,(12.1\,\textrm{MeV})$, leading to a 68% (54%) reduction on the uncertainties on $E_{\textrm{n}}(\rho)$ at $0.05\,\textrm{fm}^{-3}\,(0.1\,\textrm{fm}^{-3})$, compared with the constraints from ChPT Tew13 ; Kru13 . These results indicate that combing the QCDSR and ChPT can significantly improve the predictions on the EOS of PNM.

In Fig. 40, the EOS of SNM $E_{0}(\rho)$ with different $\Phi$ and $g$ is shown. It is interesting to find that with a smaller $\Phi^{\prime}$ (or equivalently a larger $\Phi=\Phi^{\prime}/\langle\overline{q}q\rangle_{\textrm{vac}}$) the saturation properties are closer to the empirical constraints. For example, the saturation density (binding energy at that density) is found to be about $0.18\,\textrm{fm}^{-3}$ ($-$22.9 MeV) with $\Phi^{\prime}\approx 2.70$ and $g\approx-0.817$, introducing a relative improvement about 10% (14%) on $\rho_{0}$ ($E_{0}(\rho_{0})$) compared to the QCDSR-2 prediction ($0.2\,\textrm{fm}^{-3}$ and $-$26.7 MeV). On the other hand, with a larger $\Phi^{\prime}\approx 4.00$ and correspondingly $g\approx-0.552$, the saturation density (the corresponding binding energy) is found to about $0.24\,\textrm{fm}^{-3}$ ($-$34.0 MeV). The overall uncertainty on the $E_{0}(\rho)$ is plotted in the inset of Fig. 40. Combining the discussions given above, it shows that before the EOS of PNM and the density dependence of the chiral condensates at densities around the saturation density are well determined, it is hard to make accurate predictions on the EOS of SNM. Specifically, it could be found that while keeping the EOS of PNM to be consistent with the ChPT predictions within density from zero to about $0.18\,\textrm{fm}^{-3}$ and meanwhile the symmetry energy at $\rho_{\textrm{c}}$ fixed, the uncertainty on the nuclear saturation density (binding energy) due to the $\Phi$-term in Eq. (3.137), is about $0.06\,\textrm{fm}^{-3}$ ($-$11.1 MeV), a relative uncertainty about 38% (69%) compared to $\rho_{0}\approx 0.16\,\textrm{fm}^{-3}$ ($E_{0}(\rho_{0})\approx-16\,\textrm{MeV}$), strongly indicating that the higher order terms in density in the chiral condensates also have sizable impact on the EOS of SNM. Naturally, if other uncertainties are included, e.g., the uncertainties on $\sigma_{\textrm{N}}$, the twist-four four-quark condensates discussed in the next section, and the uncertainties on the symmetry energy at $\rho_{\textrm{c}}$, etc., it may lead the corresponding uncertainties on $\rho_{0}$ and $E_{0}(\rho_{0})$ much larger. A detailed analysis on this issue is beyond the main scope of the present work.

Finally, we show in Fig. 41 the density dependence of the chiral condensates both in SNM and PNM, including the uncertainties introduced by the parameters $\Phi$ and $g$. Besides the chiral symmetry restoration pattern discussed for u and d quarks in PNM in the above paragraphs (see Fig. 38), the nonlinear density corrections in the chiral condensates (3.137) also make the restoration of the chiral symmetry in SNM to occur at an even higher density (green band in the right panel of Fig. 41). More interestingly, the d quark chiral condensate in PNM with the allowable ranges of $\Phi$ and $g$ may even increase at a critical density about $0.25\,\textrm{fm}^{-3}$ (red band), indicating that the d quark in PNM is very stable. Since the main component of neutron stars is the neutrons (roughly 2/3 of the components of a neutron stars are d quarks), the hindrance of the chiral symmetry restoration of the d quark may have important consequences on investigating the quark-involved dynamical processes in these compact objects.

IX Twist-four Four-quark Condensates, Symmetry Energy and QCDSR-3

In this section, the effects of the twist-four four-quark six-dimensional condensates on the EOS of ANM are studied and correspondingly the QCDSR-3 is constructed. The twist-four condensates effects on the symmetry energy were first studied in ref. Jeo13 . Since the discussions in this section is very similar as those done in Section VIII, here we mainly focus on the related issues on the EOS of PNM, the EOS of SNM, and the symmetry energy.

The contributions to the QCDSR Eqs. (3.101) and (3.102) from the twist-four condensates for proton are given by Jeo13 ,

[TABLE]

the five parameters, i.e., $\Omega_{1}\sim\Omega_{5}$, characterizing the twist-four condensates together with three different parameter sets are given and discussed in ref. Jeo13 . The similar contributions can be written out for the neutron by exchanging the u and d quarks in the above expressions (i.e., the parameters $\Omega_{1}\sim\Omega_{5}$). In the following, we call them the set 1, set 2, and set 3 parameter sets, respectively.

In Fig. 42, we show the effects of the twist-four condensates on the EOS of SNM and the EOS of PNM, where the other parameters are the same as those obtained in QCDSR-2. It is clearly shown from Fig. 42 that the twist-four condensates have large impact both on the $E_{0}(\rho)$ and $E_{\textrm{n}}(\rho)$, as first pointed out in ref. Jeo13 . For instance, the $E_{0}(\rho)$ at $\rho_{\textrm{c}}$ takes value about 36.3 MeV (13.3 MeV, 0.3 MeV) in the twist-four condensate parameter set 1 (set 2, set 3) compared with the QCDSR-2 prediction about $-$18.9 MeV without these condensates, expanding an uncertainty of about 55.2 MeV.

Similarly, the $E_{\textrm{n}}(\rho)$ at $\rho_{\textrm{c}}$ takes value about 72.9 MeV (49.2 MeV, 35.4 MeV) in parameter set 1 (set 2, set 3) and 12.1 MeV in QCDSR-2 without the twist-four four-quark condensates, introducing an uncertainty of about 60.8 MeV. These sizable influence on the $E_{\textrm{n}}(\rho)$ and $E_{0}(\rho)$ is found to be consistent with the findings in ref. Jeo13 . However the symmetry energy obtained by the difference between the $E_{\textrm{n}}(\rho)$ and the $E_{0}(\rho)$ (i.e., the parabolic approximation) directly from the results shown in Fig. 42 is found to be about 36.6 MeV (35.9 MeV, 35.1 MeV) and 31.0 MeV in parameter set 1 (set 2, set 3) and in QCDSR-2 at $\rho_{\textrm{c}}$. The uncertainty on the symmetry energy from these four parameter sets is found to be maximally about 5.6 MeV at the cross density. The large uncertainties on the EOS of SNM and EOS of PNM due to the twist-four four-quark condensates are roughly canceled, leading to a relative smaller impact on the symmetry energy. In Fig. 43, the density dependence of the symmetry energy with the twist-four condensates included and that obtained in QCDSR-2 are shown through Eq. (2.22). Although one uses the QCDSR-2 parameters ($f$, $\Phi$ and $g$), the density dependence of the symmetry energy including the twist-four condensates is very close to the prediction by the QCDSR-2, within a wide range of densities. Specifically, the $E_{\textrm{sym}}(\rho_{\textrm{c}})$ is found to be about 31.4 MeV (30.9 MeV, 30.5 MeV) in parameter set 1 (set 2, set 3), indicating once again that the $E_{\textrm{HO}}(\rho)$ is already non-negligible at densities smaller than $\rho_{0}$ when compared with the parabolic approximation prediction just given above. It is also very interesting to notice from Fig. 43 that the twist-four condensates tend to soften the symmetry energy at densities larger than and/or around the saturation density. In the following we readjust the parameters $\Phi$, $g$ and $f$ according to the fitting scheme (Subsection III.6), i.e., the EOS of PNM at $\rho_{\textrm{vl}}$ to be about 4.2 MeV, the symmetry energy at $\rho_{\textrm{c}}$ to be about 26.65 MeV, and the meanwhile the $E_{\textrm{n}}(\rho)$ to be fitted with the APR EOS as much as possible.

Based on the fitting scheme (Subsection III.6), we obtain $f\approx 0.360$, $\Phi^{\prime}\approx 4.50$ and $g\approx-0.50$ when the twist-four four-quark condensates are included, and this parameter set is abbreviated as QCDSR-3. In the meanwhile the Ioffe parameter is still about $-$1.22, which is unaffected by the twist-four four-quark condensates. Shown in Fig. 44 is the EOS of PNM obtained in QCDSR-3, and the corresponding predictions from QCDSR-1, QCDSR-2 as well as the nonlinear RMF model with FSUGold are also included for comparison. For reference, we note the $E_{\textrm{n}}(\rho)$ at $\rho_{\textrm{c}}$ ($\rho_{0}$) in FSUGold is about 12.8 MeV (19.4 MeV). We find that the EOS of PNM from QCDSR-3 at densities $\lesssim 0.12\,\textrm{fm}^{-3}$ is essentially the same as the one without the twist-four condensates (i.e., in QCDSR-2), which is also in nice agreement with the FSUGold prediction at these densities. And at the nuclear saturation density $\rho_{0}=0.16\,\textrm{fm}^{-3}$, the $E_{\textrm{n}}(\rho_{0})$ changes from about 17.1 MeV in QCDSR-2 to 15.9 MeV in QCDSR-3 Cai17QCDSR , introducing a relative reduction about 7%. However, the discrepancy from the APR EOS becomes eventually apparent as densities increases $\gtrsim 0.12\,\textrm{fm}^{-3}$. Since the high-twist operators have some impacts on several processes in hadronic physics Col11 ; Gre07 , the exact knowledge on the density dependence of the $E_{\textrm{n}}(\rho)$ provides a novel tool to study these operators.

We show in Fig. 45 the density dependence of the symmetry energy obtained in QCDSR-3 as well as the corresponding predictions from QCDSR-1, QCDSR-2 and FSUGold. The symmetry energy at densities $\lesssim 0.1\,\textrm{fm}^{-3}$ is improved compared with several constraints, e.g., the $E_{\textrm{sym}}(\rho)$ at $0.04\,\textrm{fm}^{-3}$ is found to be about 11.9 MeV (compared with the 10.8 MeV in QCDSR-2), which is safely within the IAS constraints. Moreover, the symmetry energy at $\rho_{0}$ now is found to be about 31.6 MeV. Furthermore, the slope parameter of symmetry energy $L$ at $\rho_{\textrm{c}}$ is now about 49.8 MeV, which is consistent with the empirical constraints about $46\pm 4.5\,\textrm{MeV}$ Zha13 . The $\gamma$ parameter at $\rho_{\textrm{c}}$ now is about 0.62, which safely fall within the band in Fig. 34. On the other hand, the $L$ parameter at the saturation density $\rho_{0}$ in QCDSR-3 is found to be about 16.4 MeV. The softening of the symmetry energy at densities $\gtrsim\rho_{0}$ is a possible signal that even higher order terms in density in the condensates beyond the $\Phi$-term need to be included effectively.

Finally, we show the density dependence of $E_{0}(\rho)$ obtained in QCDSR-3 in Fig. 46. Obviously, the saturation properties obtained in QCDSR-3 and QCDSR-2 are similar. For instance, the saturation density (binding energy) in QCDSR-3 is now given as about 0.21 $\textrm{fm}^{-3}$ ($-$29.0 MeV), which is slightly larger (deeper) than the prediction by the QCDSR-2. It is also interesting to notice that the uncertainties introduced by the twist-four four-quark condensates on the saturation properties of the SNM could be covered by the effective ranges of $\Phi$ and $g$ discussed in the last section (see Fig. 40) although the origins of the uncertainties are different, showing that the effects of the twist-four condensates on the saturation properties of the SNM are smaller than the higher order density terms in the chiral condensates. This again shows the particular importance of the $\Phi$-term in the chiral condensates in Eq. (3.137). Moreover, other characteristic quantities for the EOS of ANM are also obtained in QCDSR-3, e.g., the nucleon Landau/Dirac effective mass $M_{\textrm{L/D}}^{\ast}(\rho_{0})/M$ in SNM at $\rho_{0}$ ($\rho_{\textrm{c}}$) is found to be about 1.10/0.67 (1.21/0.74), while the high order effects in the EOS of ANM $E_{\textrm{HO}}(\rho)$ is found to be about 10.9 MeV (6.1 MeV) at $\rho_{0}$ ($\rho_{\textrm{c}}$). These features together show again that the parameter set QCDSR-3 is very similar like the QCDSR-2.

X Summary and Outlook

In this work, we have systematically investigated the EOS of isospin asymmetric nucleonic matter within the framework of QCDSR, and mainly focused on the relativistic self-energy decomposition of the nuclear symmetry energy and the EOS of PNM. Based on the fitting scheme that the EOS of PNM at $\rho_{\textrm{vl}}=0.02\,\textrm{fm}^{-3}$ and the symmetry energy $E_{\textrm{sym}}(\rho)$ at $\rho_{\textrm{c}}=0.11\,\textrm{fm}^{-3}$ are fixed at 4.2 MeV and $26.65\pm 0.2$ MeV, respectively, and the total $E_{\textrm{n}}(\rho)$ to densities about $\rho_{0}$ is fitted to the APR EOS as much as possible, several interesting results are obtained as following:

1. In the conventional QCDSR, the prediction on the nucleon mass in vacuum is not necessarily about 939 MeV Coh95 . In this work the nucleon mass in vacuum is self-consistently determined via $M_{\textrm{D}}^{\ast}(0)\equiv M_{0}^{\ast}(0)=M\equiv 939\,\textrm{MeV}$, leading to the Ioffe parameter $t$ to be about $-1.22$ (which is close to its natural value $t_{\textrm{Ioffe}}=-1$) in three QCDSR parameter sets, and more specifically the Ioffe parameter expressed in terms of the chiral condensate and the gluon condensate in vacuum is obtained, see Eq. (IV). This paves an important step to the consequent investigations on the EOS. For instance, the EOS of SNM via $E_{0}(\rho)=\rho^{-1}\int_{0}^{\rho}\mathrm{d}\rho[e_{\textrm{F}}^{\ast}(\rho)+\Sigma_{\textrm{V}}^{0}(\rho)]-M$ needs an accurate nucleon mass in vacuum.

2. The Lorentz structure based on the nucleon self-energy decomposition of the symmetry energy is carefully explored. Specifically, the first-order symmetry scalar self-energy is found to depend heavily on the nucleon sigma term $\sigma_{\textrm{N}}$. For instance, in the msQCDSR, this term is given by $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)=2\pi^{2}M_{0}^{\ast}\sigma_{\textrm{N}}\alpha\rho/(\mathscr{M}^{2}m_{\textrm{q}}e_{\textrm{F}}^{\ast}\beta)$, where $\alpha/\beta$ characterizes the isospin effects of the chiral condensates at linear order. This relation establishes a useful connection between the symmetry energy and the evolution of the quark mass since the nucleon sigma term actually characterizes the evolution of the nucleon mass as a function of the light quark mass, i.e., $\sigma_{\textrm{N}}=m_{\textrm{q}}\mathrm{d}M/\mathrm{d}m_{\textrm{q}}$. Consequently, the $\sigma_{\textrm{N}}$ term is found to largely affect $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ since $y=M_{0}^{\ast}/e_{\textrm{F}}^{\ast}\approx 1-x^{2}/2$ with $x=k_{\textrm{F}}/M_{0}^{\ast}$, where effective mass $M_{0}^{\ast}=M+\Sigma_{\textrm{S}}^{0}$ depends on the $\sigma_{\textrm{N}}$ almost linearly, and the final effect is as $\sigma_{\textrm{N}}$ increases the $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$ decreases. Moreover, the kinetic nuclear symmetry energy $E_{\textrm{sym}}^{\textrm{kin}}(\rho)=k_{\textrm{F}}^{2}/6e_{\textrm{F}}^{\ast}$ is also found to be largely affected by the $\sigma_{\textrm{N}}$ since $e_{\textrm{F}}^{\ast}=(M_{0}^{\ast,2}+k_{\textrm{F}}^{2})^{1/2}$ could be approximated as $M_{0}^{\ast}+k_{\textrm{F}}^{2}/2M_{0}^{\ast}=M+\Sigma_{\textrm{S}}^{0}+k_{\textrm{F}}^{2}/2M_{0}^{\ast}$, i.e, $E_{\textrm{sym}}^{\textrm{kin}}(\rho)$ increases as $\sigma_{\textrm{N}}$ increases. Finally, as the $\sigma_{\textrm{N}}$ increases, the total symmetry energy $E_{\textrm{sym}}(\rho)$ increases as a result of the competition between the $E_{\textrm{sym}}^{\textrm{kin}}(\rho)$ and the $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$, see the left panel of Fig. 21.

3. The vector nucleon self-energy contributes to the nuclear symmetry energy as, e.g., $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)=4\pi^{2}\rho/\mathscr{M}^{2}$, in the msQCDSR, and this part is originated from the densities of d quarks ($\langle\textrm{d}^{{\dagger}}\textrm{d}\rangle_{\rho,\delta}=(3+\delta)\rho/2$) and u quarks ($\langle\textrm{u}^{{\dagger}}\textrm{u}\rangle_{\rho,\delta}=(3-\delta)\rho/2$). The linear density dependence is almost unchanged in the fQCDSR via the very complicated QCDSR equations. Moreover, the $E_{\textrm{sym}}^{\textrm{1st,V}}(\rho)$ becomes the dominant contribution to the symmetry energy as $\rho\gtrsim 0.1\,\textrm{fm}^{-3}$.

4. Contributions to the nuclear symmetry energy due to the momentum dependence of the nucleon self-energies are found to be $E_{\textrm{sym}}^{\textrm{mom,0,S}}(\rho)\sim[\mathrm{d}\Sigma_{\textrm{S}}^{0}/\mathrm{d}|\textbf{k}|]_{k_{\textrm{F}}}<0$ and $E_{\textrm{sym}}^{\textrm{mom,0,V}}(\rho)\sim[\mathrm{d}\Sigma_{\textrm{V}}^{0}/\mathrm{d}|\textbf{k}|]_{k_{\textrm{F}}}>0$, at densities smaller or around $\rho_{0}$. Moreover, the magnitude of these terms is much smaller than the kinetic part and the first order symmetry parts, showing the corresponding momentum dependence in QCDSR is weak.

5. The dependence of the symmetry energy on the strange quark mass is found to be large. This large effect can be understood via $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)=2\pi^{2}M_{0}^{\ast}\sigma_{\textrm{N}}\alpha\rho/(\mathscr{M}^{2}m_{\textrm{q}}e_{\textrm{F}}^{\ast}\beta)$ in the msQCDSR since a larger strange quark mass induces a smaller factor $\alpha/\beta$, and consequently a smaller $E_{\textrm{sym}}^{\textrm{1st,S}}(\rho)$. The connection between the strange quark mass and the symmetry energy provides a useful bridge to understand the possible origins of the uncertainties on the $E_{\textrm{sym}}(\rho)$, see the right panel of Fig. 21 and Fig. 22.

6. A few useful approximations for the EOS of ANM and the related quantities are obtained in the QCDSR, and among which the approximation for the EOS of PNM, i.e., $E_{\textrm{n}}(\rho)\approx E_{\textrm{n}}^{\textrm{FFG}}(\rho)+({M\rho}/2{\langle\overline{q}q\rangle_{\textrm{vac}}})[(1-\xi)({\sigma_{\textrm{N}}}/{2m_{\textrm{q}}})-5]$, already has quantitative predictive power at low densities, where $E_{\textrm{n}}^{\textrm{FFG}}(\rho)$ is the FFG EOS of PNM. This low density formula clearly demonstrates how the chiral condensate goes into play in the EOS of PNM, e.g., the second term characterized by several constants ($\xi,\sigma_{\textrm{N}},m_{\textrm{q}}$ and $\langle\overline{q}q\rangle_{\textrm{vac}}$) leads to a reduction on the $E_{\textrm{n}}(\rho)$ compared to the FFG prediction.

7. The high order term in the EOS of ANM characterized by the $E_{\textrm{HO}}(\rho)\approx E_{\textrm{sym,4}}(\rho)+E_{\textrm{sym,6}}(\rho)+\cdots$ is carefully studied in the QCDSR. Specifically, $E_{\textrm{HO}}(\rho)$ is found to be sizable at densities $\lesssim\rho_{0}$, indicating the conventional parabolic approximation for the EOS of ANM is broken in QCDSR. For instance, the $E_{\textrm{HO}}(\rho)$ at the saturation density in the fQCDSR calculations with the fitting scheme adopted is found to be about 1.1$\sim$11.9 MeV, indicating the uncertainty on $E_{\textrm{HO}}(\rho)$ in QCDSR is essentially large compared with the one from phenomenological models.

8. The correlation between the symmetry energy and several quantities characterizing the quark/gluon condensates is investigated. Specifically, besides the strong dependence on the chiral condensates in vacuum $\langle\overline{q}q\rangle_{\textrm{vac}}$, the symmetry energy is also found to heavily depend on the five-dimensional mixing condensate $\langle g_{\textrm{s}}q^{{\dagger}}\sigma\mathcal{G}q\rangle_{\rho,\delta}$. This correlation provides a novel tool to explore the properties of symmetry energy and even the EOS of ANM via the knowledge on the in-medium quark/gluon condensates from, e.g., hadronic physics.

9. The effects of twist-four four-quark six-dimensional condensates on the EOS of SNM and the EOS of PNM are found to be large. For instance, these condensates induce an amount about 50 MeV to 60 MeV on the $E_{0}(\rho)$ and $E_{\textrm{n}}(\rho)$. However, their effects on the symmetry energy are almost canceled since the $E_{\textrm{sym}}(\rho)$ is roughly the difference between the $E_{\textrm{n}}(\rho)$ and $E_{0}(\rho)$, and as a result the twist-four condensates induces an uncertainty of about several MeVs on the symmetry energy.

10. The effective higher order terms in density in the chiral condensates, i.e., $\Phi(1\mp g\delta)\rho^{2}$, are found to strongly affect the EOS of PNM and the EOS of SNM. By refitting the model parameters in the presence of the $\Phi$-term, the EOS of PNM at densities around $\rho_{0}$ is found to be systematically consistent with the APR EOS, which is selected as the reference EOS in this work. Moreover, the higher-order density terms in quark condensates also leads to the stabilization of u/d chiral condensates at higher densities, which may have important implications on the QCD phase diagram under extreme conditions at low temperatures, large isospin and large baryon chemical potentials, which is essential for understanding the physical degrees of freedom in the core of neutron stars. Furthermore, using a different higher order density term such as $\Phi(1\mp g\delta)\rho^{5/3}$ gives a very similar prediction on the $E_{\textrm{n}}(\rho)$, see Fig. 31.

11. Three parameter sets of QCDSR are constructed, i.e., QCDSR-1 (naive QCDSR), QCDSR-2 and QCDSR-3, respectively. The QCDSR-1 includes only the linear approximation for the chiral condensates without the effective $\Phi$-term in the chiral condensates (3.137) and the twist-four four-quark condensates. Compared with the QCDSR-1, the QCDSR-2 additionally includes the effective $\Phi$-term in the chiral condensates (3.137) but without the twist-four four-quark condensates. The QCDSR-3 includes the linear approximation for the chiral condensates, the effective $\Phi$-term in the chiral condensates (3.137) and the twist-four four-quark condensates. The symmetry energy at $\rho_{0}\approx 0.16\,\textrm{fm}^{-3}$ in QCDSR-1, QCDSR-2 and QCDSR-3 is found to be about 45.1 MeV, 35.2 MeV and 31.6 MeV, respectively, while the corresponding slope parameter of the symmetry energy $L$ at $\rho_{\textrm{c}}$ ($\rho_{0}$) is found to be 105.9 (196.8) MeV, 64.7 (67.5) MeV, and 49.8 (16.4) MeV, respectively. The tendency of the change in $E_{\textrm{sym}}(\rho)$ and $L(\rho)$ shows that the higher order terms in density in the chiral condensates improve the density behavior of the symmetry energy compared with the empirical constraints, while the twist-four condensates soften the $E_{\textrm{sym}}(\rho)$ at densities larger than about $\rho_{\textrm{c}}$. Moreover, the saturation properties of SNM are largely improved from QCDSR-1 to QCDSR-2 or QCDSR-3, e.g., the ($\rho_{0},E_{0}(\rho_{0})$) are changed from $(0.6\,\textrm{fm}^{-3},-100.0\,\textrm{MeV})$ in QCDSR-1 to $(0.2\,\textrm{fm}^{-3},-26.7\,\textrm{MeV})$ in QCDSR-2 or $(0.21\,\textrm{fm}^{-3},-29.0\,\textrm{fm}^{-3})$ in QCDSR-3. Furthermore, the EOS of PNM obtained in three QCDSR parameter sets is consistent with each other at low densities less than about $0.08\,\textrm{fm}^{-3}$, indicating that at these low densities the naive QCDSR is well-behaved for the $E_{\textrm{n}}(\rho)$. Finally, we have better constrained the $E_{\textrm{n}}(\rho)$ in the density region from about $0.04\,\textrm{fm}^{-3}$ to $0.12\,\textrm{fm}^{-3}$ by combining the results from QCDSR and ChPT, e.g., the $E_{\textrm{n}}(0.05\,\textrm{fm}^{-3})\,(E_{\textrm{n}}(0.1\,\textrm{fm}^{-3}))$ is constrained to be between 6.8 MeV (10.0 MeV) and 7.3 MeV (12.1 MeV), leading to an uncertainty about 0.5 MeV and 2.1 MeV, respectively.

Besides the above results we have obtained from the QCDSR, a few interesting issues that are closely related to our present work should be pointed out and need further exploration in the future QCDSR calculations on dense nucleonic matter, i.e.,

1. In this work as well as in many conventional QCDSR calculations Coh95 , the four-quark condensates effects are incorporated by the effective parameter $f$, i.e., using $(1-f)\langle\overline{q}q\rangle_{\textrm{vac}}^{2}+f\langle\overline{q}q\rangle_{\rho,\delta}^{2}$ to account for the four-quark condensates at finite densities. In this work, the $f$ is found to largely influence several quantities, such as the EOS of PNM and the $E_{0}(\rho)$, and the value of $f$ is essentially determined by the $E_{\textrm{n}}(\rho)$ at the very low density $\rho_{\textrm{vl}}\approx 0.02\,\textrm{fm}^{-3}$ in the present work. From the more fundamental viewpoint, it is important to explore the density behavior of the four-quark condensates, in order to make further progress in applying QCDSR to dense nucleonic matter calculations. For instance, in ref. Tho07 , more phenomenological parameters are introduced into the QCDSR equations, and they are determined by nuclear quantities and/or other information from, e.g., hadronic physics.

2. The three-body forces (TBF) are found to be important for the saturation properties of SNM, e.g., in Brueckner–Hartree–Fock (BHF) calculations (see, e.g., ref. ZHLi06 ) and in phenomenological approach. For instance, in the SHF model, a traditional two-body force contributes a term proportional to $\rho$ to the EOS, and a $\rho^{1+\alpha}$ term emerges once the effective three-body force is considered Zha16xx with $\alpha$ the parameter characterizing the three-body force. Recently, three-body forces even the four-body forces are included in the QCDSR calculations for nucleonic matter Dru17a ; Dru17b ; Dru19 , and it is really also interesting and important to see how these many-body forces influence, for example, the $E_{\textrm{n}}(\rho)$ and/or the nuclear symmetry energy. Investigations on these problems will help us better understanding the origins of the uncertainties on the symmetry energy and/or the difficulties to produce reasonable saturation properties of the symmetric matter, and they are also important for making further progresses in the nucleonic matter calculations, such as to explore the incompressibility property of the ANM with any isospin asymmetry, and/or the single nucleon optical potential Ham90 ; Hol16xx ; LiX13 ; LiX15 . It is also interesting to investigating the three-body force and its connection to the higher order density terms in the chiral condensates.

3. The Borel transformation with the Borel mass $\mathscr{M}$ as a real parameter is the standard treatment in QCDSR to deal with the high order states including the continuum excitations. Recently, a generalization to the complex-valued Borel mass $\mathscr{M}$ was introduced in ref. Ara14 (see also ref. Iof01 ). It was demonstrated that the complex-valued sum rules approach allows one to extract the spectral function with a significantly improved resolution, and thus provides a useful tool to study more detailed structures of the hadronic spectrum Ara14 . To our purpose, it would be interesting to investigate whether the complex-valued Borel transformation could improve the calculations on the EOS of dense nucleonic matter in the high density region where the high mass-dimensional condensates and continuum effects are important.

4. Finally, the neutron matter at sub-saturation even to very low densities composed of spin-down and -up neutrons with a large s-wave scattering length shows several universal properties Tan08 , such as the simplicity of its EOS characterized by a few universal parameters CaiLi15a ; Sch05 ; Kru15 ; Kol16 ; ZhaNB17 . Moreover, the high momentum tail above the Fermi surface of the single nucleon momentum distribution function in cold PNM is also found to be very similar as that in ultra-cold atomic Fermi gases Hen14 although the magnitude of the density in the two systems differs by an amount about 25 orders Hen15 . Naturally, the cold PNM at low densities provides a perfect testing bed to explore novel ideas in the unitary region Gio08 ; Blo08 , helping to find deep physical principles behind these quantum many-body systems in this so-called unitary region Zwe12 . Recently, the sum rule approach with the help of maximum entropy (ME) method was applied to investigate the imaginary part of the particle self-energy in the unitary Fermi gas Gub15 ; Gub18 . Thus it is interesting to explore, e.g., the imaginary part of the neutron self-energies in PNM under the QCDSR+ME method, which will be extremely useful for exploring the transport properties of the PNM, or even to generalize the method to a general isospin asymmetric nucleonic matter to better understanding the quantum many-body properties of the system.

Our results in the present work have demonstrated that the QCDSR approach can be used to explore the properties of ANM in a quantitative manner, at least in lower density region. The QCDSR approach establishes a bridge connecting the EOS of ANM and the non-perturbative QCD vacuum, and thus provides a useful way to understand the properties of dense nucleonic matter from non-perturbative QCD vacuum. These studies are helpful to investigate the QCD origins about the uncertainties of nucleonic matter properties, e.g., the uncertainties of the symmetry energy. On the other hand, the exact knowledge on the EOS of ANM extracted from experiments, observations and model-independent calculations is also very useful for understanding the quark/gluon condensates in nuclear medium, which can provide important information on the chiral symmetry restoration phase transition in nuclear matter as well as the in-medium effects of hadron properties.

ACKNOWLEDGEMENTS

This work was supported in part by the National Natural Science Foundation of China under Grant Nos. 11625521 and 11805118, the Major State Basic Research Development Program (973 Program) in China under Contract No. 2015CB856904, the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning, Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education, China, and the Science and Technology Commission of Shanghai Municipality (11DZ2260700).