# Embedding short-range correlations in relativistic density functionals
through quasi-deuterons

S. Burrello<EMAIL_ADDRESS>ORCID: 0000-0002-1132-4073 Technische
Universität Darmstadt, Fachbereich Physik, Institut für Kernphysik,
Schlossgartenstraße 9, D-64289 Darmstadt, Germany  S. Typel<EMAIL_ADDRESS>darmstadt.de, ORCID: 0000-0003-3238-9973 Technische Universität Darmstadt,
Fachbereich Physik, Institut für Kernphysik, Schlossgartenstraße 9, D-64289
Darmstadt, Germany  GSI Helmholtzzentrum für Schwerionenforschung GmbH,
Theorie, Planckstraße 1, D-64291 Darmstadt, Germany

(August 27, 2024; August 27, 2024)

###### Abstract

Background

The formation of clusters at sub-saturation densities, as a result of many-
body correlations, constitutes an essential feature for a reliable
modelization of the nuclear matter equation of state (EoS). Phenomenological
models that make use of energy density functionals (EDFs) offer a convenient
approach to account for the presence of these bound states of nucleons when
introduced as additional degrees of freedom. However, in these models clusters
dissolve, by construction, when the nuclear saturation density is approached
from below, revealing inconsistencies with recent findings that evidence the
existence of short-range correlations (SRCs) even at larger densities.

Purpose

The idea of this work is to incorporate SRCs in established models for the
EoS, in light of the importance of these features for the description of
heavy-ion collisions, nuclear structure and in the astrophysical context. Our
aim is to describe SRCs at supra-saturation densities by using effective
quasi-clusters immersed in dense matter as a surrogate for correlations, in a
regime where cluster dissolution is usually predicted in phenomenological
models.

Method

Within the EDF framework, we explore a novel approach to embed SRCs within a
relativistic mean-field model with density dependent couplings through the
introduction of suitable in-medium modifications of the cluster properties, in
particular their binding energy shifts, which are responsible for describing
the cluster dissolution. As a first exploratory step, the example of a quasi-
deuteron within the generalized relativistic density functional approach is
investigated. The zero temperature case is examined, where the deuteron
fraction is given by the density of a boson condensate.

Results

For the first time, suitable parameterizations of the cluster mass shift at
zero temperature are derived for all baryon densities. They are constrained by
experimental results for the effective deuteron fraction in nuclear matter
near saturation and by microscopic many-body calculations in the low-density
limit. A proper description of well-constrained nuclear matter quantities at
saturation is kept through a refit of the nucleon meson coupling strengths.
The proposed parameterizations allow to also determine the density dependence
of the quasi-deuteron mass fraction at arbitrary isospin asymmetries. The
strength of the deuteron-meson couplings is assessed to be of crucial
importance. Novel effects on some thermodynamic quantities, such as the matter
incompressibility, the symmetry energy and its slope, are finally discerned
and discussed.

Conclusions

The findings of the present study represent a first step to improve the
description of nuclear matter and its EoS at supra-saturation densities in
EDFs by considering correlations in an effective way. In a next step, the
single-particle momentum distributions in nuclear matter can be explored using
proper wave functions of the quasi-deuteron in the medium. The momentum
distributions is expected to exhibit a high-momentum tail, as observed in the
experimental study of SRCs by nucleon knockout with high-energy electrons.

##  I Introduction

The Equation of State (EoS) of strongly interacting matter is a fundamental
ingredient in the theoretical description of compact stars. It plays also a
crucial role in many astrophysical simulations of, e.g., core-collapse
supernovae and neutron star mergers Shen  _et al._ (2011); Hempel  _et al._
(2012); Lattimer and Prakash (2016); Abbott  _et al._ (2018). Tight links have
been also established between the properties of the EoS and those of finite
nuclei, concerning both their structure and reaction dynamics Klähn  _et al._
(2006); Zheng  _et al._ (2017, 2018). A large variety of approaches has been
employed therefore in the last decades to develop reliable models for the EoS,
which are constrained by both nuclear physics experiments and astronomical
observations. For astrophysical simulations, tables of global, multi-purpose
EoSs are in particular needed, being provided by sophisticated theoretical
methods available in literature Oertel  _et al._ (2017).

One class is given by microscopic ab-initio models, which try to solve the
nuclear many-body problem using advanced methods and realistic interactions
that are constrained by scattering data and properties of few-nucleon systems,
see, e.g., Refs. Akmal  _et al._ (1998); Carlson  _et al._ (2015) Other
approaches based on the effective-field theory exploit systematic developments
from quantum chromodynamics and symmetry concepts, e.g. in chiral perturbation
theory, providing also uncertainty estimates Epelbaum  _et al._ (2009); Hammer
_et al._ (2020); Furnstahl  _et al._ (2021). However, these models fail in
properly describing the formation of clusters at densities below nuclear
saturation, where these bound states of nucleons emerge as many-body
correlations generated by the short-range nucleon-nucleon interaction. Such
low-density conditions are encountered in various systems: the debris of
heavy-ion collisions (HIC), the post-bounce evolution of core-collapse
supernovae, and the surface of nuclei Typel (2014); Tanaka  _et al._ (2021)
where the formation of clusters is a prerequisite for cluster radioactivity
and, in particular, $\alpha$-decay of heavy nuclei. The emergence of clusters
is definitely an essential feature for the modelization of a realistic EoS
Raduta and Gulminelli (2010).

Phenomenological models, whose parameters are constrained from experimental
results of HIC and astronomical observations and/or by directly fitting
properties of finite nuclei and nuclear matter near saturation, offer a
convenient alternative to microscopic models to approach this problem. A
widely used class of phenomenological approaches is based on energy density
functionals (EDFs), which are usually derived in the self-consistent mean-
field approximation with an effective in-medium interaction without a direct
connection to the nucleon-nucleon interaction in free space Bender  _et al._
(2003). There are various versions of this approach, e.g., non-relativistic
models using Skyrme or Gogny type interactions Dutra  _et al._ (2012);
Sellahewa and Rios (2014) or relativistic models based on the exchange of
mesons Dutra  _et al._ (2014). In recent years, several attempts have been
made to also directly link EDFs to microscopic ingredients Furnstahl, R. J.
(2020); Grasso (2019); Marino  _et al._ (2021). For example, a special class
of functionals inspired by effective field theories (EFTs) and bench-marked on
ab-initio predictions have been designed Yang  _et al._ (2016); Grasso  _et
al._ (2017). They were applied to finite nuclei Burrello  _et al._ (2021) and
finite temperature EoS calculations of pure neutron matter (PNM) Burrello and
Grasso (2022). Moreover, few steps were made towards the construction of a
power counting in EDF Yang  _et al._ (2017); Burrello  _et al._ (2020).

However, EDFs derived from phenomenological mean-field models fail as well,
when only nucleons are considered as basic constituents. Further progress is
only achieved if clusters are introduced as additional explicit degrees of
freedom at low densities Typel  _et al._ (2010). The dilute matter is then
depicted as an ideal mixture of nucleons and all nuclei from the table of
isotopes in thermodynamic equilibrium. Such a model is called nuclear
statistical equilibrium (NSE) and is widely used in the astrophysical context,
where it leads to a reliable description for the chemical composition of
stellar matter at sub-saturation densities Raduta and Gulminelli (2010);
Burrello  _et al._ (2015). However, models like NSE are considered valid only
as long as the interaction between the constituents can be neglected. Thus
they fail at higher densities where in-medium effects become important,
leading to the dissolution of clusters and the transition to cluster-free
nuclear matter. The dissolution of a cluster, i.e., the so-called Mott effect,
which is expected when approaching nuclear saturation density from below, can
be produced through the introduction of an excluded-volume mechanism, which is
just a simple geometric concept Raduta and Gulminelli (2010); Hempel and
Schaffner-Bielich (2010); Hempel  _et al._ (2011); Sagun  _et al._ (2014).
More microscopically, the formation and dissolution of light clusters in
nuclear matter can be treated using a quantum statistical approach with
thermodynamic Green’s functions, see, e.g., Refs. Röpke  _et al._ (1982a, b).

An alternative scenario was proposed in the last decade, when interacting
clusters were introduced as explicit degrees in relativistic density
functionals, firstly concentrating on light hydrogen and helium isotopes,
whose properties are modified in the medium Typel  _et al._ (2010). Contrary
to the chemical picture traditionally adopted in NSE-like models, where the
properties of the correlated states of nucleons are assumed to be independent
of the medium, in such a physical picture, the in-medium effects are addressed
by introducing a proper modification of the masses of the clusters, as
inspired by the quantum statistical approach. Following this idea, the
effective binding energies are expected to increase with density Typel (2020).

Many-nucleon correlations in the continuum, which still survive above the Mott
density, are then described in terms of effective resonances or quasi-
clusters. In current phenomenological approaches with the excluded-volume
mechanism or with medium-dependent mass shifts, these states are statistically
suppressed by construction beyond saturation, so that only nucleons should
remain as independent quasi-particles Pais and Typel (2017a). Mean-field type
descriptions of nuclear matter above saturation consider thus the system as a
free Fermi gas, with the usual step function in the single-particle momentum
distribution at zero temperature. Such a picture is however inconsistent with
recent experimental results from nucleon knock-out reactions on nuclei using
inelastic electron scattering Fomin  _et al._ (2012); Hen  _et al._ (2017).
These studies clearly evidence the smearing of the nucleon Fermi surface and
the emergence a high-momentum tail (HMT) in the single-nucleon momentum
distribution of cold nucleonic matter, ascribable to the existence of sizeable
nucleon-nucleon short-range correlations (SRCs) even at saturation density
Arrington  _et al._ (2012); degli Atti (2015); Hen  _et al._ (2017).

Experimental investigations assessed that SRCs pairs are formed by
approximately 20% of nucleons in various measured nuclei Egiyan  _et al._
(2006); Subedi  _et al._ (2008); Hen  _et al._ (2014). They are characterized
by large relative and small center-of-mass (c.m.) momenta. Moreover, some
results concluded that their magnitude is spin- and isospin-dependent, with a
clear dominance in the neutron-proton channel Subedi  _et al._ (2008); Hen
_et al._ (2014); Duer  _et al._ (2018, 2019), that affects the ratio of
minority and majority species in asymmetric nuclear matter, in the bulk part
and tail of the single-particle momentum distributions Duer  _et al._ (2018);
Schmookler  _et al._ (2019). Extrapolating these experimental results from
finite nuclei to infinite nuclear matter, useful information around saturation
density was then deduced. At higher densities, only numerical investigations
exist to predict the density dependence of these SRC pairs. From the analysis,
the probability for nucleons to form SRC pairs seems to have a minimum in the
neighbourhood of the saturation density, owing to the interplay between the
tensor component and the repulsive core of the nuclear force Rios  _et al._
(2014); Yang  _et al._ (2019); Li and Schulze (2016).

SRCs may also change the balance between kinetic and potential contributions
to the energy of the system. Thus their introduction is expected to
significantly affect the controversial density dependence, in particular at
supra-saturation densities, of the nuclear symmetry energy Cai and Li (2016);
Li  _et al._ (2018), which quantifies the difference between the total energy
of PNM and symmetric nuclear matter (SNM). Many theoretical and experimental
investigations are currently investigating the density dependence of this
quantity Li  _et al._ (2008); Russotto  _et al._ (2016); Zheng  _et al._
(2016); Burrello  _et al._ (2019a, b).

The purpose of the present work is to incorporate SRCs in established models
for the EoS, in light of the key-importance of these features for the
description of HIC, nuclear structure and in the astrophysical context
Weinstein  _et al._ (2011); Van Cuyck  _et al._ (2016); Guo  _et al._ (2021);
Lu  _et al._ (2021a). Our aim is to explicitly treat SRCs at supra-saturation
densities by using effective quasi-clusters immersed in dense matter as a
surrogate for correlations, through the introduction of proper in-medium
modifications of the cluster mass shifts. In this regime cluster dissolution
is usually predicted in actual realizations of phenomenological models. Within
the EDF framework, we propose thus a novel approach to embed the SRCs within a
relativistic mean-field model (RMF) with density dependent coupling Typel and
Wolter (1999) through a substantial modification of the cluster mass shift at
high densities.

Given its phenomenological nature, the adopted approach does not allow to
investigate the origin of the SRCs. Several effects coexist in the high-
density regime that can lead to a smearing of the single-particle distribution
functions and the appearance of HMTs as known from the description of heavy-
ion collisions in transport theory Buss  _et al._ (2012). They are due to
different (repulsive) components of the nucleon-nucleon interaction and can be
hardly disentangled. A deeper insight on these features would indeed require a
more microscopic treatment of these SRCs, going beyond the scope of the
present work.

As a first exploratory step, the example of a quasi-deuteron within this
generalized relativistic density functional (GRDF) is currently explored,
since two-body SRCs in the neutron-proton ${}^{3}S_{1}$ channel are much more
important than other many-body correlations. The zero temperature case is
examined, where the deuteron is represented by a boson condensate that
determines the mass fraction and leads to specific conditions to the
parameterization of the mass shift. The purpose of this work is then to
propose possible mass shift parameterizations that will be employed to
determine the density dependence of the quasi-deuteron mass fraction at
arbitrary isospin asymmetries. The final ambitious goal is to investigate the
effect of accounting for the SRCs in an effective way on the EoS, and some
related thermodynamic quantities, at supra-saturation densities. In this
context, it is worthwhile to mention that a recent study employed the concept
of a mass shift in a similar approach to model the effective interaction of a
possible heavy particle with baryon number $B=2$, the so-called sexaquark, in
nuclear matter. It is treated as a boson condensate like the deuteron and
affects the EoS of compact-star matter and thus the properties of neutron
stars Shahrbaf  _et al._ (2022).

The manuscript is structured as follows. In Section II the theoretical
formalism is illustrated and the fundamental principles and basic formulas of
the GRDF are given. The case of zero temperature is studied, where quasi-
deuteron condensation is expected, and the role of the deuteron-meson coupling
strength is discussed. Section III explores the different constraints for the
deuteron mass shifts. Section IV concentrates on the mass shift
parameterization, as suitably derived for nuclear matter at zero temperature.
Then, the corresponding density dependence of the quasi-deuteron mass fraction
is obtained. The impact on the EoS and on some general properties of nuclear
matter at arbitrary neutron-proton asymmetries is shown in Section V.
Conclusions and an outlook are finally given in Section VI. Details on the
formal derivation of some quantities, on the conversion of parameters and
analytical expressions for the mass shift parameters are furthermore collected
in four appendices.

##  II Theoretical formalism

###  II.1 Generalized Relativistic Density Functional

The GRDF is a density functional derived from a RMF with nucleons and further
degrees of freedom. Their effective interaction in the medium is described by
the exchange of mesons with density dependent couplings Typel and Wolter
(1999); Typel  _et al._ (2010). In such a model, light clusters are explicitly
introduced, allowing for a unified treatment from matter with bound nucleons
at low densities to matter possibly made of only neutrons ($n$) and protons
($p$) at high densities. All degrees of freedom are represented by quasi-
particles with self-energies that incorporate the effects of the interaction.
For sake of simplicity, only 2H nuclei (labeled as $d$ in the following) are
added to nucleons as degrees of freedom in light of their expected prevalent
importance discussed in Section I.

Let us thus consider the general case of asymmetric nuclear matter (ANM)
composed of neutrons, protons and deuterons with particle number densities
$n_{i}$ ($i=n,p,d$), baryon numbers $A_{i}$ and charge numbers $Z_{i}$. The
system is usually characterized by specifying the baryon density
$n_{b}=n_{n}+n_{p}+2n_{d}$, the isospin asymmetry $\beta=(n_{n}-n_{p})/n_{b}$
and the temperature $T$. In this section, the basic formulas of the
theoretical formalism are given for the most general case of finite
temperature. However, the zero temperature case will be considered in the
analysis performed in the following sections, leaving the analysis at finite
$T$ for future work. In this work, the masses of neutrons and protons will be
taken as equal to the average nucleon mass $m_{\mathrm{nuc}}$ so that
$m_{n}=m_{p}=m_{\mathrm{nuc}}$. The same values as given in Ref. Typel  _et
al._ (2010) are considered.

Following the framework illustrated in Ref. Pais and Typel (2017b, a), the
thermodynamic properties of nuclear matter are completely determined once the
grand canonical potential density $\tilde{\omega}(T,\\{\mu_{i}\\})$ is
specified. It depends, apart from the temperature, on the chemical potentials
$\mu_{i}$ of all constituents $i$. In the present work, three types of mesons
are considered: an isoscalar scalar $\sigma$ meson to describe the attraction
between nucleons, an isoscalar vector $\omega$ meson for their repulsion and
an isovector vector $\rho$ meson for the isospin dependence of the strong
force. The interaction between a baryon $i$ and a meson $j$
($j=\sigma,\omega,\rho$) is realized by a minimal coupling with a strength
that is given by the product of a scaling factor $\chi_{ij}$, the mass number
$A_{i}$, and a coupling $\Gamma_{j}$. The latter quantity depends on the
baryon density $n_{b}$ to describe the medium dependence of the effective
interaction. Different prescriptions might be adopted as recently discussed in
Ref. Typel and Alvear Terrero (2020). In this work the functional form of the
couplings as introduced in Ref. Typel and Wolter (1999) is used. In the
application of the GRDF to homogeneous nuclear matter, only the ratio
$\Gamma_{j}/m_{j}$ of coupling and mass $m_{j}$ of the mesons $j$ is relevant.
Hence it is convenient to introduce the coefficients

$C_{j}=\frac{\Gamma_{j}^{2}}{m_{j}^{2}}$ (1)

and their derivatives

$C_{j}^{\prime}=\frac{dC_{j}}{dn_{b}}=2\frac{\Gamma_{j}}{m_{j}^{2}}\frac{d\Gamma_{j}}{dn_{b}}$
(2)

with respect to the baryon density.

The total grand canonical potential density of the system can be written as

$\displaystyle\tilde{\omega}(T,\\{\mu_{i}\\})=$
$\displaystyle\sum_{i}\tilde{\omega}_{i}+\tilde{\omega}_{d}^{(\mathrm{cond})}+\tilde{\omega}_{\mathrm{meson}}-\tilde{\omega}_{\mathrm{meson}}^{(\mathrm{r})}-\tilde{\omega}_{\mathrm{mass}}^{(\mathrm{r})}$

containing the standard expression for the single quasi-particle (non-mesonic)
contribution

$\tilde{\omega}_{i}=-T\dfrac{g_{i}}{\sigma_{i}}\int\dfrac{d^{3}k}{(2\pi)^{3}}\ln\left[1+\sigma_{i}\exp\left(-\dfrac{E_{i}-\mu_{i}^{\ast}}{T}\right)\right]$
(4)

with the well-known integral over the momentum $k$ that appears in the quasi-
particle energy

$E_{i}=\sqrt{k^{2}+\left(m_{i}^{\ast}\right)^{2}}$ (5)

assuming natural units such that $\hbar=c=1$. In Eq. (4), the sign factor
$\sigma_{i}$ distinguishes the particle statistics ($\sigma_{i}=1$ for
fermions and $\sigma_{i}=-1$ for bosons, respectively) and $g_{i}$ is the
spin-degeneracy factor. The boson condensate term

$\tilde{\omega}_{d}^{(\mathrm{cond})}=\frac{1-\sigma_{i}}{2}n_{d}^{(\mathrm{cond})}(m^{\ast}_{d}-\mu_{d}^{\ast})$
(6)

in Eq. (II.1) with the density of the condensate, $n_{d}^{(\mathrm{cond})}$,
is only relevant for deuterons.

The effective chemical potential

$\mu_{i}^{\ast}=\mu_{i}-V_{i}.$ (7)

and the effective mass

$m_{i}^{\ast}=m_{i}+\Delta m_{i}-S_{i}$ (8)

depend on scalar

$S_{i}=\chi_{i\sigma}A_{i}C_{\sigma}n_{\sigma}$ (9)

and vector

$V_{i}=\chi_{i\omega}A_{i}C_{\omega}n_{\omega}+\chi_{i\rho}A_{i}C_{\rho}n_{\rho}+A_{i}V^{(r)}+W_{i}^{(r)}$
(10)

potentials, respectively. They are defined in terms of the different source
densities

$\displaystyle n_{\sigma}$ $\displaystyle=$
$\displaystyle\sum_{i}\chi_{i\sigma}A_{i}n_{i}^{(s)}$ (11) $\displaystyle
n_{\omega}$ $\displaystyle=$
$\displaystyle\sum_{i}\chi_{i\omega}A_{i}n_{i}^{(v)}$ (12) $\displaystyle
n_{\rho}$ $\displaystyle=$ $\displaystyle\sum_{i}\chi_{i\rho}A_{i}n_{i}^{(v)}$
(13)

where the single-particle scalar densities ($n_{i}^{(s)}$) and vector
densities ($n_{i}^{(v)}$) appear.

The mass shift $\Delta m_{i}$ in the effective mass (8) appears only for the
deuteron and is assumed to depend on the baryon density $n_{b}$. This
dependence leads to rearrangement contribution

$W_{i}^{(r)}=n_{d}^{(s)}\frac{\partial\Delta m_{d}}{\partial n_{i}^{(v)}}$
(14)

in the vector potentials (10) of nucleons and the deuteron.

The further rearrangement contribution

$V^{(r)}=\frac{1}{2}\left(C_{\omega}^{\prime}n_{\omega}^{2}+C_{\rho}^{\prime}n_{\rho}^{2}-C_{\sigma}^{\prime}n_{\sigma}^{2}\right)$
(15)

in the vector potential (10) is due to the density dependence of the couplings
$\Gamma_{j}$ and can be expressed with the coefficients $C_{j}$ and source
densities $n_{j}$ of the three mesons considered here.

Corresponding to the two rearrangement contribution in (10), there are also
two such terms in the total grand potential density (II.1): the meson term

$\tilde{\omega}_{\mathrm{meson}}^{(r)}=V^{(r)}n_{b}$ (16)

and the mass shift term

$\tilde{\omega}_{\mathrm{mass}}^{(r)}=\sum_{i}n_{i}^{(v)}W_{i}^{(r)}$ (17)

Furthermore, the meson contribution in (II.1) is given by

$\tilde{\omega}_{\mathrm{meson}}=-\frac{1}{2}\left(C_{\omega}n_{\omega}^{2}+C_{\rho}n_{\rho}^{2}-C_{\sigma}n_{\sigma}^{2}\right)$
(18)

similar in structure to (15).

The single-particle number densities can be derived from (II.1) using the
thermodynamic definitions

$\displaystyle n_{i}^{(v)}$ $\displaystyle=$
$\displaystyle-\left.\dfrac{\partial\tilde{\omega}}{\partial\mu_{i}}\right|_{T,\\{\mu_{j}\\}_{j\neq
i}}$ (19) $\displaystyle n_{i}^{(s)}$ $\displaystyle=$
$\displaystyle\left.\dfrac{\partial\tilde{\omega}}{\partial
m_{i}}\right|_{T,\\{\mu_{j}\\}}$ (20)

that give

$\displaystyle n_{i}^{(v)}$ $\displaystyle=$ $\displaystyle
g_{i}\int\dfrac{d^{3}k}{(2\pi)^{3}}d_{i}+\frac{1-\sigma_{i}}{2}n_{i}^{(\mathrm{cond})}$
(21) $\displaystyle n_{i}^{(s)}$ $\displaystyle=$ $\displaystyle
g_{i}\int\dfrac{d^{3}k}{(2\pi)^{3}}\dfrac{m_{i}^{\ast}}{E_{i}}d_{i}+\frac{1-\sigma_{i}}{2}n_{i}^{(\mathrm{cond})}$
(22)

with a thermal and a condensate contribution, once the distribution function

$d_{i}(T,k,m_{i}^{\ast},\mu_{i}^{\ast})=\left[\exp\left(\dfrac{E_{i}-\mu_{i}^{\ast}}{T}\right)+\sigma_{i}\right]^{-1}$
(23)

is defined. These expressions are consistent with the usual definitions of the
densities. The vector densities can be expressed as

$\displaystyle n_{n}^{(v)}$ $\displaystyle=$ $\displaystyle\frac{1+\beta-
X_{d}}{2}n_{b}$ (24) $\displaystyle n_{p}^{(v)}$ $\displaystyle=$
$\displaystyle\frac{1-\beta-X_{d}}{2}n_{b}$ (25) $\displaystyle n_{d}^{(v)}$
$\displaystyle=$ $\displaystyle\frac{X_{d}}{2}n_{b}$ (26)

using the baryon density $n_{b}$, the asymmetry $\beta$, and the deuteron
fraction $X_{d}$.

The deuteron fraction has to stay below a maximum value of

$X_{d}^{(\mathrm{max})}=\frac{m_{\mathrm{nuc}}}{\chi C_{\sigma}n_{b}}$ (27)

in order to ensure positive effective masses of the nucleons. This limit is
reached when the effective masses of the nucleons are zero, i.e.,
$S_{n}=S_{p}=C_{\sigma}n_{\sigma}=m_{\mathrm{nuc}}$ with the source density
$n_{\sigma}=2\chi n_{d}^{(s)}=\chi X_{d}^{\rm(max)}n_{b}$ of the $\sigma$
meson. It only has a contribution from the quasi-deuterons because the scalar
densities of the nucleons vanish for $m_{\mathrm{nuc}}^{\ast}=0$. A second
limitation of the deuteron fraction arises from the fact that for every
neutron a proton is needed, or vice versa, to form the cluster. This
translates to the condition $X_{d}\leq 1-|\beta|$ depending on the isospin
asymmetry $\beta$. So, in total, one has $0\leq
X_{d}\leq\mbox{min}\left\\{X_{d}^{(\mathrm{max})},1-|\beta|\right\\}$.

At high densities or temperatures, a mixture of deuterons, neutrons and
protons might be expected. If all three particle species have non-zero
densities, the condition of chemical equilibrium

$\displaystyle\mu_{d}=\mu_{n}+\mu_{p},$ (28)

applies between the chemical potentials of the degrees of freedom involved.
Using the definitions (7) and (8) with the potentials (9) and (10), this
relation can be written as

$\Delta m_{d}=\mu_{n}^{\ast}+\mu_{p}^{\ast}-m_{d}+S_{d}+V_{n}+V_{p}-V_{d},$
(29)

and thus an expression for the deuteron mass shift is obtained.

All thermodynamic quantities of the system can be easily obtained from the
grand canonical thermodynamic potential (II.1). For instance, the pressure is
given by

$P=-\omega(T,\\{\mu_{i}\\})$ (30)

and the free energy density $\mathcal{F}$ can be expressed as

$\mathcal{F}=\sum_{i=n,p,d}\mu_{i}n_{i}^{(v)}+\omega=\left(\mu_{b}+\frac{1-\beta}{2}\mu_{c}\right)n_{b}-P$
(31)

with the baryon chemical potential $\mu_{b}=\mu_{n}$ and the charge chemical
potential $\mu_{c}=\mu_{p}-\mu_{n}$. The entropy density finally can be
written as

$\displaystyle\mathcal{S}$ $\displaystyle=$
$\displaystyle-\left.\frac{\partial\omega}{\partial
T}\right|_{\\{\mu_{i}\\}}+\mathcal{S}_{\mathrm{cond}}$ $\displaystyle=$
$\displaystyle-\sum_{i}g_{i}\int\frac{d^{3}k}{(2\pi)^{3}}\>\left[d_{i}\ln
d_{i}+\frac{1-\sigma_{i}d_{i}}{\sigma_{i}}\ln\left(1-\sigma_{i}d_{i}\right)\right]$
$\displaystyle+n_{d}^{(\mathrm{cond})}\ln g_{d}$

with the distribution functions (23) and a contribution of the condensed
deuterons. The latter contribution arises because deuterons have spin $1$ and
thus the ground state of nuclear matter at $T=0$ contains a mixture of the
different spin substates.

###  II.2 Zero temperature and quasi-deuteron condensation

If the temperature vanishes, the vector and scalar densities of the the
nucleons $q=n,p$ can be expressed in analytical form as

$n_{q}^{(v)}=\frac{g_{q}}{6\pi^{2}}k_{q}^{3}$ (33)

with $g_{q}=2$ and

$n_{q}^{(s)}=\frac{g_{q}m_{q}^{\ast}}{4\pi^{2}}\left[k_{q}\mu_{q}^{\ast}-(m_{q}^{\ast})^{2}\ln\frac{k_{q}+\mu_{q}^{\ast}}{m_{q}^{\ast}}\right]$
(34)

with the Fermi momentum $k_{q}$ and the effective chemical potential

$\mu_{q}^{\ast}=\sqrt{k_{q}^{2}+\left(m_{q}^{\ast}\right)^{2}}\>.$ (35)

Since quasi-deuterons are bosons, they can exist only as a condensate in the
zero temperature case we are focusing on. For the deuteron there is no thermal
contribution to the density but the condensate term

$n_{d}^{(v)}=n_{d}^{(s)}=n_{d}^{(\mathrm{cond})}$ (36)

with equal vector and scalar densities and the effective chemical potential
becomes

$\mu_{d}^{\ast}=m_{d}^{\ast}\>.$ (37)

In nuclear matter at very low densities, it is advantageous to form quasi-
deuterons, which still have a positive binding energy, to gain energy as
compared to a system composed of neutrons and protons only. In SNM, all
protons and neutrons will be bound in quasi-deuterons in the low-density limit
and no free nucleons remain. The binding energy per nucleon will approach half
of the deuteron binding energy for $n_{b}\to 0$ and not zero as for
homogeneous nucleonic matter without clusters. For ANM, only quasi-deuterons
and neutrons (protons) will be the active constituents in case of positive
(negative) isospin asymmetry $\beta$.

With the help of the nucleon and deuteron densities, a simple expression for
the pressure $P$ is found. Since the the momentum integral (4) can be
calculated explicitly after partial integration, one obtains

$\displaystyle P$ $\displaystyle=$
$\displaystyle\sum_{q=n,p}\frac{1}{4}\left(\mu_{q}^{\ast}n_{q}^{(v)}-m_{q}^{\ast}n_{q}^{(s)}\right)$
$\displaystyle+\frac{1}{2}\left(D_{\omega}n_{\omega}^{2}+D_{\rho}n_{\rho}^{2}-D_{\sigma}n_{\sigma}^{2}\right)+\sum_{i=n,p,d}n_{i}^{(v)}W_{i}^{(r)}$

with

$D_{i}=C_{i}+C_{i}^{\prime}n_{b}$ (39)

for the meson couplings ($i=\sigma$, $\omega$, and $\rho$). Because the
entropy density (II.1) vanishes for $T=0$, the internal energy density
$\mathcal{E}$ is identical to the free energy density (31).

###  II.3 Coupling strength scaling factors

In Eqs. (9) - (13), the scaling factors $\chi_{ij}$ appear. These factors are
always unitary for nucleons and, in particular,

$\displaystyle\chi_{n\sigma}$ $\displaystyle=$
$\displaystyle\chi_{n\omega}=\chi_{n\rho}=1$ $\displaystyle\chi_{p\sigma}$
$\displaystyle=$ $\displaystyle\chi_{p\omega}=-\chi_{p\rho}=1$ (40)

for the three mesons. For the nucleons bound in clusters, the choice of the
scaling factors is instead widely debated Ferreira and Providencia (2012).

It is a natural choice to assume that the nucleons inside the deuteron couple
to the mesons with the same strength as the unbound nucleons. Nevertheless,
previous studies have already shown that, to take into account in-medium
effects in calculation of the EoS for astrophysical applications, a universal
scaling factor smaller than $1$ should be assumed for the cluster-meson
coupling strength Pais  _et al._ (2018, 2019). A reduced value for the
coupling of the $\sigma$-meson to different light clusters, including
deuterons, allows also a good description of the chemical equilibrium
constants determined from the NIMROD data Qin  _et al._ (2012). However,
recent Bayesian analysis Pais  _et al._ (2020a, b) have highlighted that a
larger value should be taken for the cluster $\sigma$-meson coupling, to
describe recent results from INDRA collaboration as well Bougault  _et al._
(2020). Moreover, a possible model-dependence of this result was recently
assessed Custódio  _et al._ (2020), calling for further microscopic analysis
to get more stringent constraints. On the other hand, the possible choice of
different values for the $\sigma$ and $\omega$ coupling strength factors is
known to produce a strong imbalance between the corresponding scalar and the
vector potentials. This would, in turn, reflect itself in an unrealistic
behavior of the nuclear EoS, owing to a corresponding change of the central
potential given by the difference $V_{i}-S_{i}$ of very large potentials in
lowest order non-relativistic approximation. A reasonable choice would be
therefore to assume the same scaling factor for both $\sigma$ and $\omega$
meson ($\chi_{d\sigma}=\chi_{d\omega}=\chi$) and explore the sensitivity of
our results to this ingredient ($\chi_{d\rho}=0$ because the deuteron has zero
isospin.) This is thus the strategy that will be adopted in the following.

##  III Quasi-deuteron mass shift constraints

In the GRDF, a mass shift $\Delta m_{d}$ is introduced in the effective mass
defined in Eq. (8) in order to suppress the cluster formation at supra-
saturation densities. This mass shift may generally receive several
contributions. For example, in compact star matter, where electrons are
included to fulfill the requirement of charge neutrality, a possible
contribution comes from the screening of the Coulomb potential produced by the
electronic background.

In the nuclear matter case, the mass shift is usually provided just from Pauli
blocking, which plays a crucial role in suppressing the cluster formation.
Indeed, the Pauli exclusion principle implies that a single-particle state in
momentum space would not be longer available for formation of a cluster when
it is already occupied by nucleons of the medium. The Pauli blocking of states
strongly reduces at high temperatures with increasing diffuseness of the Fermi
sphere or when the c.m. momentum of the cluster is much larger than the
typical radius of the Fermi sphere. This effect can be represented effectively
as a repulsive, medium-dependent potential or a change of the cluster binding
energy. A quantitative value can be calculated by solving the many-body
Schrödinger equation when proper potentials for the nucleons in matter are
introduced. The corresponding results, which are explicitly calculated for
various conditions of temperature, density and isospin asymmetry of the
medium, are then usually approximated by suitable parameterizations in a wide
range of thermodynamic variables Typel  _et al._ (2010); Röpke (2009, 2011,
2015).

The change of the cluster binding energy, calculated in this approach, is
valid only at sub-saturation densities and has to be extrapolated to higher
densities, in particular above the cluster dissociation (Mott) density, where
the binding energy vanishes, and only a many-body correlation in the continuum
remains. For this purpose, suitable heuristic parameterizations were
introduced within the GRDF Pais and Typel (2017b); Typel (2018). Focusing on
deuteron-like correlations, the aim of our work is to provide a unified
parameterization of the quasi-deuteron binding energy shifts, so that it can
be used as effective means to treat SRCs at supra-saturation densities. In
this context, the pioneering analysis developed in Ref. Typel (2020) already
suggests that a substantial modification of the cluster mass shift has to be
expected in the regime beyond the deuteron dissociation density, in comparison
to the traditionally adopted form.

Such a parameterization should be appropriately chosen to interpolate between
the low-density limit constrained by microscopic many-body calculations and
the high-density behavior postulated under the assumed boson condensation
condition. In the following subsections the different constraints to be
imposed to the deuteron binding-energy shifts will be discussed. As will be
seen, various scenarios will emerge depending also on the value of the
deuteron-meson scaling factors $\chi=\chi_{d\sigma}=\chi_{d\omega}$ introduced
in Eqs. (9) - (13).

###  III.1 Low-density constraint

Several parameterizations of the binding energy or mass shift of the deuteron
at sub-saturation densities have been developed to account for the Pauli
blocking effects. In particular, when neglecting the dependence of the mass
shift on the momentum of the deuteron with respect to the medium and limiting
to the zero temperature case, the simplified functional form

$\Delta m_{d}^{(\mathrm{low})}=\delta B_{d}(0)n_{d}^{(\mathrm{eff})}$ (41)

can be assumed, where the effective vector density $n_{d}^{(\mathrm{eff})}$,
as defined in Ref. Typel (2020)

$n_{d}^{(\mathrm{eff})}=\frac{2}{A_{d}}\left[N_{d}\left(n_{n}^{(v)}+n_{d}^{(v)}\right)+Z_{d}\left(n_{p}^{(v)}+n_{d}^{(v)}\right)\right]$
(42)

with $N_{d}=Z_{d}=1$ and $A_{d}=2$, turns out to be equal to the total baryon
density $n_{b}$, independent of the global isospin asymmetry $\beta$ of the
system. The quantity $\delta B_{d}$ generally regulates the temperature
dependence; its zero temperature limit which appears in Eq. (41) is $\delta
B_{d}(0)=3634.16$ MeV fm3, as given in Ref. Typel  _et al._ (2010), where
explicit expressions for finite temperature calculations are also provided.

A linear increase of the mass shift proportional to $n_{b}$ at low baryon
densities is consistent with the results obtained in Typel  _et al._ (2010);
Röpke (2009, 2011, 2015), at least for density values lying below the
dissociation or Mott density defined as

$n_{d}^{(\mathrm{diss})}=\dfrac{B_{d}}{\delta B_{d}(0)},$ (43)

where $B_{d}=m_{n}+m_{p}-m_{d}=2.225$ MeV Wang  _et al._ (2012) is the
deuteron binding energy in vacuum.

###  III.2 High-density limit

In the traditional treatment of cluster dissolution using the concept of mass
shifts, a heuristic density dependence stronger than linear in $n_{b}$ is
customarily assumed at baryon densities above the dissociation density to
prevent the clusters to reappear. A divergence of the mass shift ensures in
particular the deuteron removal from the system when approaching saturation
density, resulting in pure nucleonic matter above saturation. The possibility
of using quasi-deuterons to effectively embed nuclear SRCs at supra-saturation
densities requires thus a proper change of the usual parameterization adopted
in the GRDF, as first noticed in Ref. Typel (2020).

####  III.2.1 Deuteron mass shift

Assuming a dependence on the effective density $n_{d}^{(\mathrm{eff})}=n_{b}$,
as defined in Eq. (42), the quasi-deuteron mass shift expression derived in
Eq. (29) simplifies to

$\displaystyle\Delta m_{d}^{(\mathrm{high})}$ $\displaystyle=$
$\displaystyle\mu_{n}^{\ast}+\mu_{p}^{\ast}-m_{d}+2\chi
C_{\sigma}n_{\sigma}+2(1-\chi)C_{\omega}n_{\omega}$

because

$W_{n}^{(r)}+W_{p}^{(r)}-W_{d}^{(r)}=n_{d}^{(s)}\frac{\partial\Delta
m_{d}}{\partial n_{b}}\left(1+1-2\right)=0\>.$ (45)

Eq. (III.2.1) can be used to calculate the deuteron mass shift as a function
of the baryon density $n_{b}$, the deuteron fraction $X_{d}$ and the isospin
asymmetry $\beta$.

Figure 1:  Panel (a): Deuteron mass shift as function of the baryon density,
as determined according to Eq. (III.2.1) in the SNM case, by assuming a
unitary scaling factor $\chi$ for the deuteron-meson coupling strengths. Panel
(b): the same as in panel (a), but assuming a reduced scaling factor
$\chi=1/\sqrt{2}$. In both panels, the DD2 parameterization Typel  _et al._
(2010) of the nucleon-meson effective interaction is adopted. The inset in
panel (a) shows a zoom around the density values where the curves cross each
other. The thin blue line indicates in both panels the binding energy of the
deuteron in vacuum $B_{d}$.

An impression of the density dependence of the quasi-deuteron mass shift is
depicted in Fig. 1 for given, constant mass fractions $X_{d}$ and two values
of the scaling factor $\chi$ of the deuteron. A unitary value of $\chi$ is
considered in panel (a), whereas a reduced scaling factor of $\chi=1/\sqrt{2}$
is assumed in panel (b). The choice of the latter $\chi$ value will be
explained below. In both panels, for sake of simplicity, the SNM case is
considered. Moreover, the DD2 parameterization Typel  _et al._ (2010) of the
nucleon-meson effective interaction is adopted.

Fig. 1 shows that the deuteron mass shift $\Delta m_{d}$ depends quite
strongly on the value assumed for the scaling factor $\chi$. This is deduced
by comparing the results of panel (a) with the corresponding ones of panel
(b), especially in the high-density regime, we are focusing on. Indeed, for
large densities, the effective chemical potentials of the nucleons in Eq.
(III.2.1) are dominated by the Fermi momenta $k_{q}$ as defined in Eq. (33)
and thus exhibit a $n_{b}^{1/3}$ behavior. Since the meson couplings
$C_{\sigma}$ and $C_{\omega}$ approach constants at high densities, the
mesonic contributions to the mass shift are determined by the behavior of the
source densities $n_{\sigma}$ and $n_{\omega}$. For $\chi\neq 1$, the source
density $n_{\omega}$ rules the dominating term. Since the deuteron fraction
has to vanish for $n_{b}\to\infty$ due the constraint (27), a linear
asymptotic increase of $\Delta m_{d}$ with the baryon density is expected.
Indeed, such a dependence is observed in panel (b) of Fig. 1 for
$\chi=1/\sqrt{2}$. A softer increase of the mass shift with the baryon density
is seen in panel (a) because the last term in Eq. (III.2.1) does not
contribute for $\chi=1$ and the high-density behavior is driven by the
$\sigma$ meson term. However, the source density $n_{\sigma}$ may
asymptotically receive contributions only from the scalar densities of the
nucleons because the deuteron fraction asymptotically approaches zero. An
asymptotic dependence proportional to $n_{b}^{2/3}$ is then expected for the
$\sigma$ meson term and the deuteron mass shift of Eq. (III.2.1).

In the (unrealistic) case when all nucleons are bound inside the deuteron
($X_{d}=1$), independent of the density, both the scalar and vector densities
of the nucleons vanish. The mass shift defined in Eq. (III.2.1) assumes the
following form

$\Delta
m_{d}^{(\mathrm{high})}(X_{d}=1)=B_{d}(0)+2\chi(1-\chi)(C_{\omega}-C_{\sigma})n_{b}\>,$
(46)

i.e., a linear dependence on $n_{b}$ is generally predicted, for any finite
value of the scaling factor $\chi$, as shown by the dotted line of panel (b).
An exception is the case with $\chi=1$, when the deuteron mass-shift coincides
with the deuteron binding energy, see dotted line in panel (a).

Furthermore, Fig. 1 shows that the largest deuteron fraction generally
corresponds to the lowest mass shift, with a clear ordering in panel (b).
However, in the case when a unitary scaling factor is adopted, such a
statement actually holds only below the region where the lines cross each
other, which is observed above $0.45$ fm-3. This region is also emphasized in
the inset of panel (a) of Fig. 1 to evidence the fact that there is no single
crossing point between the different curves.

The emergence of this crossing can be easily understood when looking at the
explicit form of Eq. (III.2.1). Indeed, for $\chi=1$, the vector contribution
vanishes, and there exists a delicate interplay between the remaining terms.
At low densities, the leading role is played by the effective chemical
potentials, which reduce when the deuteron mass fraction increases. However,
at larger densities, the importance of the term involving the source density
of the scalar meson is enhanced. Since the corresponding contribution
increases with the deuteron fraction, a crossing among the curves is observed
at a certain baryon density $n_{b}^{\mathrm{cross}}$ and an inversion of their
previous ordering is expected at higher densities. On the other hand, for
smaller values of the scaling factor, as the one considered in panel (b) of
Fig. 1, the vector term plays also a role. The source density of the $\omega$
meson field does not change with the deuteron mass fraction at constant
$n_{b}$. Since the vector term dominates asymptotically over the scalar one,
in light of the power of density involved there, for scaling factor values
small enough, such as the one considered in panel (b) of Fig. 1, no crossing
is observed and the same ordering for all allowed densities is preserved.

Moreover, it is worthwhile to notice in panel (a) that the dotted-dashed line,
which corresponds to the deuteron mass shift obtained for a fixed deuteron
mass fraction $X_{d}=0.5$, ends at a density around $0.84$ fm-3. Beyond this
value, the deuteron mass fraction would exceed its maximum allowed value, see
Eq. (27). The same occurs also for the two largest deuteron mass fraction
values considered in both panels of Fig. 1, although at smaller densities. It
is also worth noticing that Eq. (27) implies that any finite asymptotic value
for the deuteron mass fraction is not allowed. As a result, the cluster is
forced to dissolve asymptotically, at least as $n_{b}^{-1}$, for any finite
value of $\chi$.

Figure 2:  Panel (a): Difference between the deuteron mass shift $\Delta
m_{d}$, as determined according to Eq. (III.2.1) and its value in the case
with zero deuteron mass fraction $\Delta m_{d}^{(0)}$, as function of the
deuteron mass fraction $X_{d}$ for different, but constant, values of the
baryon number density $n_{b}$. The SNM case is considered and a unitary
scaling factor $\chi$ for the deuteron-meson coupling strengths are assumed.
Panel (b): the same as in panel (a), but assuming a reduced scaling factor
$\chi=1/\sqrt{2}$. In both panels, the DD2 parameterization Typel  _et al._
(2010) of the nucleon-meson effective interaction is adopted.

The observation of the line crossing in panel (a) of Fig. 1 motivates to study
the dependence of the deuteron mass shift on the deuteron mass fraction for
constant baryon density $n_{b}$ and constant asymmetry $\beta$. A deeper
insight into this behavior can be achieved, by looking at Fig. 2, where
$\Delta m_{d}$ is plotted as a function of the deuteron mass fraction for
different baryon density values. In order to facilitate the comparison among
the curves, in Fig. 2 the differences with respect to the deuteron mass shift
evaluated in the limiting case without deuteron (denoted as $\Delta
m_{d}^{(0)}$) are actually considered. These differences turn out to be
systematically lower than zero, at least for the reduced value of the scaling
factor considered in panel (b). A change of sign is instead observed for
$\chi=1$ at larger densities, i.e., beyond the crossing point observed in the
panel (a) of Fig. 1 and already discussed above. Quite interestingly, one
observes in both panels that the (negative) slope of the curves strongly
increases for those lines that approach $X_{d}=1$. This represents the ideal,
but unrealistic, case where no free nucleons exist in the system also at
supra-saturation densities. The reason behind this behavior will be clarified
in the following, when the mass shift derivatives are investigated in detail.

####  III.2.2 Mass shift derivatives

In Appendix A, an explicit, general expression for the mass fraction
derivative of the mass shift is derived. For sake of simplicity, the case of
SNM is considered here, where the Fermi momenta, effective masses and
effective chemical potentials of the nucleons are identical, i.e.,
$k_{\mathrm{nuc}}=k_{n}=k_{p}$,
$m_{\mathrm{nuc}}^{\ast}=m_{n}^{\ast}=m_{p}^{\ast}$,
$\mu_{\mathrm{nuc}}^{\ast}=\mu_{n}^{\ast}=\mu_{p}^{\ast}$. Then the simple
form

$\displaystyle\left.\frac{\partial\Delta m_{d}^{(\mathrm{high})}}{\partial
X_{d}}\right|_{n_{b},\beta=0}=$
$\displaystyle\left[\frac{2C_{\sigma}}{1+f_{\mathrm{nuc}}C_{\sigma}}\left(\chi-\frac{m_{\mathrm{nuc}}^{\ast}}{\mu_{\mathrm{nuc}}^{\ast}}\right)^{2}-\frac{\pi^{2}}{k_{\mathrm{nuc}}}\frac{1}{\mu_{\mathrm{nuc}}^{\ast}}\right]n_{b}$

with the factor

$f_{\mathrm{nuc}}=3\left(\frac{n_{\mathrm{nuc}}^{(s)}}{m_{\mathrm{nuc}}^{\ast}}-\frac{n_{\mathrm{nuc}}^{(v)}}{\mu_{\mathrm{nuc}}^{\ast}}\right)$
(48)

and the total nucleon densities
$n_{\mathrm{nuc}}^{(v)}=n_{n}^{(v)}+n_{p}^{(v)}$ and
$n_{\mathrm{nuc}}^{(s)}=n_{n}^{(s)}+n_{p}^{(s)}$ is obtained. The derivative
(III.2.2) is the difference of two positive contributions and depending on the
choice of $\chi$ there can be a zero at a certain baryon density.

Figure 3:  Panel (a): Mass fraction derivative of the deuteron mass shift as
function of the baryon density, as determined according to Eq. (III.2.2) for
different, but constant values of the deuteron mass fraction $X_{d}$. The SNM
case is considered and a unitary scaling factor $\chi$ for the deuteron-meson
coupling strengths is assumed. Panel (b): The same as in panel (a), but
assuming a reduced scaling factor $\chi=1/\sqrt{2}$. In both panels, the DD2
parameterization Typel  _et al._ (2010) of the nucleon-meson effective
interaction is adopted. The thin blue line indicates the zero of the deuteron
mass fraction derivative of the mass shift. The inset in panel (a) shows a
zoom around the density where the derivative vanishes.

The dependence of the derivative (III.2.2) on the baryon density $n_{b}$ is
depicted in Fig. 3 panel (a) for $\chi=1$ and in panel (b) for
$\chi=1/\sqrt{2}$. The same constant values of $X_{d}$ as in Fig. 1 are
considered. Only the limit case with $X_{d}=1$ is not shown because for this
deuteron mass fraction, $k_{\rm nuc}$ vanishes and a (negative) divergent mass
shift derivative is obtained. Moreover, as already observed in Fig. 2, this
result is independent of the adopted value for the scaling factor.

Even though we are interested here in the investigation of the high-density
limit, it is instructive to study both the low- and high-density behaviors of
the derivative (III.2.2). One observes that, for $k_{\mathrm{nuc}}\to 0$, the
derivative is dominated by the negative, diverging contribution inside the
brackets and the derivative approaches thus zero from negative values, as
shown in both panels of Fig. 3. For $k_{\mathrm{nuc}}\to\infty$, the total
nucleonic vector and scalar densities scale as
$n_{\mathrm{nuc}}^{(v)}\sim[2/(3\pi^{2})]k_{\mathrm{nuc}}^{3}$ and
$n_{\mathrm{nuc}}^{(s)}\sim(1/\pi^{2})m_{\mathrm{nuc}}^{\ast}k_{\mathrm{nuc}}^{2}$,
respectively. Furthermore, $n_{b}\sim n_{\rm nuc}^{(v)}$,
$m_{\mathrm{nuc}}^{\ast}\sim\pi^{2}/(C_{\sigma}k_{\mathrm{nuc}}^{2})$, and
$\mu_{\mathrm{nuc}}^{\ast}\sim k_{\mathrm{nuc}}$, so that
$f_{\mathrm{nuc}}\sim k_{\mathrm{nuc}}^{2}/\pi^{2}$ and the simple asymptotic
form

$\left.\frac{\partial\Delta m_{d}^{(\mathrm{high})}}{\partial
X_{d}}\right|_{n_{b},\beta=0}\sim\left(2\chi^{2}-1\right)\frac{2}{3}k_{\mathrm{nuc}}$
(49)

remains. For $\chi<1/\sqrt{2}$ the derivative in the high-density limit is
negative as for $k_{\mathrm{nuc}}\to 0$ and no zero at finite baryon densities
is expected. In contrast, for $\chi=1$, the derivative approaches
asymptotically a positive value and a zero at a certain baryon value,
$n_{b}^{\rm cross}$, appears, as depicted in Fig. 3, panel (a). The curves
cross the zero line at similar values above $0.45$ fm-3 with a weak dependence
on the deuteron fraction, highlighted in the inset. Below the crossing an
ordering of the lines as observed in panel (b) of Fig. 1 follows, whereas an
inversion will appear above $n_{b}^{\mathrm{cross}}$. For $\chi=1/\sqrt{2}$,
however, the derivative of the mass shift with respect to the deuteron mass
fraction of Eq. (III.2.2) is always non-positive, as evidenced in the panel
(b) of Fig. 3. Then, the zero line is only asymptotically approached and
$\chi=1/\sqrt{2}$ constitutes the largest value for which the ordering of the
lines with respect to $X_{d}$, depicted in Fig. 1, panel (b) persists, for all
baryon densities.

Figure 4:  Panel (a): Baryon density derivative of the deuteron mass shift as
function of the baryon density, as determined according to Eq. (50) for
different, but constant values of the deuteron mass fraction $X_{d}$. The SNM
case is considered and a unitary scaling factor $\chi$ for the deuteron-meson
coupling strengths is assumed. Panel (b): The same as in panel (a), but
assuming a reduced scaling factor $\chi=1/\sqrt{2}$. In both panels, the DD2
parameterization Typel  _et al._ (2010) of the nucleon-meson effective
interaction is adopted.

In the following, also the derivative of the mass shift with respect to the
baryon density will be studied. It can be calculated explicitly from Eq.
(III.2.1). The general case for arbitrary values of $\beta$ is treated in
Appendix B. Again, the result for the simplified case of SNM is given here. It
can be written for $\chi=\chi_{d\omega}=\chi_{d\sigma}$ as

$\left.\frac{\partial\Delta m_{d}^{(\mathrm{high})}}{\partial
n_{b}}\right|_{\beta=0}=\mathcal{W}_{d}^{\mathrm{SNM}}-\mathcal{Z}_{d}^{\mathrm{SNM}}\left.\frac{\partial(n_{b}X_{d})}{\partial
n_{b}}\right|_{\beta=0}$ (50)

with the quantities

$\displaystyle\mathcal{Z}_{d}^{\mathrm{SNM}}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{\mu_{\mathrm{nuc}}^{\ast}k_{\mathrm{nuc}}}+2\left(1-\chi\right)^{2}C_{\omega}$
$\displaystyle-\frac{2C_{\sigma}}{1+f_{\mathrm{nuc}}C_{\sigma}}\left(\chi-\frac{m_{\mathrm{nuc}}^{\ast}}{\mu_{\mathrm{nuc}}^{\ast}}\right)^{2}$

and

$\displaystyle\mathcal{W}_{d}^{\mathrm{SNM}}=\frac{\pi^{2}}{\mu_{\mathrm{nuc}}^{\ast}k_{\mathrm{nuc}}}+2(1-\chi)\left(C_{\omega}+C_{\omega}^{\prime}n_{\omega}\right)$
$\displaystyle+\frac{2}{1+f_{\mathrm{nuc}}C_{\sigma}}\left(C_{\sigma}^{\prime}n_{\sigma}+C_{\sigma}\frac{m_{\mathrm{nuc}}^{\ast}}{\mu_{\mathrm{nuc}}^{\ast}}\right)\left(\chi-\frac{m_{\mathrm{nuc}}^{\ast}}{\mu_{\mathrm{nuc}}^{\ast}}\right)$

that contain again the factor (48). The dependence of the derivative (50) for
the SNM case is depicted in the two panels of Fig. 4 for the selected scaling
factors of $\chi=1$ (left) and $\chi=1/\sqrt{2}$ (right). In the zero-density
limit, a divergent behavior is observed, owing to the contribution originating
from the density derivatives of the effective chemical potentials that lead to
the terms proportional to $k_{\mathrm{nuc}}^{-1}$ in
$\mathcal{Z}_{d}^{\mathrm{SNM}}$ and $\mathcal{W}_{d}^{\mathrm{SNM}}$. The
divergence of (50) for $n_{b}\to 0$ will disappear only for constant
$X_{d}=1$. The relative importance of the terms proportional to
$k_{\mathrm{nuc}}^{-1}$ strongly reduces with increasing density, so that a
rise driven by the $\sigma$ meson contribution is observed at larger $n_{b}$
values, until a maximum is reached. Then, a continuous decrease of the density
derivative of the deuteron mass shift emerges, which asymptotically approaches
zero or a constant value for the unitary or the reduced value of $\chi$,
respectively. An interplay among the different involved terms takes place,
analogously to the one illustrated to describe the results of Fig. 1. As a
consequence, the ordering of the curves in Fig. 1 is reflected in Fig. 4. As
discussed before, for any finite deuteron fraction value $X_{d}$, the
corresponding curve ends again at the baryon density where
$X_{d}^{(\mathrm{max})}$, as a function of $n_{b}$, attains this value.

##  IV Mass shift parameterization

The condensation condition allows to calculate, for a given mass-fraction
function, the quasi-deuteron mass shift and its derivatives at high-densities
through Eqs. (III.2.1), (III.2.2) and (50), respectively. However, the density
dependence of the mass fraction $X_{d}$ is not known a priori, in particular
at supra-saturation densities. It should originate from microscopic
calculations in a similar manner as the fractions of light clusters are
determined by their mass shifts in the low-density domain, see Ref. Typel  _et
al._ (2010). Since calculations of the deuteron mass shift using proper
interactions and many-body methods are only available at very low densities,
it is necessary to resort to exemplary forms of the mass shift in the full
range of baryon densities to study the properties of the system. Instead of
calculating the mass shift and its derivatives for a given mass fraction
function, the main aim is then to choose a density dependent parameterization
of the mass shift and to determine the deuteron fraction, not only for
symmetric but also asymmetric matter. Although such a function is not yet
available from microscopic models, the proposed form should comply with the
available constraints.

###  IV.1 Piecewise mass-shift parameterization

Figure 5:  Panel (a): Energy per nucleon as function of the baryon density as
determined in the SNM case, by assuming a unitary scaling factor $\chi$ for
the deuteron-meson coupling strengths. A simple piecewise parameterization, as
given by Eq. (IV.1) and interpolating between the low- and high-density
constraints for the mass shift is adopted. Different, but constant, values of
the deuteron mass fraction $X_{d}^{\rm(high)}$ at high density are considered.
Panel (b): The same as in panel (a), but assuming a reduced scaling factor
$\chi=1/\sqrt{2}$. In both panels, the DD2 parameterization Typel  _et al._
(2010) of the nucleon-meson effective interaction is adopted.

The most simple choice for the deuteron mass shift function is the piecewise
parameterization

$\displaystyle\Delta m_{d}(n_{b})=$ $\displaystyle\mbox{min}\left\\{\Delta
m_{d}^{({\rm low})}(n_{b}),\Delta m_{d}^{({\rm high})}(n_{b},X_{d})\right\\}$

that combines the low-density form (41) with the high-density function
(III.2.1) assuming, e.g., a constant deuteron fraction in the high-density
region as discussed in the previous subsection. Although Eq. (IV.1) provides a
continuous function, the same does not apply for its derivatives. In addition,
a constant deuteron fraction at high densities is not compatible with the
asymptotic constraint which imposes a vanishing mass fraction $X_{d}$ for
$n_{b}\to\infty$.

Despite these shortcomings, it is informative to investigate the effect of
finite deuteron fractions on the energy per nucleon $E/A$. Fig. 5 depicts the
density behavior of $E/A$ in SNM as obtained with the piecewise mass shift
parameterization (IV.1) for different, but constant mass fractions in the
high-density regime. For comparison, the energy per nucleon obtained from the
deuteron free case of the standard DD2 parameterization is also shown. One
observes that an extra binding is predicted in the energy per nucleon around
saturation density, with respect to the case without deuterons. This result is
independent of the adopted choice of the scaling factor $\chi$ of the
deuteron-meson couplings. Moreover, Fig. 5 shows that the overbinding persists
much beyond the saturation density, especially in panel (b), where a reduced
value of $\chi$ is considered. However, the stiffness of these curves in the
high-density regime increases with the deuteron mass fraction. As a result, a
crossing among the curves is generally observed at larger densities in case
when $\chi\geq 1/\sqrt{2}$, except for deuteron fraction values for which the
corresponding curves end at lower densities.

Both insets of Fig. 5 highlight that the correct low-density limit is
reproduced because the low-density constraint is properly taken into account.
As a result, differently from the standard DD2 model without clusters, the
curves tend to half the deuteron binding energy in the limit $n_{b}\to 0$.

The overbinding observed in Fig. 5 in a broad range of densities around
saturation seems to be surprising in light of the large positive mass shift of
the deuteron depicted in Fig. 1. The increase of the binding energy per
nucleon is a result of two main effects. When nucleons are replaced by quasi-
deuterons at a given baryon density, the value of $n_{\omega}$ does not change
when $\chi=1$. At the same time, the source density $n_{\sigma}$ increases
because the scalar and vector densities of the condensed deuterons are
identical whereas $n_{q}^{(s)}/n_{q}^{(v)}<1$ for nucleons. Hence, on the one
hand, a stronger attraction from the $\sigma$ meson is induced. This effect
will be smaller for $\chi<1$ when nucleons are replaced by deuterons. Here,
the source density of the $\omega$ meson will reduce but the corresponding
decrease of the $\sigma$ meson source density is also less strong. On the
other hand, nucleons at the Fermi surface with energies close to the nucleon
chemical potential

$\mu_{q}^{(0)}=\sqrt{\left(k_{q}^{(0)}\right)^{2}+\left(m_{q}-S_{q}^{(0)}\right)^{2}}+V_{q}^{(0)}$
(54)

of the deuteron-free system, indicated by a superscript $(0)$, are replaced by
deuterons with energy

$E_{d}=\mu_{d}=m_{d}-S_{d}+V_{d}+\Delta m_{d}=\mu_{n}+\mu_{p}\>.$ (55)

Since $k_{q}<k_{q}^{(0)}$, $S_{q}>S_{q}^{(0)}$, and $V_{q}\approx
V_{q}^{(0)}$, the chemical potentials of the nucleons $\mu_{q}$ are lowered,
i.e., $\mu_{q}<\mu_{q}^{(0)}$, even for a positive deuteron mass shift $\Delta
m_{d}$. This corresponds to a stronger binding of the system.

###  IV.2 Specific features in the determination of deuteron mass fractions
from mass shifts

In principle, it would be sufficient to use the condensation condition
(III.2.1) to find $X_{d}$ for a given $\Delta m_{d}^{(\mathrm{high})}$. In
practice, however, it is found that, in the high-density regime, Eq. (III.2.1)
can have multiple solutions when a meson coupling scaling factor $\chi\geq
1/\sqrt{2}$ is considered. As a consequence, Eq. (III.2.1) can not be inverted
uniquely in all cases. The correct solution has to be selected such that a
continuous function of the density is obtained for the deuteron mass fraction.

In this context, further insight can be obtained from the explicit functional
form of the mass shift derivative with respect to the density, given by Eq.
(50). It is, in fact, a first-order differential equation for $X_{d}$. The
expression can be analyzed most easily for SNM. It is convenient to write Eq.
(50) in the form

$\displaystyle\left.\frac{\partial(n_{b}X_{d})}{\partial
n_{b}}\right|_{\beta=0}=$
$\displaystyle\frac{1}{\mathcal{Z}_{d}^{\mathrm{SNM}}}\left(\mathcal{W}_{d}^{\mathrm{SNM}}-\left.\frac{\partial\Delta
m_{d}}{\partial n_{b}}\right|_{\beta=0}\right)$

with the functions $\mathcal{Z}_{d}^{\mathrm{SNM}}$ and
$\mathcal{W}_{d}^{\mathrm{SNM}}$ as defined in Eqs. (III.2.2) and (III.2.2),
respectively. A special role is played by $\mathcal{Z}_{d}^{\mathrm{SNM}}$
that can be written as

$\mathcal{Z}_{d}^{\mathrm{SNM}}=2\left(1-\chi\right)^{2}C_{\omega}-\frac{1}{n_{b}}\left.\frac{\partial\Delta
m_{d}}{\partial X_{d}}\right|_{n_{b},\beta=0}$ (57)

with the help of the mass fraction derivative defined in Eq. (III.2.2). The
superscript (high) is no longer added to the mass shift because a mass shift
parameterization is considered now in the whole density range. For $\chi=1$,
the term from the $\omega$ meson does not contribute. The remaining term in
Eq. (57) develops a zero at a certain density $n_{b}^{\mathrm{cross}}$ as
discussed in section III.2 close to Eq. (49). Thus the density derivative
(IV.2) develops a pole at $x_{b}^{\mathrm{cross}}$ and a continuous solution
of the differential equation can only be obtained if the term in parentheses
in Eq. (IV.2) vanishes at the same density. For $\chi\leq 1/\sqrt{2}$,
however, the mass-fraction derivative (III.2.2) is negative, implying a
positive $\mathcal{Z}_{d}^{\mathrm{SNM}}$ because also the contribution of the
$\omega$ meson is positive. Then, a continuous solution of the differential
equation (IV.2) can be found for all densities. The discussion developed above
justifies the choice of the reduced scaling factor. Thus, two different values
of the scaling factor, namely $\chi=1$ and $\chi=1/\sqrt{2}$, will be
considered in the following analysis. It is worthwhile to notice that the
value of $\chi=1/\sqrt{2}$ is significantly smaller than the universal scaling
factor for the cluster-meson coupling strength proposed in some recent works
Pais  _et al._ (2018, 2019). Our choice complies, however, with the aim to
consider two extreme values of $\chi$ with two distinct paths, bearing in mind
that any intermediate behavior may also occur.

###  IV.3 Saturation constraints

The overbinding observed around saturation density in both panels of Fig. 5
implies that a proper refit of the nucleon-meson couplings to the saturation
properties of SNM is mandatory if one wants to keep nuclear matter quantities
around the saturation point $n_{0}$ well constrained. The actual deuteron
fraction at saturation, $X_{d,0}$, which has to be specified to fix the
couplings, can be imposed from recent experimental investigations of SRCs by
extrapolating results of nuclei to infinite nuclear matter. They assess that
SRCs pairs amount to approximately $20\%$ of the nucleon density Egiyan  _et
al._ (2006); Subedi  _et al._ (2008); Hen  _et al._ (2014). Here, as in the
following, the index $0$ on the quantities indicates the values at saturation.

Table 1: Values at saturation density of nucleon-meson coupling strengths $\Gamma_{j,0}$ ($j=\sigma,\omega,\rho)$), deuteron mass shift $\Delta m_{d,0}$ and its slope, for the two different scaling factors $\chi$ considered in this work. The mass shift and its slope are expressed in MeV and MeV fm3, respectively, whereas all other quantities are dimensionless.  $\chi$ | $\Gamma_{\sigma,0}$ | $\Gamma_{\omega,0}$ | $\Gamma_{\rho,0}$ |  $\Delta m_{d,0}$ [MeV] |  $\left.\frac{d\Delta m_{d}}{dn}\right|_{n_{0},\beta=0}$ [MeV fm3]
---|---|---|---|---|---
1 | $10.580042$ | $13.217226$ | $3.556424$ | $104.92$ | $813.98$
$1/\sqrt{2}$ | $10.919963$ | $13.719324$ | $3.400187$ | $58.23$ | $570.80$

In order to reproduce the properties at saturation in SNM, e.g., of the DD2
model, the binding energy per nucleon $B_{0}$ and the effective nucleon mass
$m_{\mathrm{nuc},0}^{\ast}$ at the saturation density should be obtained also
in the model with a finite deuteron fraction $X_{d,0}$. This will be realized
by rescaling the meson-nucleon couplings assuming no change in their density
dependence as given in a reference parameterization. The values of $n_{0}$,
$B_{0}$, $m_{\mathrm{nuc},0}$ and $X_{d,0}$ together with the pressure
$P_{0}=0$ MeV fm-3 give four conditions that allow to determine the rescaled
couplings $\Gamma_{\sigma,0}$, $\Gamma_{\omega,0}$, the deuteron mass shift
$\Delta m_{d,0}$ and its derivative $\left.d\Delta
m_{d}/dn_{b}\right|_{n_{0}}$ at saturation. Owing to the rescaling of the
$\sigma$ and $\omega$ coupling strengths, also the energy per nucleon of PNM
would be modified. A proper rescaling also of the $\rho$-nucleon coupling
strength is hence in order if one wants to keep the symmetry energy at
saturation $J_{0}$ unaltered, with respect to the selected reference
parameterization. Here, the symmetry energy is calculated in the parabolic
approximation as the difference between the energies per nucleon in PNM and
SNM. This recipe gives finite values at sub-saturation densities in models
with clusters or liquid-gas phase transition differently than the original
definition using second derivatives of the energy per nucleon in SNM with
respect to the isospin asymmetry, see, e.g., Ref. Typel  _et al._ (2014). The
full conversion procedure is illustrated in detail in Appendix C. The actual
values of these quantities are given in Table 1 for the two considered values
of the deuteron coupling scaling factor $\chi$ and a deuteron fraction
$X_{d,0}=0.2$. The standard DD2 model is chosen as reference parameterization
with $n_{0}=0.149065$ fm-3, $B_{0}=16.0224$ MeV,
$m_{\mathrm{nuc},0}^{\ast}=0.562544\leavevmode\nobreak\ m_{\mathrm{nuc}}$ and
$J_{0}=32.73$ MeV. The $\sigma$, $\omega$, and $\rho$ nucleon-meson couplings
are given by $\Gamma_{\sigma,0}=10.686681$, $\Gamma_{\omega,0}=13.342362$ and
$\Gamma_{\rho,0}=3.626940$ in this case.

Quite interestingly, one observes that, in the case with a unitary scaling
factor $\chi$, a small reduction of the three nucleon-meson coupling strengths
is obtained at saturation, with respect to the original DD2 parameterization.
On the other hand, these quantities turn out to be larger for both $\sigma$\-
and $\omega$-meson, when the smaller $\chi$ value is considered. This change
reflects the balance between the couplings and scaling factors to achieve the
same strength of the effective interaction at saturation density. Concerning
the mass shifts, larger values and stronger density slopes are predicted at
saturation when assuming that the nucleons bound inside the deuterons couple
to the mesons with the same strength as the unbound nucleons, see also Fig. 1.

A possible mass shift parameterization will be proposed in the following
section. It will be constrained at saturation density, in the low-density
limit by microscopic many-body calculations and at high-density by an assumed
mass-fraction behavior that respects the condensation condition and the
constraint on the maximum deuteron fraction.

###  IV.4 Unified mass shift parameterization

Different forms of the mass shift density parameterization might be employed
to interpolate among the low- and the high-density constraints discussed in
the previous section, while keeping the required saturation properties. A
unified form, adopted for the SNM case, might be then employed also for ANM to
obtain predictions at arbitrary isospin asymmetries. The mass shift
parameterization introduced in this work combines two limiting dependencies on
the density that reproduce by construction the linear increase in the zero-
density limit of the dilute region and the high-density asymptotic behavior.
Taking into account also the two constraints imposed at saturation, such a
parameterization has to depend on at least four parameters and satisfy $\Delta
m_{d}(0)=0$. However, some additional parameters should enter in the proposed
parameterization to guarantee a smooth transition between the different
density regimes and leave some freedom in the high-density behavior of the
deuteron fraction.

A possible choice, among others, is the function

$\Delta m_{d}(x)=\Delta m_{d,1}(x)+\Delta m_{d,2}(x)+\Delta m_{d,3}(x)$ (58)

depending on $x=n_{b}/n_{0}$ with three contributions

$\displaystyle\Delta m_{d,1}(x)$ $\displaystyle=$
$\displaystyle\dfrac{ax}{1+bx}$ (59) $\displaystyle\Delta m_{d,2}(x)$
$\displaystyle=$ $\displaystyle cx^{\eta+1}\left[1-\tanh{(ex)}\right]$ (60)
$\displaystyle\Delta m_{d,3}(x)$ $\displaystyle=$ $\displaystyle
fx^{\gamma}\tanh{(gx)}$ (61)

with $\gamma=1$ or $2/3$ for $\chi=1/\sqrt{2}$ or $1$, respectively, and seven
coefficients $a$, $b$, $c$, $\eta$, $e$, $f$, $g$, which allow to comply with
the constraints by a proper choice.

The different terms in Eqs. (58) - (61) are chosen so that the density
derivative of the mass shift is given by

$\left.\frac{d\Delta m_{d}}{dn_{b}}\right|_{n_{b}=0}=\frac{a}{n_{0}}$ (62)

at vanishing $n_{b}$ and the high-density behavior is dominated by the third
term $\Delta m_{d,3}$. The second contribution acts mainly in an intermediate
density range. The coefficient $a$ is determined as $\delta B_{d}(0)n_{0}$ by
the limiting form of the deuteron mass shift parameterization from microscopic
calculations, c.f., Eq. (41). There is no a priori constraint for the
parameter $b$. Here, it is set to $b=a/B_{d}=n_{0}/n_{d}^{(\mathrm{diss})}$ so
that $\lim_{n_{b}\to\infty}\Delta m_{d,1}=B_{d}$ with the deuteron binding
energy $B_{d}$ and the dissociation density $n_{d}^{(\mathrm{diss})}$ defined
in Eq. (43). In the asymptotic limit, for $\chi=1/\sqrt{2}$, the mass shift
approaches a linear function in the baryon density and a slope determined by
the ratio of $f$ and the saturation density $n_{0}$. On the other hand, the
asymptotic form

$\Delta m_{d}(n_{b})\sim f\left(\frac{n_{b}}{n_{0}}\right)^{2/3}$ (63)

is expected for the deuteron mass shift, in the case when $\chi=1$.

The coefficients $c$ and $\eta$, whose analytical expressions are given in
Appendix D, are finally determined by the constraints introduced on the mass
shift and its derivative at saturation. The remaining parameters $e$ and $g$
are free and allow to tune the relative role of the different contributions
$\Delta m_{d,2}$ and $\Delta m_{d,3}$ in Eq. (58). Only a tiny sensitivity of
the results was assessed by varying the parameter $e$. Thus, this parameter
was kept fixed to $1$ and only different values for the parameter $g$ were
considered. Different choices of $g$ permit indeed to produce alternative
supra-saturation scenarios while keeping the same asymptotic behavior.

####  IV.4.1 Deuteron mass shift and mass fraction for $\chi=1$

Table 2: Values of the parameters in the deuteron mass shift parameterization defined by Eqs. (58) - (61) for six different sets. They are obtained by employing the DD2 nucleon-meson effective interaction, with properly rescaled meson coupling strengths at saturation. The parameters $a$, $c$ and $f$ are expressed in MeV, $b$, $\eta$, $e$, $g$ and $\gamma$ are dimensionless. The first three sets refer to the case with deuteron-meson coupling scaling factor $\chi=1$, the others to $\chi=1/\sqrt{2}$.  | $a$ | $b$ | $c$ | $\eta$ | $e$ | $f$ | $g$ | $\gamma$
---|---|---|---|---|---|---|---|---
DD2 - d1 | $541.726060$ | $243.472387$ | $-83.230901$ | $3.491787$ | $1.0$ | $214.368137$ | $0.65$ | $2/3$
DD2 - d2 | $541.726060$ | $243.472387$ | $-98.923123$ | $3.200967$ | $1.0$ | $214.368137$ | $0.67632$ | $2/3$
DD2 - d3 | $541.726060$ | $243.472387$ | $-140.309501$ | $2.715545$ | $1.0$ | $214.368137$ | $0.75$ | $2/3$
DD2 - $\chi$d1 | $541.726060$ | $243.472387$ | $99.677247$ | $1.656159$ | $1.0$ | $181.113975$ | $0.18$ | $1$
DD2 - $\chi$d2 | $541.726060$ | $243.472387$ | $70.476986$ | $1.230947$ | $1.0$ | $181.113975$ | $0.22$ | $1$
DD2 - $\chi$d3 | $541.726060$ | $243.472387$ | $41.777908$ | $0.257252$ | $1.0$ | $181.113975$ | $0.26$ | $1$

Figure 6:  Panel (a): Deuteron mass shift as function of the baryon density,
as determined according to the parameterization proposed in Eqs. (58) - (61).
The SNM case is considered and a scaling factor $\chi=1$ is assumed for the
deuteron-meson coupling strengths. Three different set of parameters were
employed. The red shaded areas evidence the region of allowed mass shift or
deuteron mass fraction values. Panel (b): Deuteron mass fraction $X_{d}$ as
function of the baryon density, as determined by employing the same sets of
parameters considered in panel (a). In both panels, the DD2 nucleon-meson
effective interaction, with properly rescaled meson coupling strengths at
saturation, is adopted. Two curves with symbols (not shown in panel (a)) are
considered in panel (b), as a result of slightly varying the $g$ parameter of
the set labeled as DD2-d2 (see text for more details).

Let us consider first the case with $\chi=1$. For such a scaling factor, three
different sets of parameters will be considered. They are determined according
to the parameterization proposed in Eqs. (58) - (61) with $\gamma=2/3$ and are
labeled in the following as DD2-d1, DD2-d2 and DD2-d3. The values of the
parameters for these mass shift parameterizations are listed in the first
three lines of Table 2. They were obtained by employing the DD2 nucleon-meson
effective interaction with properly rescaled meson coupling strengths at
saturation as given in Table 1. Panel (a) of Fig. 6 shows the deuteron mass
shift as function of the baryon density for these three different sets. The
red shaded area corresponds to the range of possible mass shift values that
are explored by assuming, for each density, deuteron fractions within the
range $[0,\mbox{min}\\{1,X_{d}^{\rm(max)}\\}]$. First of all, Fig. 6
highlights the validity of the proposed parameterization to comply with the
constraints imposed on the deuteron mass-shift. The three lines lie indeed
within the red area for the whole range of displayed baryon densities.

However, as shown in Table 2, a rather small range of the parameter $g$ may be
actually explored, owing to the shrinkage of the red area in the density
region around the crossing points, which were observed in Fig. 1, panel (a).
Despite their proximity, when these parameterizations are employed in Eq.
(IV.2) to determine the density behavior of the deuteron fraction, extremely
different outcomes are obtained, at least when approaching the region of the
crossings. As clearly depicted in Fig. 6, panel (b), the three adopted
parameterizations correctly reproduce the low-density limit, corresponding to
the situation in which the matter is entirely clusterized ($X_{d}=1$), and
account also for the constraints at saturation.

Nonetheless, highly diverse scenarios are predicted at high-densities. The
dashed curve, which is the lowest one in the region around the shrinkage
observed in panel (a), corresponds to a deuteron mass fraction which tends to
exceed, at a certain density, the maximum allowed value for
$X_{d}^{\rm(max)}$. We remind that the density region beyond this point would
be characterized by a negative value of the Dirac effective mass of the
nucleons. To exclude this opportunity, the curve in Fig. 6, panel (b) is thus
stopped beyond that density. Owing to the unfeasible high-density scenario,
this set of parameter will not be further employed in the following.
Conversely, the DD2-d3 parameterization, plotted as the dotted line, lying
above the other curves around $0.5$ fm-3 in panel (a), predicts a sudden
disappearance of the deuterons, when approaching the pole of Eq. (IV.2).
Actually, even in this case, the curve is not plotted beyond the region of the
crossing, to exclude the unrealistic situation in which the clusters reappear
at very high densities. As already anticipated in Section IV.2, a quasi-
continuous solution of Eq. (IV.2) might be arranged for all densities with a
fine tuning of the parameters. The corresponding curve is plotted in both
panels of Fig. 6 as the full line. In this case, a smooth behavior is
apparently recovered for the density behavior of the deuteron mass fraction.
However, a strong sensitivity persists in correspondence of the pole, as
manifested by the two thin lines with symbols plotted in panel (b). These
curves are obtained with mass shifts functions that are found by varying the
$g$ parameter only by 1% of the set labeled as DD2-d2. Owing to the presence
of the pole, the DD2-d2 set of parameters will be set aside hereafter too. A
more refined method to find a continuous function would be to consider a more
general form of the mass shift parameterization with respect to the one
proposed in Eqs. (58)-(61). By enlarging the number of the involved
parameters, it would be possible to constrain the mass-shift and its slope
values such that the pole will be definitely washed out. Nonetheless, none of
such possible parameterizations would allow to accomplish our aim to extend
the predictions to ANM. Since the position of the pole evolves with the
asymmetry, a divergence of the mass fraction would emerge once again as soon
as the asymmetry of the matter is changed, despite its removal in the SNM
case. The same holds, obviously, for any deuteron-meson coupling scaling
factor $\chi>1/\sqrt{2}$. However, one should bear in mind that, for values of
the scaling factor smaller than $1$, the pole is expected to appear at higher
densities. Since for $\chi=1$ the pole emerges already at rather large
densities (around $3n_{0}$), one expects that for more realistic values of
$\chi$, its position will be located much beyond the range that is relevant in
the applications of the model. Then its emergence could be neglected in
practice. For $\chi=1$, only the DD2-d3 parameterization will be employed in
the following, when the general properties of both SNM and ANM matter will be
investigated.

####  IV.4.2 Deuteron mass shift and mass fraction for $\chi=1/\sqrt{2}$

As discussed in the previous sections, a possible way out to overcome the
issue of the pole might be to assume a smaller scaling factor which is, at
maximum, equal to $1/\sqrt{2}$. For such a scaling factor, three different
sets of parameters, labeled as DD2-$\chi$d1, DD2-$\chi$d2, and DD2-$\chi$d3,
are proposed here. The values of the parameters in Eqs. (58) - (61) for these
sets are listed again in Table 2.

One observes a strong sensitivity to the $g$ parameter of the deduced values
for $c$ and $\eta$. Moreover, when increasing $g$, both parameters to account
for the saturation constraints decrease. A further increase of $g$ beyond a
maximum value $g_{\rm max}\approx 0.26$ is excluded, since it would imply a
negative value for $d$ and thus a dominant role of $\Delta m_{d,2}$ in the
zero-density limit, where a pole could even emerge. However $\Delta m_{d,1}$
returns already, by construction, the correct low-density trend.

Figure 7:  Panel (a): Deuteron mass shift as function of the baryon density,
as determined according to the parameterization proposed in Eqs. (58) - (61).
The SNM case is considered and a scaling factor $\chi=1/\sqrt{2}$ is assumed
for the deuteron-meson coupling strengths. Three different set of parameters
were employed. The red shaded area evidences the region between the (upper)
curve, related to the deuteron-free case, i.e. $X_{d}=0$, and the (lower) one,
obtained by assuming
$X_{d}=\mbox{min}\left\\{1,X_{d}^{(\mathrm{max})}\right\\}$. Panel (b):
Deuteron mass fraction $X_{d}$ as function of the baryon density as determined
by employing the same sets of parameters of panel (a). In both panels, the DD2
nucleon-meson effective interaction, with properly rescaled meson coupling
strengths at saturation, is adopted.

Further insights may be achieved by looking at both panels of Fig. 7. The
density behaviors of these mass shifts are shown in Fig. 7 panel (a). The same
parameterizations are then employed in panel (b) of Fig. 7 to determine the
corresponding density behavior of the deuteron mass fraction $X_{d}$. In the
two panels, the red shaded area evidences the allowed region between the
(upper) curve, related to the deuteron-free case, i.e. $X_{d}=0$, and the
(lower) one, obtained by assuming
$X_{d}=\mbox{min}\left\\{1,X_{d}^{(\mathrm{max})}\right\\}$.

First of all, looking at panel (a), one notices that the black curves lie
always within the red shaded area up to very large baryon densities, so
validating the choice of the adopted parameter sets. In light of the
constraints imposed in the extremely dilute regime and at saturation, the
full, dashed and dottes black curves remain rather close up to $n_{0}$. All
the curves also converge to the line characterized by $X_{d}=0$, in the
asymptotic limit. Some differences emerge instead in the high-density
behavior, around and beyond $3n_{0}$.

The observed differences in the mass shifts are then reflected in the density
behavior of the deuteron mass fraction, which is plotted in panel (b).
Independent on the parameterization, in the zero-density limit, the matter is
completely clusterized and the deuteron mass fraction $X_{d}$ is equal to 1.
With increasing density, a considerable reduction of $X_{d}$ is observed and a
local minimum emerges in a density region around the saturation density. There
the value $X_{d}=0.2$ is reached, as required in agreement with the
experimental evidences concerning the emergence of SRCs pairs. At supra-
saturation densities, several scenarios take place. The curves never overshoot
the line indicating the maximum allowed value for $X_{d}$, which ensures a
non-negative value for the Dirac effective mass of the nucleons. Then, at
higher densities, a decreasing trend is observed for all curves, which
converge each other, approaching zero asymptotically.

It is worthwhile to notice that alternative scenarios, similar to the ones
displayed by the dotted and the dashed lines in panel (b) of Fig. 6, would be
possible also in the case with $\chi=1/\sqrt{2}$. Solutions with deuteron mass
fraction values which tend to exceed the maximum allowed at a certain density
or abruptly vanishing might accidentally occur, when considering mass shift
parameterizations which cross the lowest or the highest border, respectively,
of the red shaded area shown in panel (a) of Fig. 7. Differently than the case
with $\chi=1$, these solutions are not connected to the emergence of any pole.
The three sets of parameters proposed in panel (b) of Fig. 7 avoid these
scenarios. Solutions with the disappearance of the clusters at a certain
baryon density would be likewise acceptable, even though none of the three
chosen sets of parameters for $\chi=1/\sqrt{2}$ provides a similar result.
Indeed, for this class of solutions, the density behavior of the deuteron mass
fraction would closely resemble the one obtained with the DD2-d3
parameterization. The main difference will be only the possible wider density
range with a non-vanishing deuteron mass fraction values before the cluster is
suppressed. For $\chi=1/\sqrt{2}$, we will consider three set of parameters,
such that the mass shifts are characterized by similar smooth trends as
functions of the baryon density, but different sizes of the deuteron mass
fractions in the supra-saturation density regime. In such a way, we could
assess and isolate the role of this ingredient.

The sets of parameters DD2-d3, DD2-$\chi$d1, DD2-$\chi$d2, and DD2-$\chi$d3
listed in Table 2 return a smooth behavior of the mass fraction for the whole
range of baryon densities and for any asymmetry. Thus, they will be employed
in the next section to study various properties, both for SNM and ANM.

##  V Thermodynamic quantities

Once the density behavior of the mass shift and of the corresponding deuteron
mass fraction is determined, it is interesting to see how the embedding of
quasi-clusters at supra-saturation densities affects some general
thermodynamic quantities.

###  V.1 SNM: EoS and incompressibility

Figure 8:  Energy per nucleon as a function of the baryon density, as
determined by employing four selected set of parameters listed in Table 2, for
the mass shift parameterization given in Eqs. (58)-(61). For comparison, the
curve obtained for the DD2 parameterization in the deuteron-free case is also
shown, as the red full line. The inset shows the differences with respect to
the predictions of the standard DD2 parameterization in a larger baryon
density range.

Let us focus in this section on the results for SNM. The density dependence of
the energy per nucleon $E/A$ is plotted in Fig. 8. The sets of parameters
DD2-d3, DD2-$\chi$d1, DD2-$\chi$d2, DD2-$\chi$d3 listed in Table 2 are
employed. The standard DD2 parameterization, describing the deuteron free
case, is also shown as reference for comparison. The inset of Fig. 8 depicts
in particular the differences between the energy per nucleon derived with each
set and the chosen reference. First of all, one notices that, in light of the
fit performed at saturation, all the curves remain rather close in the low-
density regime. However, the parameterizations with deuterons differ from the
standard DD2 result in the zero-density limit, approaching one half of the
deuteron binding energy in vacuum. Moreover, remarkable differences also
emerge in the high-density behavior of the energy per nucleon. Despite the fit
performed at saturation, as in Fig. 5, a stronger binding is generally
observed in the neighbourhood of $n_{0}$ in the parameterizations accounting
for the presence of deuterons. This stronger attraction persists also to all
densities in the case of the DD2-d3 set, despite the disappearance of the
deuterons which is expected to occur at $n_{b}$ around 0.45 fm-3 with this
parameterization (see Fig. 6). The reduction in the energy per nucleon
observed at higher densities with the DD2-d3 is thus only driven by the
changed balance between the scalar and vector components, which is a result of
the rescaling of the meson coupling strengths at saturation.

A different scenario manifests itself with the parameterizations derived with
a reduced scaling factor $\chi$. There, no systematic increased binding, as
compared to the DD2 case, is observed. A delicate interplay takes place
between the stronger attraction, which is produced by the presence of the
deuterons, and the repulsion determined by the increased stiffness of the EoS.
This is a result of the modification introduced in the strengths of the
$\sigma$ and $\omega$ meson couplings. A global repulsive contribution is seen
at large baryon densities, while a significant reduction of the energy per
nucleon might be observed up to approximately $3n_{0}$, depending on the value
reached by the deuteron mass fraction in correspondence of the local maximum
observed in Fig. 7. However, as highlighted in the inset of Fig. 8, the three
black curves converge in the asymptotic limit, where a smooth disappearance of
the clusters was depicted in panel (b) of Fig. 6.

Figure 9:  Incompressibility $K$ as a function of the baryon density as
determined by employing four selected set of parameters listed in Table 2 and
for the mass shift parameterization given in Eqs. (58)-(61). For comparison,
the curve obtained for the DD2 parameterization in the deuteron-free case is
also shown as the red full line. The inset shows the differences with respect
to the predictions of the standard DD2 parameterization in a larger baryon
density range.

The incompressibility characterizes the curvature of the energy per nucleon.
It is defined here as

$K(n_{b})=9n_{b}^{2}\dfrac{\partial^{2}(E/A)}{\partial n_{b}^{2}}$ (64)

through a second derivative with respect to the baryon density.111The original
definition of $K$ uses a second derivative with respect to the Fermi momentum
and gives different results as compared to the definition used here but can
not be used at finite temperatures. However, both definitions coincide at
saturation. We recall that the incompressibility was not constrained within
the approach adopted in this work. A constraint on the incompressibility would
translate to a constraint on the second density derivative of the mass shift
at saturation. In addition, it would require, as a further input at
saturation, the knowledge of the density derivative of the deuteron mass
fraction. Although some numerical analyses have suggested that SRC pairs have
a minimum in the neighbourhood of the saturation, owing to the interplay
between the tensor component and the repulsive core of the nuclear force Rios
_et al._ (2014); Yang  _et al._ (2019); Li and Schulze (2016), we preferred to
prescind from applying such a constraint. Instead, the predictions for the
density behavior of $K$ are numerically extracted from the energy per nucleon.
They are plotted in Fig. 9 for the four selected parameterizations of Table 2
considered before. The related inset shows the differences with respect to the
standard DD2 reference, which is also shown in the main plot for comparison.
As a general feature, one observes that a softening of the EoS is recovered in
the region beyond saturation, up to a density around $3n_{0}$. Furthermore,
the size of this effect depends on the magnitude of the deuteron mass
fraction.

Table 3:  Values of the incompressibility $K_{0}$, in MeV, at saturation density as derived according to Eq. (64), for four selected parameterizations employed in this work. The result for the DD2 is also given for comparison.  | DD2 | DD2-d3 | DD2-$\chi$d1 | DD2-$\chi$d2 | DD2-$\chi$d3
---|---|---|---|---|---
$K_{0}$ [MeV] | $242.7$ | $199.6$ | $185.3$ | $207.3$ | $240.3$

As clearly evidenced in Table 3, the predictions for the incompressibility at
saturation lie within or below the range of values generally assumed for this
quantity Blaizot (1980); Youngblood  _et al._ (1999); Shlomo  _et al._ (2006);
Stone  _et al._ (2014) for all the parameterizations here employed. Strong
differences are observed in the high-density region, depending on the value
for the scaling factor $\chi$. A much stiffer EoS is envisaged in particular
at very large baryon densities, when a scaling factor $\chi=1/\sqrt{2}$ is
assumed, as a consequence of the significant change introduced in the balance
between the scalar and the vector components in this case.

Fig. 9 also exhibits another aspect to be discussed. The blue curve, which
corresponds to the DD2-d3 parameterization, reveals the emergence of a
discontinuity at a baryon density around 0.45 fm-3. This striking feature
signals the abrupt disappearance of the cluster, which was observed in Fig. 6,
panel (b) and already discussed before. It is worthwhile to mention that a
discontinuity in the matter incompressibility or in any other quantity related
to the second derivative of a thermodynamic potential is the signature for the
possible emergence of a second-order phase transition. One notices, by the
way, that this feature is in complete analogy to the disappearance of pairing
correlations that was observed in previous works at low density Burrello  _et
al._ (2014, 2016). Another discountinuity would emerge moreover at larger
densities, if the calculation with the DD2-d3 parameterization is not stopped
when the cluster dissolves, so that a further reappearance of the deuterons at
higher density is allowed. The inset of Fig. 9 shows that no discontinuity is
instead observed when the deuteron mass fraction smoothly decreases with the
density. The latter situation may however occur only for the parameterizations
characterized by a reduced value of the scaling factor $\chi$.

###  V.2 Predictions for ANM

Figure 10:  Panel(a): Deuteron mass fraction $X_{d}$ as function of the baryon
density, as determined by employing the DD2-$\chi$d1 parameterization. The
results for SNM ($\beta=0$) are compared with the corresponding ones deduced
with different values of the asymmetry $|\beta|$. Panel (b): Energy per
nucleon as a function of the baryon density, as determined by employing the
same parameterization and the same asymmetry values as in panel (a). The inset
shows the difference of the energy per nucleon between the DD2-$\chi$d1 and
DD2 parameterizations for the same asymmetry values, in a wider range of
baryon densities.

####  V.2.1 Deuteron mass fraction and EoS

Let us finally concentrate on the predictions for ANM. In Fig. 10, panel (a),
the density dependence of the deuteron mass fraction is plotted for different
values of the isospin asymmetry $|\beta|$. The DD2-$\chi$d1 parameterization
is employed for sake of illustration. The adopted parameterization allows one
to get the largest value for the deuteron mass fraction around the local
maximum, which was observed beyond 0.5 fm-3 in panel (b) of Fig. 7. In such a
way, the effect of embedding the quasi-deuterons at supra-saturation densities
is better emphasized. However, similar results, at least from a qualitative
point of view, would be obtained with the other sets of parameters accounting
for the presence of the deuterons. The red line indicates the maximum allowed
deuteron mass fraction values, which are compatible with a non negative value
of the Dirac effective mass of the nucleons. Let us recall that this curve
corresponds to $X_{d}=\mbox{min}\\{1-|\beta|,X_{d}^{\rm(max)}\\}$, so that it
evolves with $|\beta|$. Then only the border for SNM ($\beta=0$) is plotted in
panel (a) of Fig. 10 to avoid to overload the figure.

As a quite interesting result, Fig. 10 highlights that, although the mass
shift function defined in Eqs. (58)-(61) has no explicit dependence on the
isospin asymmetry, the corresponding deuteron mass fraction evolves with
$|\beta|$, giving rise to a continuous overall reduction when increasing the
neutron-proton asymmetry of the matter. Moreover, the smooth transition to the
cluster-free matter realized in SNM with the DD2-$\chi$d1 is preserved also in
the ANM case.

The density behavior of the energy per nucleon is depicted in panel (b) of
Fig. 10, for the same parameterization and the same asymmetry values
considered in panel (a). The inset of panel (b) displays moreover the
difference in the energy per nucleon, as determined with the DD2-$\chi$d1 and
the DD2 parameterization. Although not clearly visible, differently than in
the deuteron-free DD2 case, the zero-density limit does not approach zero,
except for the PNM case ($\beta=1$) in which deuterons are obviously not
formed. This feature will be more visible below, when studying the symmetry
energy.

The results shown in the inset of panel (b) help to disentangle the effect
induced on the stiffness of the EoS, owing to the rescaling of the meson
couplings at saturation and the changes ascribable to the presence of the
deuterons. The curve related to $|\beta|=1$, being characterized by $X_{d}=0$,
demonstrates that, apart from a tiny enhancement of the attraction below
saturation (not clearly visible in the figure), a much more repulsive PNM EoS
is produced for DD2-$\chi$d1 beyond saturation as compared to DD2. This is the
result of the modification induced on the effective interaction by changing
the meson coupling strengths. However, for the curves characterized by smaller
asymmetry values, such a repulsion is counterbalanced by the attraction
produced by the deuterons. An interplay analogous to the one discussed in the
SNM case takes place. In such a way, a reduction of the energy per nucleon of
ANM might be observed with the DD2-$\chi$d1 parameterization, in the
intermediate density region beyond saturation. This region may actually extend
up to very large densities, close to $n_{b}=0.9$ fm-3 in SNM. On the other
hand, in the asymptotic limit the quasi-deuterons tend to dissolve, so that
their extra-binding vanishes and the black curves depicted in the inset
converge to the $|\beta|=1$ one.

Figure 11:  Left panels: Deuteron mass fraction as function of the isospin
asymmetry $\beta$, for four selected parameterizations accounting for the
presence of the deuterons considered in this work. Right panels: Energy per
nucleon as a function of $\beta$ as determined by employing the same
parameterizations as in the left panels. For comparison, the curve obtained
for the DD2 parameterization in the deuteron-free case is also shown. Three
different values of the total baryon density are considered:
$n_{b}^{\rm(low)}=10^{-4}$ fm-3 (panels (a) and (d)), $n_{0}$ (panels (b) and
(e)) and $n_{b}^{\rm(high)}=10^{0}$ fm-3 (panels (c) and (f)).

A further insight into the dependence of the deuteron mass fraction and the
energy per nucleon on the isospin asymmetry might be achieved by looking at
Fig. 11. Three different values of the total baryon density are considered:
$n_{b}^{\rm(low)}=10^{-4}$ fm-3 (panels (a) and (d)), $n_{0}$ (panels (b) and
(e)) and $n_{b}^{\rm(high)}=10^{0}$ fm-3 (panels (c) and (f)). The standard
DD2 parameterization is also plotted in the right panels.

First of all, one observes that all the quantities are symmetric with respect
to the SNM ($\beta=0$) case. At the lowest density value considered in Fig.
11, $n_{b}^{\rm(low)}$, the deuteron fraction $X_{d}$ equals the maximum
allowed fraction and behaves thus like $X_{d}=1-|\beta|$, for the
parameterizations accounting for the presence of the deuterons. As a result, a
characteristic triangular shape of $X_{d}$ and $E/A$ is observed in panels (a)
and (d). The energy per nucleon does not follow the standard parabolic law
which is predicted in the deuteron-free DD2 case and reaches smaller values in
SNM. It approaches half of the deuteron binding energy in vacuum in the zero-
density limit. The energy per nucleon determined with the parameterizations
accounting for the deuterons coincides with the DD2 result only for matter
composed exclusively of neutrons or protons.

Secondly, a different picture is observed at the saturation density $n_{0}$.
There, by varying the isospin asymmetry, the deuteron mass fraction departs
from the value $X_{d}=0.2$ imposed for $\beta=0$ (see panel (b)). A mild
dependence of the deuteron mass fraction on $\beta$ is assessed for the
parameterization characterized by a reduced scaling factor $\chi$. A larger
sensitivity exists in the case of the DD2-d3 parameterization. Let us recall
that, for this parameterization, the clusters disappear for SNM at density
around 0.45 fm-3. Panel (b) shows that the clusters may dissolve already at
saturation density for finite $\beta$ values, at least in the case of the
DD2-d3 parameterization. On the other hand, the presence of the deuteron
persists at $n_{0}$ for all asymmetry values except $|\beta|=1$, for
parameterizations with $\chi=1/\sqrt{2}$. However, the parabolic dependence of
the energy per nucleon on the isospin asymmetry of the DD2 parameterization is
perfectly reproduced with all the considered parameterizations. This result is
clearly shown in panel (e). It originates from the requirement to keep the
energy per nucleon at saturation constrained, both for SNM and for matter
composed exclusively of neutrons or protons.

Thirdly, it is interesting to discuss what happens at the highest density
value, $n_{b}^{\rm(high)}$, considered in Fig. 11, panels (c) and (f). Here,
different results are obtained among the parameterizations characterized by a
scaling factor $\chi=1/\sqrt{2}$. However, for this density value, the
deuterons survive at all asymmetries only with the DD2-$\chi$d1
parameterization, for which the largest value was already predicted in the SNM
case. For the other two parameterizations with $\chi=1/\sqrt{2}$, the clusters
dissolve already for $|\beta|$ values smaller than 1. On the other hand, since
$n_{b}^{\rm(high)}$ lies beyond the density at which the cluster dissolution
is predicted in SNM, the deuteron mass fraction identically vanishes in case
of the DD2-d3 parameterization. The corresponding asymmetry dependence of the
energy per nucleon, which is depicted in panel (f), is then driven only by the
modification in the coupling strengths which was needed to keep the saturation
properties well constrained. As in Fig. 8, a slightly larger attraction is
foreseen with the DD2-d3 parameterization, with respect to the DD2 reference
case. The opposite happens instead when the parameterizations with a reduced
value of the scaling factor are considered. In this case, a stronger repulsion
is envisaged, partially mitigated, at least for small asymmetry values, by the
stronger binding provided by the presence of deuterons. Quite interestingly,
one observes that the parameterizations plotted by black curves always
converge when approaching $|\beta|=1$, where $X_{d}=0$. A change of the
parameterization, which implies an according change of the deuteron fraction
at supra-saturation densities, affects the curvature of the dependence of the
energy per nucleon on $\beta$ and thus the symmetry energy. The latter
quantity will be studied in detail below.

####  V.2.2 Symmetry energy and its slope

Figure 12:  Panel (a): Symmetry energy $J$ as function of the baryon density,
as determined through Eq. (65) for four selected parameterizations employed in
this work. The inset shows a zoom at sub-saturation densities. Panel (b):
Slope $L$ of the symmetry energy as function of the baryon density, as
determined through Eq. (67), for the same parameterizations as in panel (a).
The inset shows a zoom around the saturation density $n_{0}$.

In the present work, the symmetry energy $J$ is calculated as the difference
between the energies per nucleon in PNM and SNM

$J(n_{b})=\left.\dfrac{E}{A}\right|_{\beta=1}(n_{b})-\left.\dfrac{E}{A}\right|_{\beta=0}(n_{b})\>.$
(65)

The quantity obtained through this equation is identical to the symmetry
energy calculated from the usual definition

$J(n_{b})=\frac{1}{2}\left.\frac{\partial^{2}(E/A)}{\partial\beta^{2}}\right|_{\beta=0}$
(66)

using a second derivative of the energy per nucleon with respect to the
asymmetry, if $E/A$ follows a quadratic dependence on $\beta$. The density
dependence of $J$ is plotted in Fig. 12, panel (a). The inset of panel (a)
shows a zoom at sub-saturation densities. Once again, the inset highlights the
dissimilar behavior in the zero density limit, when the presence of the
clusters is taken into account or not. It reflects the differences existing in
the very dilute regime of the SNM EoS with respect to the deuteron-free case.
If clustering is taken into account, the symmetry energy approaches indeed
half of the deuteron binding energy in the zero-density limit in contrast to
the simple description without explicit two-particle correlations. Some
differences emerge among the parameterizations which account for the presence
of the deuterons below $n_{0}$. Apart from the constraint at saturation, no
restrictions have been imposed on the density behavior of $J$.

Concerning the behavior beyond saturation density, despite the presence of the
deuterons, at least up to $0.45$ fm-3, the blue curve remains close to the
result of the standard DD2 parameterization. Huge differences are instead
observed in the supra-saturation density region when the parameterizations
with a reduced value of the deuteron-meson coupling scaling factor is
considered. Furthermore, the size of the effect depends quite strongly on the
mass fraction of the deuterons. The differences vanish in the asymptotic
limit, where the deuterons disappear. The black and blue curves approach the
result of the standard DD2 parameterization as the contribution of the $\rho$
meson to the symmetry energy vanishes due to the suppression of its coupling,
leaving the imbalance of the Fermi momenta of the nucleons as the main
contribution to $J$.

Finally, it is interesting to look at the density derivative of the curves
plotted in panel (a). The slope $L$ of the symmetry energy is numerically
calculated here as

$L(n_{b})=3n_{b}\dfrac{dJ}{dn_{b}}$ (67)

and the result is shown in Fig. 12, panel (b). Except for the different
behavior in the zero-density limit, the blue curve roughly coincides with the
standard DD2 in the whole range of densities. A small kink is only observed
for the DD2-d3 parameterization around 0.45 fm-3. This kink is related to the
disappearance of the deuterons in the SNM case. We recall that this feature
was also responsible for the emergence of the discontinuity in the matter
incompressibility discussed before. A huge double oscillation around the
result obtained with the standard DD2 parameterization is observed for the
parameterizations with a reduced scaling factor $\chi$. The magnitude of this
oscillation depends again on the deuteron mass fraction predicted at supra-
saturation densities, thus reflecting the result shown in panel (a) of Fig.
12. The inset of panel (b) displays the predictions for the symmetry energy
slope around saturation. In spite of the extremely large differences in the
high-density regime, reasonable values are obtained at saturation density with
all the parameterizations considered in this work. These values, which are
collected in Table 4, lie within the range usually assumed for the slope of
the symmetry energy, see, e.g., Li  _et al._ (2021) and references therein.
The three black curves, crossing each other at saturation, naturally provide
the same value. Alternative scenarios manifest for the high-density behavior
of the symmetry energy and its slope. The stiffness of the EoS in the supra-
saturation density regime turns out to be strongly affected by the value of
the scaling factor. It is less dependent on the mere presence of the
deuterons. Thus, as a general feature, one concludes that smaller values of
the scaling factor correspond to higher stiffness values in the density region
beyond saturation relevant for the applications of the model.

Table 4:  Values of the slope of the symmetry energy $L_{0}$, in MeV, at saturation density as derived according to Eq. (67), for four selected parameterizations employed in this work. The result for the DD2 is also given for comparison.  | DD2 | DD2-d3 | DD2-$\chi$d1 | DD2-$\chi$d2 | DD2-$\chi$d3
---|---|---|---|---|---
$L_{0}$ [MeV] | $57.94$ | $56.49$ | $67.50$ | $67.50$ | $67.50$

##  VI Conclusions and outlook

In this paper, we have proposed and explored a novel approach to embed SRCs
within the GRDF, a well-established phenomenological EDF based on nucleon and
cluster degrees of freedom. In such a way, we aim to overcome the
inconsistencies between recent experimental evidences, which brought to light
the existence of SRCs, and the predictions of phenomenological models derived
from mean-field approaches without explicit correlations at densities around
saturation. Previous generalisations of these EDFs represent many-body
correlation by cluster which dissolve, by construction, when the nuclear
saturation density is approached from below. Within an extended relativistic
mean-field model with density dependent couplings, the idea of this work was
to effectively account for the existence of SRCs through proper in-medium
modifications of the cluster properties around saturation and above. They are
considered as quasi-particles with density-dependent binding energies. Quasi-
deuterons immersed in dense matter are used as surrogate for correlations in
this first exploratory step. For the time being, the zero temperature case,
where the deuteron fraction is determined by the density of a boson
condensate, was addressed.

Suitable parameterizations of the cluster mass shift were derived, for the
first time, for all baryon densities. The proposed mass shift functions comply
with the available constraints and were employed to determine the density
dependence of the quasi-deuteron mass fraction at arbitrary isospin
asymmetries, thus for symmetric as well as for asymmetric nuclear matter. They
were constrained by microscopic many-body calculations in the low-density
limit, by specifying the actual deuteron fraction at saturation and by
assuming a deuteron mass fraction behavior that respects the boson
condensation condition at higher densities. The effective deuteron fraction
around saturation was specified by extrapolating the experimental results on
SRC pairs in nuclei to infinite nuclear matter. Further constraints were
moreover imposed at supra-saturation densities by the maximum allowed deuteron
fraction, compatible with a non-negative value of the Dirac effective mass of
the nucleons.

A proper description of well-constrained nuclear matter quantities at
saturation required a refit of the nucleon-meson coupling strengths. An
important role of the coupling scaling factor $\chi$ was revealed. It rules
the coupling strength of the mesons with nucleons bound in the clusters. Such
a scaling factor plays actually a primary role in the whole analysis developed
in this work.

The natural choice was supposed to be that the nucleons inside the deuterons
couple to the mesons with the same strength as the unbound nucleons. However,
with this choice, the deuteron mass shift is not a monotonic function of the
deuteron mass fraction for all baryon densities. As a result, the relation
between the mass shift and the deuteron mass fraction can not be inverted
uniquely in all cases. The same holds for any scaling factor $\chi$ larger
than $1/\sqrt{2}$. Thus, as possible extreme values of two well distinct
behaviors, two different values of the scaling factor, namely $\chi=1$ and
$\chi=1/\sqrt{2}$ were chosen, in the calculations performed in this work. The
latter value is however significantly smaller than the universal scaling
factor for the cluster-meson coupling strength. A value smaller than $\chi=1$
was proposed in previous calculations of the EoS to take into account in-
medium effects and to get a good description of the chemical equilibrium
constants determined from recent experimental data.

As a general feature, our analysis shows that, for $\chi=1$, the only possible
smooth solution for the density dependence of the deuteron mass fraction
implies a sudden disappearance of the clusters at a density below the one
corresponding to the emergence of a pole. In correspondence of this density, a
discontinuity in the matter incompressibility emerges, analogous to the one
observed at low density, owing to the disappearance of the pairing
correlations and indicating the emergence of a second order phase transition.
The analogy between the behavior of pairing and SRCs deserves however further
investigation. For $\chi=1$, the density where the pole emerges is located
around three times the saturation density. However, for more realistic values
of the scaling factor, the pole is expected to appear at much higher
densities, thus much beyond the range that is relevant in applications of the
model. When the scaling factor value $\chi=1/\sqrt{2}$ is considered,
alternative solutions exist, permitting smooth functions of the density
dependence of the deuteron mass fraction for all densities. Three different
parameterizations, providing such a smooth behavior and characterized by
different maximum deuteron fraction values at supra-saturation densities, were
proposed.

Striking effects on some thermodynamic quantities are recognized, owing to the
presence of the quasi-deuterons in the neighbourhood of saturation and at
supra-saturation densities. In particular, a softening of the SNM EoS is
systematically observed with respect to the standard DD2 parameterization,
which does not include deuteron-like correlations. However, the stronger
attraction, which is produced by the presence of the deuterons, might be
counterbalanced by the repulsion driven by the modified coupling strengths.
This delicate interplay is additionally tuned by the value of the scaling
factor, which determines then alternative scenarios for the high-density
behavior of the symmetry energy and its slope. In general, one concludes that
smaller values of the scaling factor correspond to higher stiffness of the EoS
in the supra-saturation density regime.

Last, but not least, it is worthwhile to recall that our analysis permits to
also recover the correct low-density limit of the EoS. Indeed, at zero-
density, both the energy per nucleon of SNM and the symmetry energy tend to be
equal to one half of the deuteron binding energy in vacuum, in contrast to the
predictions of standard mean-field models without cluster correlations.

The findings of the present study represent a first step to improve the
description of nuclear matter and its EoS at supra-saturation densities in
EDFs by considering correlations in an effective way. In a next step, the
single-particle momentum distributions can be explored using proper wave
functions of the quasi-deuteron in the medium. They have to be derived
consistently with the interaction used in the model and will lead to
prediction of the cluster mass shifts and fractions that can be compared to
the suggested forms of the present work. The many-body wave function of a
cluster contains correlated nucleons with a specific momentum distribution.
Then an imprint on the single-nucleon momentum distribution in nuclear matter
is expected, such that a high-momentum tail develops even at zero temperature,
as observed in the experimental study of SRCs by nucleon knockout with high-
energy electrons.

The present approach can be generalized to finite temperatures, where a
further change of the single-nucleon momentum distribution arises owing to the
thermal change in the distribution functions. Also a momentum dependence of
the mass shift and a more involved dependence on the isospin asymmetry might
be considered in a future work, together with the effect of including heavier
clusters and to investigate their relative importance.

As a perspective, we finally aim at investigating the effect of SRCs on
neutron stars in the EDF framework, similarly to what was done in some prior
studies, see, e.g., Lu  _et al._ (2021b), but our approach is to replace
heuristic parameterizations of the momentum distributions with more
microscopically founded descriptions.

More in general, we aim at achieving a more comprehensive description of
correlations and clustering phenomena, which represents still a challenge from
a theoretical point of view, despite the importance of these features in the
widest scope of astrophysical applications and for general aspects of
reactions dynamics in heavy-ion collisions.

## Acknowledgments

The authors thank Maria Colonna and Gerd Röpke for their comments and
suggestions on this work. S. B. acknowledges support from the Alexander von
Humboldt foundation.

##  Appendix A Mass fraction derivative of the deuteron mass shift

The deuteron mass shift (III.2.1) is in general a function of the baryon
density $n_{b}$, the asymmetry $\beta$, the deuteron fraction $X_{d}$, and the
temperature T. In this section, the derivative of $\Delta m_{d}^{\rm(high)}$
with respect to $X_{d}$ is derived for constant $n_{b}$ and $\beta$ at $T=0$.
In a first step, the derivative of the nucleon effective chemical potential
with respect to $X_{d}$ is expressed as

$\left.\frac{\partial\mu_{q}^{\ast}}{\partial
X_{d}}\right|_{n_{b},\beta}=\frac{1}{\mu_{q}^{\ast}}\left(k_{q}\left.\frac{\partial
k_{q}}{\partial X_{d}}\right|_{n_{b},\beta}+m_{q}^{\ast}\left.\frac{\partial
m_{q}^{\ast}}{\partial X_{d}}\right|_{n_{b},\beta}\right)$ (68)

using Eq. (35). The derivative of the Fermi momentum of a nucleon $q=n,p$ is
found with help of the relation

$\left.\frac{\partial n_{q}^{(v)}}{\partial
X_{d}}\right|_{n_{b},\beta}=-\frac{n_{b}}{2}=\frac{3n_{q}^{(v)}}{k_{q}}\left.\frac{\partial
k_{q}}{\partial X_{d}}\right|_{n_{b},\beta}$ (69)

for the vector density (19). Using Eqs. (68) and (69), the derivative of the
scalar density of the nucleons can be written as

$\left.\frac{\partial n_{q}^{(s)}}{\partial
X_{d}}\right|_{n_{b},\beta}=-\frac{n_{b}}{2}\frac{m_{q}^{\ast}}{\mu_{q}^{\ast}}+f_{q}\left.\frac{\partial
m_{q}^{\ast}}{\partial X_{d}}\right|_{n_{b},\beta}$ (70)

with the factor

$f_{q}=3\left(\frac{n_{q}^{(s)}}{m_{q}^{\ast}}-\frac{n_{q}^{(v)}}{\mu_{q}^{\ast}}\right)$
(71)

after several steps of recasting the individual contributions. Eq. (70)
contains again the derivative of the effective mass that assumes the simple
form

$\left.\frac{\partial m_{q}^{\ast}}{\partial
X_{d}}\right|_{n_{b},\beta}=-C_{\sigma}\left.\frac{\partial
n_{\sigma}}{\partial X_{d}}\right|_{n_{b},\beta}$ (72)

because $C_{\sigma}$ depends only on $n_{b}$. With the derivatives

$\left.\frac{\partial n_{d}^{(v)}}{\partial
X_{d}}\right|_{n_{b},\beta}=\left.\frac{\partial n_{d}^{(s)}}{\partial
X_{d}}\right|_{n_{b},\beta}=\frac{n_{b}}{2}$ (73)

of the deuteron densities, the derivative of the $\sigma$ meson source density
(11) is found as

$\displaystyle\left.\frac{\partial n_{\sigma}}{\partial
X_{d}}\right|_{n_{b},\beta}$ $\displaystyle=$
$\displaystyle\frac{n_{b}}{1+\left(f_{n}+f_{p}\right)C_{\sigma}}\>\mathcal{U}_{d}$
(74)

with

$\displaystyle\mathcal{U}_{d}$ $\displaystyle=$
$\displaystyle\chi_{d\sigma}-\frac{m_{\rm
nuc}^{\ast}}{2}\left(\frac{1}{\mu_{n}^{\ast}}+\frac{1}{\mu_{p}^{\ast}}\right)\>$
(75)

whereas

$\left.\frac{\partial n_{\omega}}{\partial
X_{d}}\right|_{n_{b},\beta}=-n_{b}(1-\chi_{d\omega})$ (76)

for the source density of the $\omega$ meson. Finally, the derivative of the
mass shift with respect to the mass fraction is obtained in the compact form

$\displaystyle\left.\frac{\partial\Delta m_{d}^{\mathrm{(high)}}}{\partial
X_{d}}\right|_{n_{b},\beta}$ $\displaystyle=$
$\displaystyle\left[\frac{2C_{\sigma}}{1+(f_{n}+f_{p})C_{\sigma}}\mathcal{U}_{d}^{2}-\frac{\pi^{2}}{2\mu_{n}^{\ast}k_{n}}-\frac{\pi^{2}}{2\mu_{p}^{\ast}k_{p}}\right]n_{b}$

with a contribution from the $\sigma$ meson and kinetic terms.

##  Appendix B Density derivative of deuteron mass shift

In this section, the derivative of $\Delta m_{d}^{\rm(high)}$ with respect to
$n_{b}$ is derived for an arbitrary function $X_{d}(n_{b})$ and constant
$\beta$ at $T=0$. It requires again several steps. First, the derivatives of
the source densities have to be determined. For the $\omega$ meson one finds

$\left.\frac{\partial n_{\omega}}{\partial
n_{b}}\right|_{\beta}=1-\left(1-\chi_{d\omega}\right)Y_{d}$ (78)

with the quantity

$Y_{d}=\left.\frac{\partial(n_{b}X_{d})}{\partial
n_{b}}\right|_{\beta}=X_{d}+n_{b}\left.\frac{\partial X_{d}}{\partial
n_{b}}\right|_{\beta}$ (79)

that contains the derivative of the deuteron mass fraction. For the $\sigma$
meson, the calculation is more involved. Here, the relations

$\displaystyle\left.\frac{\partial n_{\sigma}}{\partial n_{b}}\right|_{\beta}$
$\displaystyle=$ $\displaystyle\left(f_{n}+f_{p}\right)\left.\frac{\partial
m_{\mathrm{nuc}}^{\ast}}{\partial n_{b}}\right|_{\beta}+\chi_{d\sigma}Y_{d}$
$\displaystyle+\frac{m_{\mathrm{nuc}}^{\ast}}{2\mu_{n}^{\ast}}\left(1+\beta-
Y_{d}\right)+\frac{m_{\mathrm{nuc}}^{\ast}}{2\mu_{p}^{\ast}}\left(1-\beta-
Y_{d}\right)$

with the factor $f_{q}$ defined in Eq. (71) and

$\left.\frac{\partial m_{\rm nuc}^{\ast}}{\partial
n_{b}}\right|_{\beta}=-C_{\sigma}\left.\frac{\partial n_{\sigma}}{\partial
n_{b}}\right|_{\beta}-C_{\sigma}^{\prime}n_{\sigma}$ (81)

for the derivative of the effective nucleon mass can be combined to obtain the
form

$\displaystyle\left.\frac{\partial n_{\sigma}}{\partial n_{b}}\right|_{\beta}$
$\displaystyle=$
$\displaystyle\left[1+\left(f_{n}+f_{p}\right)C_{\sigma}\right]^{-1}$
$\displaystyle\left[-\left(f_{n}+f_{p}\right)C_{\sigma}^{\prime}n_{\sigma}+\chi_{d\sigma}Y_{d}\right.$
$\displaystyle\left.+\frac{m_{\mathrm{nuc}}^{\ast}}{2\mu_{n}^{\ast}}\left(1+\beta-
Y_{d}\right)+\frac{m_{\mathrm{nuc}}^{\ast}}{2\mu_{p}^{\ast}}\left(1-\beta-
Y_{d}\right)\right]$

with an explicit dependence on $Y_{d}$. In the next step, the derivative of
the mass shift assumes the form

$\displaystyle\left.\frac{\partial\Delta m_{d}^{(\mathrm{high})}}{\partial
n_{b}}\right|_{\beta}$ $\displaystyle=$
$\displaystyle\frac{k_{n}}{\mu_{n}^{\ast}}\left.\frac{\partial k_{n}}{\partial
n_{b}}\right|_{\beta}+\frac{k_{p}}{\mu_{p}^{\ast}}\left.\frac{\partial
k_{p}}{\partial
n_{b}}\right|_{\beta}+\left(\frac{1}{\mu_{n}^{\ast}}+\frac{1}{\mu_{p}^{\ast}}\right)m_{\rm
nuc}^{\ast}\left.\frac{\partial m_{\rm nuc}^{\ast}}{\partial
n_{b}}\right|_{\beta}$ $\displaystyle+2(1-\chi_{d\omega})\left(C_{\omega}\
\left.\frac{\partial n_{\omega}}{\partial
n_{b}}\right|_{\beta}+C_{\omega}^{\prime}n_{\omega}\right)$
$\displaystyle+2\chi_{d\sigma}\left(C_{\sigma}\left.\frac{\partial
n_{\sigma}}{\partial
n_{b}}\right|_{\beta}+C_{\sigma}^{\prime}n_{\sigma}\right)$

with

$\left.\frac{\partial k_{n}}{\partial
n_{b}}\right|_{\beta}=\frac{k_{n}}{3n_{n}^{(v)}}\left(\frac{1+\beta-
Y_{d}}{2}\right)$ (84)

and

$\left.\frac{\partial k_{p}}{\partial
n_{b}}\right|_{\beta}=\frac{k_{p}}{3n_{p}^{(v)}}\left(\frac{1-\beta-
Y_{d}}{2}\right)\>.$ (85)

Using the expressions (78) and (B), the final result can be expressed in
compact form as

$\left.\frac{\partial\Delta m_{d}^{(\mathrm{high})}}{\partial
n_{b}}\right|_{\beta}=\mathcal{W}_{d}-\mathcal{Z}_{d}Y_{d}$ (86)

with the auxiliary quantities

$\displaystyle\mathcal{Z}_{d}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{2\mu_{n}^{\ast}k_{n}}+\frac{\pi^{2}}{2\mu_{p}^{\ast}k_{p}}+2\left(1-\chi_{d\omega}\right)^{2}C_{\omega}$
$\displaystyle-\frac{2C_{\sigma}}{1+\left(f_{n}+f_{p}\right)C_{\sigma}}\>\mathcal{U}_{d}^{2}$
$\displaystyle\mathcal{W}_{d}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{2\mu_{n}^{\ast}k_{n}}\left(1+\beta\right)+\frac{\pi^{2}}{2\mu_{p}^{\ast}k_{p}}\left(1-\beta\right)$
$\displaystyle+2(1-\chi_{d\omega})\left(C_{\omega}+C_{\omega}^{\prime}n_{\omega}\right)$
$\displaystyle+\frac{2}{1+\left(f_{n}+f_{p}\right)C_{\sigma}}\>\mathcal{U}_{d}$
$\displaystyle\left[C_{\sigma}^{\prime}n_{\sigma}+C_{\sigma}\frac{m_{\mathrm{nuc}}^{\ast}}{2}\left(\frac{1+\beta}{\mu_{n}^{\ast}}+\frac{1-\beta}{\mu_{p}^{\ast}}\right)\right]$

and $\mathcal{U}_{d}$ as given in (75).

##  Appendix C Conversion of parameters at saturation

In order to find the coupling strengths $\Gamma_{\sigma,0}$ and
$\Gamma_{\omega,0}$ as well as the deuteron mass shift and its density
derivative at saturation, a step-by-step procedure can be followed. These
quantities are determined as soon as the saturation density $n_{0}$, the
binding energy per nucleon $B_{0}$, the effective nucleon mass
$m^{\ast}_{\mathrm{nuc},0}$ and deuteron fraction $X_{d,0}$ of SNM are
specified.

In a first step, the scalar and vector densities

$n_{d,0}^{(s)}=n_{d,0}^{(v)}=n_{0}\frac{X_{d,0}}{2}$ (89)

of the deuteron and the total vector density

$n_{\mathrm{nuc},0}^{(v)}=n_{0}\left(1-X_{d,0}\right)$ (90)

of the nucleons are immediately obtained from $n_{0}$ and $X_{d,0}$ in SNM.
Then the Fermi momentum

$k_{\mathrm{nuc},0}=\left[\frac{6\pi^{2}}{g_{\mathrm{nuc}}}n_{\mathrm{nuc},0}^{(v)}\right]^{1/3}$
(91)

with degeneracy factor $g_{\mathrm{nuc}}=4$ and the effective chemical
potential

$\mu_{\mathrm{nuc},0}^{\ast}=\sqrt{k_{\mathrm{nuc},0}^{2}+\left(m^{\ast}_{\mathrm{nuc},0}\right)^{2}}$
(92)

allow to calculate the scalar density

$\displaystyle n^{(s)}_{\mathrm{nuc},0}$ $\displaystyle=$
$\displaystyle\frac{g_{\mathrm{nuc}}m^{\ast}_{\mathrm{nuc},0}}{4\pi^{2}}\left[k_{\mathrm{nuc},0}\mu_{\mathrm{nuc},0}^{\ast}-\left(m^{\ast}_{\mathrm{nuc},0}\right)^{2}\ln\frac{k_{\mathrm{nuc},0}+\mu_{\mathrm{nuc},0}^{\ast}}{m^{\ast}_{\mathrm{nuc},0}}\right]$

using the effective nucleon mass $m_{\mathrm{nuc},0}^{\ast}$. Then the source
densities

$n_{\sigma,0}=n_{\mathrm{nuc}}^{(s)}+2\chi n_{d,0}^{(s)}$ (94)

and

$n_{\omega,0}=n_{\mathrm{nuc}}^{(v)}+2\chi n_{d,0}^{(v)}$ (95)

with the deuteron-meson coupling scaling factor $\chi$ are found and the
pressure contribution

$p_{\mathrm{nuc},0}=\frac{1}{4}\left[\mu_{\mathrm{nuc},0}^{\ast}n_{\mathrm{nuc},0}^{(v)}-m^{\ast}_{\mathrm{nuc},0}n_{\mathrm{nuc},0}^{(s)}\right]$
(96)

of the nucleons can be calculated immediately.

In the next step, the effective nucleon mass determines the scalar potential

$S_{\mathrm{nuc},0}=m_{\mathrm{nuc}}-m^{\ast}_{\mathrm{nuc},0}$ (97)

of the nucleons and thus the scalar coupling

$C_{\sigma,0}=\frac{S_{\mathrm{nuc},0}}{n_{\sigma,0}}$ (98)

and finally the coupling strength

$\Gamma_{\sigma,0}=m_{\sigma}\sqrt{C_{\sigma,0}}$ (99)

of the $\sigma$ meson. The binding energy per nucleon $B_{0}$ gives the
chemical potential

$\mu_{\mathrm{nuc},0}=m_{\mathrm{nuc}}-B_{0}$ (100)

at saturation and then the vector potential

$V_{\mathrm{nuc},0}=\mu_{\mathrm{nuc},0}-\mu_{\mathrm{nuc},0}^{\ast}$ (101)

of the nucleons. The latter quantity can be expressed in general as

$V_{\mathrm{nuc},0}=C_{\omega,0}n_{\omega,0}+C_{\rho,0}n_{\rho,0}+U_{0}^{(r)}$
(102)

with the auxiliary quantity

$U_{0}^{(r)}=V_{0}^{(r)}+W_{0}^{(r)}$ (103)

that also appears in the total pressure

$\displaystyle P_{0}=p_{\mathrm{nuc},0}$
$\displaystyle+\frac{1}{2}\left(C_{\omega,0}n_{\omega,0}^{2}+C_{\rho,0}n_{\rho,0}^{2}-C_{\sigma,0}n_{\sigma,0}^{2}\right)+U_{0}^{(r)}n_{0}\>.$

For SNM, however, the $\rho$-meson contribution does not appear, since
$n_{\rho}$ is identically zero. The two equations (102) and (C) allow to solve
for the $\omega$ coupling

$\displaystyle C_{\omega,0}$ $\displaystyle=$
$\displaystyle\left(2n_{\omega,0}n_{0}-n_{\omega,0}^{2}\right)^{-1}$
$\displaystyle\left(2p_{\mathrm{nuc},0}+2V_{\mathrm{nuc},0}n_{0}-C_{\sigma,0}n_{\sigma,0}^{2}\right)$

and further the coupling strength

$\Gamma_{\omega,0}=m_{\omega}\sqrt{C_{\omega,0}}$ (106)

using $P_{0}=0$. With known $C_{\sigma,0}$ and $C_{\omega,0}$, their
derivatives $C_{\sigma,0}^{\prime}$ and $C_{\omega,0}^{\prime}$ can be
determined using the same functional density dependence of the couplings as in
the reference parameterization. Thus also the rearrangement contribution

$V_{0}^{(r)}=\frac{1}{2}\left(C_{\omega,0}^{\prime}n_{\omega,0}^{2}+C_{\rho,0}^{\prime}n_{\rho,0}^{2}-C_{\sigma,0}^{\prime}n_{\sigma,0}^{2}\right)$
(107)

is given.

Finally, from Eqs. (102) and (103) one finds

$W_{0}^{(r)}=V_{\mathrm{nuc},0}-C_{\omega,0}n_{\omega,0}-V_{0}^{(r)}$ (108)

and the deuteron mass shift derivative

$\left.\frac{d\Delta
m_{d}}{dn_{b}}\right|_{n_{0}}=\frac{W_{0}^{(r)}}{n_{d,0}^{(s)}}$ (109)

at saturation. The deuteron mass shift itself is determined as

$\displaystyle\Delta m_{d,0}$ $\displaystyle=$ $\displaystyle
B_{d}+2\left(\mu_{\mathrm{nuc},0}^{\ast}-m_{\mathrm{nuc},0}^{\ast}\right)$
$\displaystyle+2(1-\chi_{d\omega})C_{\omega}n_{\omega}-2(1-\chi_{d\sigma})C_{\sigma}n_{\sigma}$

from the condensation condition with the binding energy of the deuteron in
vacuum $B_{d}$.

The rescaling of the $\sigma$ and $\omega$ coupling strengths induces a
modification of the energy per nucleon in PNM at saturation and thus of the