On gravitational preheating

Oleg Lebedev and Jong-Hyun Yoon

Department of Physics and Helsinki Institute of Physics,
Gustaf Hällströmin katu 2a, FI-00014 Helsinki, Finland

Abstract

We consider dark matter production during the inflaton oscillation epoch. It
is conceivable that renormalizable interactions between dark matter and
inflaton may be negligible. In this case, the leading role is played by higher
dimensional operators generated by gravity and thus suppressed by the Planck
scale. We focus on dim-6 operators and study the corresponding particle
production in perturbative and non–perturbative regimes. We find that the dark
matter production rate is dominated by non–derivative operators involving
higher powers of the inflaton field. Even if they appear with small Wilson
coefficients, such operators can readily account for the correct dark matter
abundance.

###### Contents

  1. 1 Introduction
  2. 2 The set–up
  3. 3 Perturbative dark matter production
    1. 3.1 $\phi^{2}s^{2}$ interaction
    2. 3.2 $\phi^{2}(\partial_{\mu}s)^{2}$ interaction
    3. 3.3 $\phi^{4}s^{2}$ interaction
    4. 3.4 Relative efficiency
  4. 4 Resonant dark matter production via dim–6 operators
    1. 4.1 On resonant production via $\phi^{2}s^{4}$
  5. 5 Lattice simulations: reproducing the correct DM abundance
  6. 6 Conclusion

##  1 Introduction

The nature of dark matter (DM) remains an outstanding mystery of modern
physics. The null DM direct detection results motivate one to explore the
possibility that dark matter has feeble interactions, in which case it does
not reach thermal equilibrium with the environment. Therefore, its abundance
is sensitive to the production mechanism. One of such mechanisms is provided
by gravity, which can efficiently produce particles in non–adiabatic
environments.

In the absence of any non–gravitational couplings, the expansion of the
Universe is itself a source of particle production [1],[2]. For example, the
equation of motion (EOM) for a momentum mode $\chi_{k}$ of a free scalar in
the Friedmann Universe with the metric $ds^{2}=a^{2}(\eta)\,(d\eta^{2}-d{\bf
x}^{2})$ reads [3]

$\chi_{k}^{\prime\prime}+\omega_{k}^{2}\chi_{k}=0\;,$ (1)

where

$\omega_{k}^{2}=k^{2}-{a^{\prime\prime}\over a}(1-6\xi)+m_{\chi}^{2}a^{2}\;,$
(2)

$a$ is the scale factor, $\xi$ is the non–minimal coupling to gravity [4] and
the prime denotes differentiation with respect to conformal time $\eta$. Time
variation of $\omega_{k}$ is non–adiabatic if
$\omega^{\prime}_{k}/\omega_{k}^{2}\gtrsim 1$, which implies particle creation
due to expansion. For low $k$ and conformal coupling $\xi=1/6$, this is
equivalent to $a^{\prime}/a^{2}=H\gtrsim m_{\chi}$ such that particles lighter
than the Hubble rate $H$ are constantly created. The effect can be even
stronger for non–conformal $\xi$. The accumulated abundance of $\chi$ can
constitute dark matter [3, 5], depending on $m_{\chi}$ and its
self–interaction.

Such particle production can be viewed in terms of the scalar field condensate
$\langle\chi^{2}\rangle$. Light scalars are subject to quantum fluctuations of
order $H$ [6] so that a semi–classical field $\chi$ experiences a random walk.
As a result, a significant $\langle\chi^{2}\rangle$ can accumulate, for
example, by the end of inflation and play the role of dark matter [7],[8].
Again, this mechanism is purely gravitational. The consequent dark matter
distribution is not correlated with the inflaton fluctuations, therefore this
possibility is subject to strict isocurvature constraints [8].

In this work, we focus on other aspects of particle production due to
gravitational effects. Specifically, gravity is believed to generate couplings
between different sectors of the theory as long as these are consistent with
gauge symmetries. The corresponding operators may be non–renormalizable and
thus suppressed by the Planck scale. Nevertheless, they can play an important
role in dark matter production. This was recently emphasised in [9],[10],
where the effects due to tree level graviton exchange were considered.

We study dark matter production during the inflaton oscillation phase, which
sets in immediately after inflation and creates a non–adiabatic environment
[11, 12, 13]. Using the effective field theory approach, we focus on the
leading gravity–induced dim–6 operators assuming that the renormalizable
couplings between the inflaton and dark matter vanish. If DM is feebly
interacting, its eventual abundance is determined by the number of DM quanta
produced at this ‘‘preheating’’ stage. To this end, we identify the dominant
operator and study whether it can be responsible for the correct dark matter
abundance.

##  2 The set–up

Consider the possibility that the renormalizable couplings between the
inflaton $\phi$ and dark matter $s$ are zero or negligibly small. Then, the
$\phi-s$ interaction can be described by a series of higher dimensional
operators generated by gravity and thus suppressed by the Planck scale $M_{\rm
Pl}$. Let us assume for simplicity that these operators exhibit an approximate
$\phi\rightarrow-\phi$ symmetry such that the lowest operator dimension is
six:111The dim-5 operator $\phi^{3}s^{2}$ may in general be present, in which
case it would dominate dark matter production.

$\displaystyle\Delta{\cal L}_{6}={C_{1}\over M_{\rm
Pl}^{2}}\;(\partial_{\mu}\phi)^{2}s^{2}+{C_{2}\over M_{\rm
Pl}^{2}}\;(\phi\partial_{\mu}\phi)(s\partial^{\mu}s)+{C_{3}\over M_{\rm
Pl}^{2}}\;(\partial_{\mu}s)^{2}\phi^{2}-{C_{4}\over M_{\rm
Pl}^{2}}\;\phi^{4}s^{2}-{C_{5}\over M_{\rm Pl}^{2}}\;\phi^{2}s^{4}\;,$ (3)

where we have replaced the covariant derivatives with the partial ones. The
inflaton field with mass $m_{\phi}$ is assumed to have either (locally)
quadratic or quartic potential, while the dark matter mass $m_{s}$ is taken to
be negligible compared to the typical scales of the problem. Some of the above
operators such as $(\partial_{\mu}\phi)^{2}s^{2}$ and
$(\partial_{\mu}s)^{2}\phi^{2}$, along with ${m_{\phi,s}^{2}\over M_{\rm
Pl}^{2}}\,\phi^{2}s^{2}$, are generated by the tree level graviton exchange
[9]. Others can be generated at loop level and non–perturbatively. Since
gravity is non–renormalizable, their coefficients should be treated as
arbitrary input parameters. The above interactions are responsible for dark
matter production after inflation, in particular, during the inflaton
oscillation phase. Depending on the $C_{i}$ coefficients, the production
mechanism can be perturbative or non–perturbative (resonant).

To mention but one example, some of the above operators operators appear
automatically in theories with non–minimal couplings of scalars to gravity
[4]. In particular, $(\partial_{\mu}s)^{2}\phi^{2}$,
$(\partial_{\mu}\phi)^{2}s^{2}$ and ${s^{2}}\,V(\phi)$ are induced already at
tree level by the metric transformation from the Jordan frame to the Einstein
frame [14]. The (unsuppressed) operators $\phi^{2}s^{4}$ and $\phi^{4}s^{2}$
appear in these models at 1–loop via the graviton loop. Their coefficients are
proportional to the product of the non–minimal couplings and the loop factor,
and thus expected to be significant (in the absence of symmetry arguments). In
general, it is a challenging task to estimate the relative size of the
different operators since this can only be done reliably within UV complete
gravity theories.

On shell, two of the derivative operators can be eliminated via integration by
parts:

$\displaystyle(\partial_{\mu}\phi)^{2}s^{2}\rightarrow(\partial_{\mu}s)^{2}\phi^{2}+m_{\phi}^{2}\phi^{2}s^{2}\;,$
(4)
$\displaystyle(\phi\partial_{\mu}\phi)(s\partial^{\mu}s)\rightarrow-{1\over
2}(\partial_{\mu}s)^{2}\phi^{2}\;,$ (5)

where we have neglected the dark matter mass and the Hubble rate, which,
during preheating is small compared to the particle energy. (We consider the
Hubble–induced effects in Section 4). Focussing on dark matter pair
production, we can thus restrict ourselves to the operators

${\cal O}_{3}={1\over M_{\rm
Pl}^{2}}\;(\partial_{\mu}s)^{2}\phi^{2}~{}~{},~{}~{}{\cal O}_{4}={1\over
M_{\rm Pl}^{2}}\;\phi^{4}s^{2}\;,$ (6)

amended with the renormalizable interaction

${\cal O}_{\rm renorm}={m_{\phi}^{2}\over M_{\rm Pl}^{2}}\;\phi^{2}s^{2}\;.$
(7)

Although this term is renormalizable, the coupling strength is highly
suppressed: for typical inflaton masses it is below $10^{-10}$. The operator
$\phi^{2}s^{4}$ produces a final state with 4 DM quanta. The corresponding
reaction rate is similar to that of the derivative operator, since the final
state phase space gives analogous energy dependence. Therefore, we will not
discuss this operator separately within the perturbative regime, while its
non–perturbative analysis will be presented in Sec. 4.

Clearly, operators ${\cal O}_{3}$ and ${\cal O}_{4}$ exhibit qualitatively
different behaviour in regard to dark matter production. Indeed, ${\cal
O}_{3}$ involves the particle energy which is of the order of the (effective)
inflaton mass, while in ${\cal O}_{4}$ this dependence is replaced by the
inflaton field value. The latter is not far from the Planck scale in typical
models, thus

$\phi\gg E_{\phi}$ (8)

and one expects much more efficient DM production from ${\cal O}_{4}$. In what
follows, we make this argument more quantitative.

##  3 Perturbative dark matter production

An oscillating classical background can lead to particle production [11, 12,
13]. After inflation, $\phi$ oscillates coherently in either $\phi^{2}$ or
$\phi^{4}$ potential, depending on the inflationary model. As a result, the
$\phi-s$ couplings induce dark matter pair production. If the corresponding
coupling is small, the process can be described perturbatively. Below, we
consider the contributions of the 3 basic operators to this reaction. We treat
the Hubble expansion adiabatically such that the time dependence can be
inserted in the inflaton oscillation amplitude at the end of the calculation.
Also, we treat the produced dark matter particles as free and neglect
backreaction. These approximations are justified at small inflaton–DM
couplings.

###  3.1 $\phi^{2}s^{2}$ interaction

Consider the 4–point interaction

$-\Delta{\cal L}_{\rm renorm}={1\over 4}\lambda_{\phi s}\,\phi^{2}s^{2}\;,$
(9)

where $\lambda_{\phi s}\sim m_{\phi}^{2}/M_{\rm Pl}^{2}$. Let us expand the
inflaton field as

$\phi^{2}(t)=\sum_{n=-\infty}^{\infty}\zeta_{n}e^{-in\omega t}\;,$ (10)

where the coefficients $\zeta_{n}$ are time–independent. Creation of a
two–particle DM state with momenta $p,q$ from the vacuum is described by the
amplitude (in Peskin–Schroeder conventions [15])

$-i\int_{-\infty}^{\infty}dt\langle f|V(t)|i\rangle=-i\,{\lambda_{\phi s}\over
2}\,(2\pi)^{4}\delta({\bf{p}}+{\bf{q}})\sum_{n=1}^{\infty}\zeta_{n}\delta(E_{p}+E_{q}-n\omega)\;,$
(11)

with $V(t)$ given by Eq. 9. The corresponding invariant amplitude for the
$n$-th inflaton mode decay is ${\cal M}_{n}=-\lambda_{\phi s}\zeta_{n}/2$. The
resulting reaction rate for DM pair production per unit volume is

$\Gamma=\sum_{n=1}^{\infty}\Gamma_{n}=\sum_{n=1}^{\infty}{1\over 2}\int|{\cal
M}_{n}|^{2}d\Pi_{n}={\lambda_{\phi s}^{2}\over
64\pi}\sum_{n=1}^{\infty}|\zeta_{n}|^{2}\sqrt{1-\left({2m_{s}\over
n\omega}\right)^{2}}\;\theta(n\omega-2m_{s})\;.$ (12)

Here we have kept the DM mass for generality. The inflaton decay rate
$\Gamma_{\phi}$ can be computed using energy conservation,
$\rho_{\phi}\Gamma_{\phi}=\langle E\rangle\Gamma$, where $\rho_{\phi}$ and
$\langle E\rangle$ are the inflaton energy density and the average energy of
the decay products, respectively. Hence,

$\Gamma_{\phi}={\lambda_{\phi s}^{2}\omega\over
64\pi\rho_{\phi}}\;\sum_{n=1}^{\infty}n|\zeta_{n}|^{2}\sqrt{1-\left({2m_{s}\over
n\omega}\right)^{2}}\;\theta(n\omega-2m_{s})\;.$ (13)

In the massless limit, one thus recovers the result of [16].

###  3.2 $\phi^{2}(\partial_{\mu}s)^{2}$ interaction

The calculation proceeds as above, except the final state receives an
additional momentum–dependent factor in the amplitude:

$p\cdot q={1\over 2}(p+q)^{2}={1\over 2}(E_{p}+E_{q})^{2}={1\over
2}n^{2}\omega^{2}\;,$

where the DM mass has been neglected. So, effectively in the amplitude for the
quartic interaction, one replaces $\lambda_{\phi s}\rightarrow 4C_{3}/M_{\rm
Pl}^{2}\,p\cdot q$, which in the $m_{s}\rightarrow 0$ limit leads to

$\Gamma={C_{3}^{2}\,\omega^{4}\over 16\pi M_{\rm
Pl}^{4}}\,\sum_{n=1}^{\infty}n^{4}|\zeta_{n}|^{2}\;.$ (14)

###  3.3 $\phi^{4}s^{2}$ interaction

The calculation is similar to that for the $\phi^{2}s^{2}$ case. $\phi^{4}$
can be expanded as

$\phi^{4}(t)=\sum_{n=-\infty}^{\infty}\hat{\zeta}_{n}e^{-in\omega t}\;,$ (15)

where

$\hat{\zeta}_{n}=\sum_{m=-\infty}^{\infty}\zeta_{n-m}\zeta_{m}\;.$ (16)

Replacing $\lambda_{\phi s}\rightarrow 4C_{4}/M_{\rm Pl}^{2}$, in the massless
DM limit we get

$\Gamma={C_{4}^{2}\over 4\pi M_{\rm
Pl}^{4}}\,\sum_{n=1}^{\infty}|\hat{\zeta}_{n}|^{2}\;.$ (17)

###  3.4 Relative efficiency

Let us estimate the relative particle production efficiency of the different
operators. If the inflaton potential is quadratic, $V(\phi)={1\over
2}m_{\phi}^{2}\phi^{2}$, we have

$\phi(t)=\phi_{0}\cos m_{\phi}t\;,$ (18)

where $\phi_{0}$ is the oscillation amplitude. (The fact that $\phi_{0}$
decreases slowly in time, $\phi_{0}\propto 1/(m_{\phi}t)$, is insignificant
for our purposes). In the quartic case, $V(\phi)={1\over
4}\lambda_{\phi}\phi^{4}$, the inflaton field is given by the Jacobi cosine,

$\phi(t)=\phi_{0}\,{\rm
cn}\left(\sqrt{\lambda_{\phi}}\phi_{0}\,t,{1\over\sqrt{2}}\right)={\sqrt{\pi}\Gamma\left({3\over
4}\right)\over\Gamma\left({5\over
4}\right)}\,\phi_{0}\,\sum_{n=1}^{\infty}\left(e^{i(2n-1)\omega
t}+e^{-i(2n-1)\omega t}\right)\,{e^{-(\pi/2)(2n-1)}\over 1+e^{-\pi(2n-1)}}\;,$
(19)

where

$\omega={1\over 2}\sqrt{\pi\over 6}{\Gamma\left({3\over
4}\right)\over\Gamma\left({5\over 4}\right)}\,m_{\phi}^{\rm eff}$ (20)

and

$m_{\phi}^{\rm eff}=\sqrt{3\lambda_{\phi}}\phi_{0}\;.$ (21)

For many purposes, the above sum can be approximated by the first term with
$n=1$.

The relative efficiency of ${\cal O}_{3}$ and ${\cal O}_{4}$ is given by

${\Gamma\left[{\cal O}_{3}\right]\over\Gamma\left[{\cal
O}_{4}\right]}\sim{C_{3}^{2}\over
C_{4}^{2}}\,{\omega^{4}\over\phi_{0}^{4}}\;.$ (22)

Clearly, the reaction rate due to $\phi^{2}(\partial_{\mu}s)^{2}$ is much
suppressed compared to that of $\phi^{4}s^{2}$. For the quadratic inflaton
potential, the suppression factor is

$m_{\phi}^{4}/\phi_{0}^{4}\sim 10^{-20}\;,$ (23)

assuming the typical values $\phi_{0}\sim 1$ and $m_{\phi}\sim 10^{-5}$ in
Planck units. In the quartic case,
$\omega^{4}/\phi_{0}^{4}\sim\lambda_{\phi}^{2}<10^{-20}$ for typical
$\lambda_{\phi}<10^{-10}$.

The contribution of the $\phi^{2}s^{2}$–operator of the form (7) is similarly
suppressed by $m_{\phi}^{4}/\phi_{0}^{4}$. It is also clear that, due to the
phase space integral, the rate of the ${\cal O}_{5}$–induced process
$\phi\phi\rightarrow ssss$ contains an additional factor of
$E_{\phi}^{4}/\phi_{0}^{4}$ compared to the pair production rate from ${\cal
O}_{4}$. Hence we conclude that the $\phi^{4}s^{2}$ interaction dominates DM
production, unless there is a large hierarchy in the Wilson coefficients, e.g.
$C_{3}/C_{4}\sim 10^{10}$.

The above calculation also tells us that higher dimensional operators

${\tilde{C}_{6}\over M_{\rm Pl}^{4}}\phi^{6}s^{2}+{\tilde{C}_{8}\over M_{\rm
Pl}^{6}}\phi^{8}s^{2}+...$ (24)

are important. Indeed, their contributions are only suppressed by
$\phi_{0}^{4}/M_{\rm Pl}^{4}$, etc. relative to that of ${\cal O}_{4}$. If the
inflaton amplitude is not far away from the Planck scale, this suppression is
not very significant.

##  4 Resonant dark matter production via dim–6 operators

Perturbative calculations ignore the Bose enhancement of the amplitude due to
the presence of identical states. Depending on the coupling, this enhancement
can be very significant and lead to resonant production [12, 17, 18]. Such a
regime can be described semiclassically by analyzing the equations of motion
for the DM field $s$. In what follows, we compare the corresponding resonant
particle production via operators ${\cal O}_{3}$ and ${\cal O}_{4}$.

To derive the EOM for $s$ in the presence of the inflaton background
$\phi(t)$, consider the action

$S=\int d^{4}x\sqrt{|g|}\left({1\over 2}\,{\cal
K}(\phi)\,g^{\mu\nu}\partial_{\mu}s\,\partial_{\nu}s-{\cal V}\right)\;,$ (25)

where the kinetic function ${\cal K}(\phi)$ depends on $\phi$ only and ${\cal
V}$ is the $s-dependent$ part of the scalar potential. Here $g={\rm
det}\,g_{\mu\nu}$ and the Friedmann metric is $g_{\mu\nu}={\rm
diag}(1,-a^{2},-a^{2},-a^{2})$. Since $\phi$ and $a$ are functions of time
only, variation of the action with respect to $s$ yields

$\ddot{s}-{1\over a^{2}}\,\partial_{i}\partial_{i}s+\left({\dot{\cal
K}\over{\cal K}}+3H\right)\dot{s}+{{\cal V}^{\prime}_{s}\over{\cal K}}=0\;.$
(26)

Let us expand $s(t,{\bf x})$ in spacial Fourier modes $s_{k}(t)$, where $k$ is
the comoving 3–momentum. If $V$ is quadratic in $s$, the different momentum
modes decouple and we have

$\ddot{s}_{k}+\left({\dot{\cal K}\over{\cal
K}}+3H\right)\dot{s}_{k}+\left({k^{2}\over a^{2}}+{{\cal
V}^{\prime\prime}_{s}\over{\cal K}}\right)s_{k}=0\;.$ (27)

Now suppose that the renormalizable potential vanishes and only one dim–6
operator is present at a time, such that either ${\cal K}=1$ or ${\cal V}=0$.

(a) Consider first the case of a non–trivial kinetic function ${\cal K}$ and
${\cal V}=0$. The first order time derivatives can be eliminated by the change
of variables

$s_{k}=a^{-3/2}{\cal K}^{-1/2}X_{k}\;,$ (28)

such that

$\ddot{X}_{k}+\left[{1\over 4}{\dot{\cal K}^{2}\over{\cal K}^{2}}-{1\over
2}{\ddot{\cal K}\over{\cal K}}-{3\over 2}H{\dot{\cal K}\over{\cal K}}+{9\over
4}wH^{2}+{k^{2}\over a^{2}}\right]X_{k}=0\;,$ (29)

where the equation of state coefficient is

$w=-\left(1+{2\dot{H}\over 3H^{2}}\right)\;.$ (30)

In the limit of a constant inflaton amplitude, the coefficients are periodic
in time such that the above EOM belongs to the class of Hill’s equations.
Depending on the parameters, the solution $X_{k}$ can grow exponentially in
time signifying particle production.

Let us now specialize to operator ${\cal O}_{3}$,

${\cal K}=1+{2C_{3}\over M_{\rm Pl}^{2}}\,\phi^{2}\;.$ (31)

The effective field theory expansion makes sense if $C_{3}\phi^{2}/M_{\rm
Pl}^{2}\ll 1$. Consider the $locally$ quadratic inflaton potential such that
$\phi(t)=\phi_{0}(t)\cos m_{\phi}t$ with $\phi_{0}\propto 1/(m_{\phi}t)$. In
this case, the Universe is matter dominated, $w=0$ and

$H={m_{\phi}\phi_{0}\over\sqrt{6}M_{\rm Pl}}\;.$ (32)

For $C_{3}<1$ and $\phi_{0}<M_{\rm Pl}$, the term $(\dot{\cal K}/{\cal
K})^{2}$ in the square brackets is insignificant, as is $wH^{2}$. The terms
$H\dot{\cal K}/{\cal K}$ and $\ddot{\cal K}/{\cal K}$ are similar in magnitude
initially, but the former is cubic in $\phi_{0}$ and thus decreases faster in
time. Keeping just the $\ddot{\cal K}/{\cal K}$ term, we may approximate

$\ddot{X}_{k}+\left[{2C_{3}m_{\phi}^{2}\phi_{0}^{2}\over M_{\rm Pl}^{2}}\,\cos
2m_{\phi}t+{k^{2}\over a^{2}}\right]X_{k}=0\;.$ (33)

This has the form of the Mathieu equation,

$X_{k}^{\prime\prime}+\left[A+2q\,\cos 4z\right]\,X_{k}=0\;,$ (34)

where $z=m_{\phi}t/2$, the prime stands for differentiation with respect to
$z$, and

$q={4C_{3}\phi_{0}^{2}\over M_{\rm Pl}^{2}}~{}~{},~{}~{}A={4k^{2}\over
a^{2}m_{\phi}^{2}}\;.$ (35)

A large $q$ generally implies fast amplitude growth and efficient particle
production via broad parametric resonance [17]. In our case, however, $q\ll 1$
and the resonance is narrow. As a result, no efficient particle production is
possible, especially in view of the redshifting of the produced particle
momenta [17],[19].

(b) Let us now consider the effect of ${\cal O}_{4}$, so we take ${\cal K}=1$
and

${\cal V}={C_{4}\over M_{\rm Pl}^{2}}\;\phi^{4}s^{2}\;.$ (36)

For the quadratic inflaton potential, the corresponding $X_{k}$ satisfies

$\ddot{X}_{k}+\left[{2C_{4}\phi_{0}^{4}\over M_{\rm
Pl}^{2}}\,\cos^{4}m_{\phi}t+{k^{2}\over a^{2}}\right]X_{k}=0\;.$ (37)

Using $\cos^{4}x={1\over 8}\cos 4x+{1\over 2}\cos 2x+{3\over 8}$, one can
bring the EOM into the form of the Whittaker–Hill equation,

$X_{k}^{\prime\prime}+\left[A+2p\cos 2z+2q\cos 4z\right]X_{k}=0\;,$ (38)

where now $z=m_{\phi}t$ and

$q={C_{4}\phi_{0}^{4}\over 8m_{\phi}^{2}M_{\rm
Pl}^{2}}~{}~{},~{}~{}p={C_{4}\phi_{0}^{4}\over 2m_{\phi}^{2}M_{\rm
Pl}^{2}}~{}~{},~{}~{}A={k^{2}\over a^{2}m_{\phi}^{2}}+{3C_{4}\phi_{0}^{4}\over
4m_{\phi}^{2}M_{\rm Pl}^{2}}\;.$ (39)

The efficiency of particle production is characterized by $p$ and $q$, whose
large values (depending on $A$) generally lead to broad resonance [20],[21].
We see that this regime is easily achieved in the presence of ${\cal O}_{4}$.
In particular, $p,q\gg 1$ is consistent with $C_{4}\ll 1$ and sub–Planckian
$\phi_{0}$ values, as long as $\phi_{0}\gg m_{\phi}$.

We therefore conclude that ${\cal O}_{4}$ is much more efficient in particle
production than ${\cal O}_{3}$. The same conclusion applies to the quartic
inflaton potential: the analysis proceeds analogously up to the replacement of
the inflaton mass with the effective inflaton mass $m_{\phi}^{\rm
eff}\sim\sqrt{\lambda_{\phi}}\phi_{0}$.

The resonance efficiency is determined by the ratio of the inflaton–induced DM
mass to the inflaton effective mass. In the list (3), this ratio can be large
only for operator ${\cal O}_{4}$. Indeed, $(\partial_{\mu}\phi)^{2}s^{2}$
induces the DM mass of order $m_{\phi}$ at best. The term
$(\phi\partial_{\mu}\phi)(s\partial^{\mu}s)=(\partial_{\mu}\phi^{2})(\partial^{\mu}s^{2})/4$
can be rewritten as a mass term for $s$ by integrating by parts. The resulting
DM mass scale is determined by $m_{\phi}$ or $H$, leading to the same
conclusion. The operator $\phi^{2}s^{4}$ does not induce any mass in our
approximation and the corresponding EOM is not of Hill’s type (see
below).222The induced mass term appears when the variance $\langle
s^{2}\rangle$ becomes significant. Finally, the renormalizable term
$\phi^{2}s^{2}$ of the form (7) does not lead to broad resonance since the
corresponding $q\lesssim 1$. It is thus clear that ${\cal O}_{4}$ dominates
particle production.

Similar conclusions apply to higher dimensional operators $\phi^{6}s^{2}$,
etc. As long as $\phi_{0}$ is not much below the Planck scale, the effective
$q$ parameter can be much greater than one, signifying efficient particle
production. Therefore, the results are sensitive to the presence of operators
of this type.

The amount of dark matter produced during preheating is difficult to estimate
analytically. The reason is that the parameters of the Whittaker–Hill equation
evolve in time making the resonance stochastic, which is further complicated
by the non–trivial 3-D stability band structure [21]. Depending on the size of
$C_{4}$, tangible backreaction and rescattering effects [22],[23] can also
take place. Thus, to make reliable predictions, we have to resort to lattice
simulations.

###  4.1 On resonant production via $\phi^{2}s^{4}$

Unlike for other operators considered in this work, resonant particle
production via ${\cal O}_{5}=\phi^{2}s^{4}$ is not described by the Hill’s
equation. In this case, the system is non–linear already to leading order and
thus difficult to handle analytically. In this subsection, we discuss some of
its properties relevant to the subject of our paper.

The EOM for the DM field $s$ in the presence of ${\cal O}_{5}$ reads

$\ddot{s}-{1\over
a^{2}}\,\partial_{i}\partial_{i}s+3H\,\dot{s}+{4C_{5}\phi^{2}(t)\over M_{\rm
Pl}^{2}}\,s^{3}=0\;.$ (40)

Clearly, the different momentum modes do not decouple in this case. As a
representative example, let us focus on the zero mode of $s$ which normally
plays a major role in particle production. Omitting the gradient term and
introducing a rescaled field $X$ ($s=a^{-3/2}X$), we get

$\ddot{X}+{9\over 4}wH^{2}\,X+{4C_{5}\phi^{2}(t)\over M_{\rm
Pl}^{2}\,a^{3}}\,X^{3}=0\;,$ (41)

where $w$ is the coefficient of the equation of state of the system. Let us
now specialize to the quadratic inflaton potential such that $w=0$,

${a\over a_{0}}=\left({m_{\phi}t\over m_{\phi}t_{0}}\right)^{2/3}\;,$ (42)

and $\phi(t)\simeq 1.85\,\varphi_{0}{\cos\,m_{\phi}t\over m_{\phi}t}$. Here
the initial condition is chosen such that $\phi(t_{0})=\varphi_{0}$ with
$m_{\phi}t_{0}=1$. The above EOM should be supplemented by boundary
conditions. The magnitude of a light field is given by quantum fluctuations,
such that we may take $s\sim H$ and $\dot{s}\sim H^{2}$ initially, where
$H=m_{\phi}\varphi_{0}/(\sqrt{6}M_{\rm Pl})$. Introducing

$z=m_{\phi}t~{}~{},~{}~{}Y=X/m_{\phi}\;,$ (43)

we get

$Y^{\prime\prime}+\kappa\,{\cos^{2}z\over z^{4}}\,Y^{3}=0\;,$ (44)

where the prime denotes differentiation with respect to $z$. The
representative boundary conditions can be chosen as $Y(1)=1,Y^{\prime}(1)=1$.
The coupling $\kappa$ is given by

$\kappa\simeq{14C_{5}\varphi_{0}^{2}\over M_{\rm Pl}^{2}\,a_{0}^{3}}\;,$ (45)

where typically $\varphi_{0}\sim M_{\rm Pl}$, $a_{0}\sim 1$ and the effective
field theory approach is expected to be valid for $\kappa\lesssim{\cal O}(1)$.

Figure 1:  Solutions to Eq. 44 for $\kappa=1,100,2010$.

The oscillating term in Eq. 44 falls faster with time then the analogous
coefficient in the Mathieu equation, hence the duration of the resonance is
shorter in the present case. The equation exhibits a simple asymptotic
solution for $z\gg 1$ (or $\kappa\ll 1$),

$Y\propto z\;.$ (46)

On the other hand, for $\kappa\gg 1$, $Y(z)$ varies much faster than
$\cos^{2}z/z^{4}$ does such that the latter can be treated adiabatically. In
this case, Eq. 44 takes the form $Y^{\prime\prime}+c\,Y^{3}=0$ with
$c=\kappa\,\cos^{2}z/z^{4}$, whose solution is a Jacobi cosine. Given the
amplitude of oscillations $Y_{0}$, locally we have

$Y(z)\simeq Y_{0}\,{\rm cn}\left(\sqrt{c}\,Y_{0}z,{1\over\sqrt{2}}\right)$
(47)

in the convention of [18]. The oscillation frequency $\sqrt{c}Y_{0}$ changes
non–adiabatically around the inflaton zero crossings, implying particle
production.

Numerical solutions to Eq. 44 are presented in Fig. 1. For $\kappa=1$, the
solution quickly takes on the asymptotic form (46), which corresponds to
constant $s$ in our approximation ($m_{s}\sim 0$). Even for $\kappa=100$,
there is no significant amplitude growth, while for $\kappa=2010$ the solution
exhibits truly resonant behaviour.333We find that the solution in this regime
is quite sensitive to the numerical method used for solving the differential
equation.  We thus find that $\kappa\gtrsim 10^{3}$ is necessary for efficient
DM production.

The discussion above concerns the zero mode. Normally it serves as an
indicator whether or not resonant particle production takes place. To reaffirm
it, we have performed lattice simulations of the full inflaton–DM system in
the regime $C_{5}\phi^{2}/M_{\rm Pl}^{2}\lesssim 1$ and indeed found no
tangible particle production.

We therefore conclude that no efficient particle production is induced by
operator ${\cal O}_{5}$ as long as $C_{5}\phi^{2}/M_{\rm Pl}^{2}\lesssim{\cal
O}(1)$.

##  5 Lattice simulations: reproducing the correct DM abundance

In this Section, we focus on the leading operator ${C_{4}\over M_{\rm
Pl}^{2}}\,\phi^{4}s^{2}$ and compute the resulting DM relic abundance in the
resonant regime, $C_{4}\phi_{0}^{4}/M_{\rm Pl}^{2}\gg m_{\phi}^{2}$, assuming
$m_{\phi}\gg m_{s}$. As explained above, the analytical approach to resonant
particle production has significant limitations, given the complexity of the
Whittaker–Hill equation as well as backreaction and rescattering effects.
Therefore, we resort to lattice simulations using the numerical tool
CosmoLattice [24, 25].

A realistic framework must also account for the Standard Model particle
production. The simplest way to incorporate reheating is to include a small
Higgs–inflaton coupling following [26, 27],

$V_{\phi h}=\sigma_{\phi h}\phi H^{\dagger}H\;,$ (48)

which would lead to late–time decay of the inflaton into the Higgs pairs (for
$m_{\phi}>2m_{h}$). As long as $\sigma_{\phi h}$ is sufficiently small in
Planck units, resonant dark matter production is not affected by this
coupling. The dark matter abundance is expressed in terms of

$Y={n\over s_{\rm SM}}~{}~{},~{}~{}s_{\rm SM}={2\pi^{2}\over
45}\,g_{*s}\,T^{3}\;,$ (49)

where $n$ is the DM number density, $s_{\rm SM}$ is the Standard Model entropy
density at temperature $T$ and $g_{*s}$ is the effective number of SM degrees
of freedom contributing to the entropy. We are interested in very weakly
interacting dark matter such that it never reaches thermal equilibrium with
the environment. After preheating ends, the total number of the DM quanta
remains constant. Since the SM entropy is also conserved, $Y$ can be computed
at the reheating stage. The observed value is [28]

$Y_{\infty}=4.4\times 10^{-10}\;\left({{\rm GeV}\over m_{s}}\right)\;,$ (50)

which sets a constraint on the model parameters.

Reheating occurs when

$H_{R}\simeq\Gamma_{\phi\rightarrow hh}\;,\quad\Gamma_{\phi\rightarrow
hh}=\frac{\sigma_{\phi h}^{2}}{8\pi m_{\phi}},$ (51)

where $H_{R}$ is the Hubble rate at reheating and $\Gamma_{\phi\rightarrow
hh}$ takes into account 4 Higgs d.o.f. at high energies. The reheating
temperature is given by

$H_{R}=\sqrt{\frac{\pi^{2}g_{*}}{90}}\,{T_{R}^{2}\over M_{\rm Pl}}\;,$ (52)

where $g_{*}$ is the effective number of the Standard Model degrees of freedom
contributing to the energy density. Combining $T_{R}$ with the dark matter
density $n$ computed on the lattice, one determines $Y$ according to (49).

The resulting abundance is sensitive to the energy balance between the
inflaton and dark matter, which affects the expansion history. Depending on
the coupling and the inflaton potential, DM can contribute a significant
fraction up to 50% to the total energy density at the end of preheating. We
therefore parametrize

$\rho_{e}(s)=\delta\,\rho_{e}(\phi)\;,$ (53)

where $\rho_{e}(s)$, $\rho_{e}(\phi)$ are the DM and inflaton energy
densities, respectively, evaluated at the end of the simulation. At weak
coupling, in the quadratic inflaton potential we have $\delta\sim 0$, while,
at strong coupling, $\delta$ can reach a value close to 1. The Universe
evolution proceeds in stages: first, it can be dominated by radiation; later,
when the energy per quantum becomes comparable to the inflaton mass, it
evolves as non–relativistic matter; finally, the Universe reheats and becomes
radiation–like. Denoting the corresponding scale factors as $a_{e}$ (end of
the simulation), $a_{*}$ (transition), $a_{R}$ (reheating), we have

$a_{e}\stackrel{{\scriptstyle\rm{rel}}}{{\longrightarrow}}a_{*}\stackrel{{\scriptstyle\rm
nrel}}{{\longrightarrow}}a_{R}\;,$ (54)

such that the Hubble rate evolves as $H\sim a^{-3(w+1)/2}$ with $w=1/3$ and
$w=0$ during the two periods, respectively. After the transition point
$a_{*}$, $\rho(s)$ becomes negligible and at $a_{R}$ the inflaton energy
density $\rho(\phi)$ converts into SM radiation. Thus,

$H_{R}=\frac{H_{e}}{\sqrt{1+\delta}}\,\frac{a_{e}^{2}}{a_{*}^{2}}\,\frac{a_{*}^{3/2}}{a_{R}^{3/2}}\;,$
(55)

where $H_{e}$ is the Hubble rate at the end of the simulation. We note that
the first stage of the radiation–like expansion may collapse to a point, i.e.
$a_{e}=a_{*}$. This is the case for the quadratic inflaton potential at weak
inflaton–DM coupling.

Solving for $a_{R}$, we find $\sigma_{\phi h}$ required by the correct DM
abundance in terms of the simulation output:

$\sigma_{\phi h}\simeq 1.6\times 10^{-8}\;\sqrt{m_{\phi}\,M_{\rm
Pl}^{3}}\;{H_{e}^{2}\over(1+\delta)\;n_{e}}\,{a_{e}\over
a_{*}}~{}~{}\left({{\rm GeV}\over m_{s}}\right)$ (56)

for $g_{*}\simeq 107$ and $M_{\rm Pl}$ being the reduced Planck mass. The
values of $H_{e},n_{e},\delta$ at the end of preheating are computed by
CosmoLattice, while $a_{e}/a_{*}$ can be determined by tracking the equation
of state of the system.

This formula can be simplified further if we define $a_{*}$ according to

$\langle E_{e}(\phi)\rangle\;{a_{e}\over a_{*}}\simeq m_{\phi}\;,$ (57)

where the average energy of the inflaton quantum at the end of the simulation
is $\langle E_{e}(\phi)\rangle=\rho_{e}(\phi)/n_{e}(\phi)$. This definition of
$a_{*}$ is more practical in that it does not require tracking the equation of
state of the system over a long period, which is computationally challenging.
We then get

$\sigma_{\phi h}\simeq 5\times 10^{-9}\;{m_{\phi}^{3/2}\over M_{\rm
Pl}^{1/2}}~{}{n_{e}(\phi)\over n_{e}(s)}~{}~{}\left({{\rm GeV}\over
m_{s}}\right)\;,$ (58)

which only requires the particle densities as an output of the simulations.
Here $n_{e}(\phi)$ includes the inflaton quanta with zero momentum and
typically $n_{e}(\phi)\gg n_{e}(s)$ unless the coupling is relatively strong.

The number densities are computed via the $k$–mode occupation numbers $n_{k}$.
For the dark matter field, we have

$n_{k}={\omega_{k}\over
2}\,\left({|\dot{X}_{k}|^{2}\over\omega_{k}^{2}}+|X_{k}|^{2}\right)-{1\over
2}\;,$ (59)

where $\omega_{k}^{2}(t)={2C_{4}\phi_{0}^{4}\over M_{\rm
Pl}^{2}}\,\cos^{4}m_{\phi}t+{k^{2}\over a^{2}}$ in the quadratic inflaton
potential. Here $X_{k}$ is a solution to the EOM with the boundary condition
given by quantum fluctuations. The resulting number density is then given by

$n(s)={1\over(2\pi a)^{3}}\int d^{3}k\;{n_{k}}\,.$ (60)

On the lattice, the momentum spectrum is discrete which allows one to treat
the zero mode separately. Analogous formulae apply to the inflaton field and
the quartic potential. It is important to remember that the EOM for the
different momentum modes decouple at weak couplings only. The lattice approach
allows us to incorporate the couplings among the $k$–modes of the inflaton and
DM, thereby accounting for backreaction and rescattering. The latter can have
a crucial impact on the dynamics of the system (see [29, 30] for recent
examples).

Figure 2:  $C_{4}$ vs $\sigma_{\phi h}$ required for the correct DM abundance
in the quadratic inflaton potential ($m_{\phi}=10^{13}$ GeV, $m_{s}=1$ GeV,
$\varphi_{0}\simeq M_{\rm Pl}$). The simulations are performed with
CosmoLattice [24].

Our numerical results are presented in Fig. 2. The correct relic density is
produced in a wide range of $C_{4}$ between $10^{-7}$ and $10^{-4}$. The lower
bound comes from requiring semiclassical behaviour, that is, the occupation
numbers must be sufficiently large. The upper bound has technical nature: the
simulation becomes unstable. We observe that the curve tends to flatten out at
larger couplings. This is expected from quasi–equilibrium: as $n_{e}(\phi)$
approaches $n_{e}(s)$, the $\sigma_{\phi h}$ coupling becomes constant within
our approximation [27]. Although such flattening is clearly visible, we find
that quasi–equilibrium has not yet been reached at $C_{4}\sim 10^{-4}$. We
estimate the required $C_{4}$ to be of order $10^{-3}$, yet the simulations in
this range become less reliable.

Given these results, we can now determine under what circumstances the
renormalizable coupling ${1\over 4}\lambda_{\phi s}\phi^{2}s^{2}$ becomes
unimportant. According to Fig. 4 of [27], $\lambda_{\phi s}<10^{-8}$ does not
make any significant contribution to the dark matter abundance in the
parameter range of interest. This can be understood intuitively since
$\lambda_{\phi s}/4<C_{4}\varphi^{2}_{0}/M_{\rm Pl}^{2}$ in this case (see
also [20]). On the other hand, for $\lambda_{\phi s}>10^{-7}-10^{-6}$ or
$C_{4}>10^{-3}$, the inflaton–DM system reaches quasi–equilibrium and the DM
abundance becomes independent of these couplings, with the required
$\sigma_{\phi h}\sim 10^{-17}M_{\rm Pl}$.

Subsequently, when the inflaton coherence is lost, the operator ${\cal O}_{4}$
relinquishes its privileged role in DM production. In quasi–equilibrium, the
scattering processes $\phi\phi\rightarrow ss$ can become comparably
significant, depending on the corresponding $C_{i}$. Since all of our
operators are Planck–suppressed, such processes are slower than the Hubble
rate and thus do not lead to inflaton thermalization [31] nor significant DM
production. At weak couplings, the inflaton field may remain semi–classical
during reheating, in which case the effects of the SM thermal bath can be
included along the lines of [32],[33].

Throughout this work we assume that other sources of dark matter are
subdominant. In particular, we neglect the DM coupling to the Higgs such that
no tangible freeze–in contribution from the Higgs thermal bath appears. This
approximation is justified if this coupling is below $10^{-11}$ [34].
Furthermore, as explained in the Introduction, there is also a truly
gravitational source of DM: the Universe expansion. However, in the presence
of ${\cal O}_{4}$, the effective mass of $s$ during inflation can be larger
than the Hubble rate due to super–Planckian inflaton values, which suppresses
$\langle s^{2}\rangle$ and makes this production mechanism inefficient.

In our analysis, we have relied on the effective field theory expansion in the
Einstein frame, which is expected to be meaningful during preheating.
Gravitational dark matter production can also be encoded in the DM non–minimal
coupling to gravity [35], which corresponds to a specific choice of higher
dimensional operators in our approach. A related option is provided by
gravity–induced inflaton decay in Starobinsky–like models [36].

##  6 Conclusion

In this work, we have studied perturbative and non–perturbative dark matter
production during the inflaton oscillation epoch. We focus on the regime where
the renormalizable interactions between the inflaton and dark matter are
negligible. To determine the leading contributions, we resort to the effective
field theory expansion in the inverse Planck mass. Such higher dimensional
operators are expected to be generated by perturbative or non–perturbative
gravitational effects. In the absence of quantum gravity theory, their
coefficients are unknown and therefore treated as arbitrary input parameters.

We have focussed on Planck–suppressed dim–6 operators and studied their
relative importance in the perturbative and resonant regimes. We find that
operators of the form $\phi^{n}s^{2}$ ($n\geq 4$) by far dominate particle
production. They can generate the correct (non–thermal) dark matter abundance
even for small values of the corresponding Wilson coefficients. Therefore, the
phenomenological frameworks describing dark matter production are sensitive to
the presence of such operators, which reinforces the importance of
gravitational effects in this context.

Acknowledgements. The authors wish to thank the Finnish Computing Competence
Infrastructure (FCCI) for supporting this project with computational and data
storage resources.

## References

  * [1] L. Parker, Phys. Rev. 183, 1057-1068 (1969); A. A. Grib and S. G. Mamaev, Yad. Fiz. 10, 1276-1281 (1969).
  * [2] S. G. Mamaev, V. M. Mostepanenko and A. A. Starobinsky, Zh. Eksp. Teor. Fiz. 70, 1577-1591 (1976); A. A. Grib, S. G. Mamaev and V. M. Mostepanenko, Gen. Rel. Grav. 7, 535-547 (1976).
  * [3] V. Kuzmin and I. Tkachev, Phys. Rev. D 59, 123006 (1999).
  * [4] N. A. Chernikov and E. A. Tagirov, Ann. Inst. H. Poincare Phys. Theor. A 9, 109 (1968).
  * [5] D. J. H. Chung, P. Crotty, E. W. Kolb and A. Riotto, Phys. Rev. D 64, 043503 (2001).
  * [6] A. A. Starobinsky and J. Yokoyama, Phys. Rev. D 50, 6357-6368 (1994).
  * [7] K. Enqvist, S. Nurmi, T. Tenkanen and K. Tuominen, JCAP 08, 035 (2014).
  * [8] T. Markkanen, A. Rajantie and T. Tenkanen, Phys. Rev. D 98, no.12, 123532 (2018).
  * [9] Y. Mambrini and K. A. Olive, Phys. Rev. D 103, no.11, 115009 (2021).
  * [10] S. Clery, Y. Mambrini, K. A. Olive and S. Verner, [arXiv:2112.15214 [hep-ph]].
  * [11] A. D. Dolgov and D. P. Kirilova, Sov. J. Nucl. Phys. 51, 172-177 (1990).
  * [12] L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. Lett. 73, 3195-3198 (1994).
  * [13] Y. Shtanov, J. H. Traschen and R. H. Brandenberger, Phys. Rev. D 51, 5438-5455 (1995).
  * [14] D. S. Salopek, J. R. Bond and J. M. Bardeen, Phys. Rev. D 40, 1753 (1989).
  * [15] M. E. Peskin and D. V. Schroeder, ‘‘An Introduction to quantum field theory,’’ Addison-Wesley, 1995.
  * [16] K. Ichikawa, T. Suyama, T. Takahashi and M. Yamaguchi, Phys. Rev. D 78, 063545 (2008).
  * [17] L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56, 3258-3295 (1997).
  * [18] P. B. Greene, L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56, 6175-6192 (1997).
  * [19] V. Mukhanov, ‘‘Physical Foundations of Cosmology,’’ Cambridge University Press, 2005.
  * [20] J. F. Dufaux, G. N. Felder, L. Kofman, M. Peloso and D. Podolsky, JCAP 07, 006 (2006).
  * [21] K. Enqvist, M. Karciauskas, O. Lebedev, S. Rusak and M. Zatta, JCAP 11, 025 (2016).
  * [22] S. Y. Khlebnikov and I. I. Tkachev, Phys. Rev. Lett. 77, 219-222 (1996).
  * [23] T. Prokopec and T. G. Roos, Phys. Rev. D 55, 3768-3775 (1997).
  * [24] D. G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, [arXiv:2102.01031 [astro-ph.CO]].
  * [25] D. G. Figueroa, A. Florio, F. Torrenti and W. Valkenburg, JCAP 04, 035 (2021).
  * [26] O. Lebedev, Prog. Part. Nucl. Phys. 120, 103881 (2021).
  * [27] O. Lebedev, F. Smirnov, T. Solomko and J. H. Yoon, JCAP 10, 032 (2021).
  * [28] P. A. R. Ade et al. [Planck], Astron. Astrophys. 594, A13 (2016).
  * [29] O. Lebedev and J. H. Yoon, Phys. Lett. B 821, 136614 (2021).
  * [30] J. Kost, C. S. Shin and T. Terada, Phys. Rev. D 105, no.4, 043508 (2022).
  * [31] O. Lebedev, T. Nerdi, T. Solomko and J. H. Yoon, [arXiv:2110.05488 [hep-ph]].
  * [32] W. Y. Ai, M. Drewes, D. Glavan and J. Hajer, JHEP 11, 160 (2021).
  * [33] Z. L. Wang and W. Y. Ai, [arXiv:2202.08218 [hep-ph]].
  * [34] O. Lebedev and T. Toma, Phys. Lett. B 798, 134961 (2019).
  * [35] M. Fairbairn, K. Kainulainen, T. Markkanen and S. Nurmi, JCAP 04, 005 (2019).
  * [36] Q. Li, T. Moroi, K. Nakayama and W. Yin, JHEP 09, 179 (2021).


# Density Dependence of the Symmetry Energy
in the Post PREX-CREX Era

Brendan T. Reed<EMAIL_ADDRESS>Theoretical Division, Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA  Department of Astronomy,
Indiana University, Bloomington, Indiana 47405, USA  F. J. Fattoyev
<EMAIL_ADDRESS>Department of Physics & Astronomy, Manhattan
College, Riverdale, NY 10471, USA  C. J. Horowitz<EMAIL_ADDRESS>Facility for Rare Isotope Beams, Michigan State University, East Lansing,
Michigan 48824, USA  Center for Exploration of Energy and Matter and
Department of Physics, Indiana University, Bloomington, IN 47405, USA  J.
Piekarewicz<EMAIL_ADDRESS>Department of Physics, Florida State
University, Tallahassee, FL 32306, USA

###### Abstract

The recently published CREX results suggest a rather peculiar picture for the
density dependence of the symmetry energy. Whereas PREX favors a large neutron
skin thickness in 208Pb, thereby suggesting a stiff equation of state, CREX
suggests instead a much softer equation of state. This discrepancy has caused
a large spur in the theoretical community since no model has been able to
simultaneously reproduce within $1\sigma$ the PREX and CREX results. Motivated
by a novel correlation between a CREX observable and a combination of bulk
symmetry energy parameters, we calibrate three new covariant energy density
functionals that reproduce binding energies and charge radii of spherical
nuclei—and also accommodate the constraints imposed by PREX and CREX. Given
that these models suggest a stiff equation of state at high densities,
predictions for neutron star properties are also discussed.

###### pacs:

21.60.Jz, 21.65.Ef, 24.10.Jv,

##  I Introduction

After almost a decade since first conceived, the Calcium Radius Experiment
(CREX) at Thomas Jefferson National Accelerator Facility (JLab) was recently
completed Adhikari et al. (2022). CREX followed on the footsteps of the
successful Lead Radius Experiment (PREX) Abrahamyan et al. (2012); Horowitz et
al. (2012); Adhikari et al. (2021) that aimed to constrain the equation of
state (EOS) of neutron-star matter in the vicinity of nuclear saturation
density by measuring the neutron skin thickness of 208Pb. Such a powerful
connection between atomic nuclei and neutron stars is encoded in the nuclear
symmetry energy, which quantifies the energy cost of converting symmetric
nuclear matter into pure neutron matter. In the vicinity of nuclear saturation
density, $\rho_{0}\\!\approx\\!0.15\,{\rm fm}^{-3}$, the density dependence of
the symmetry energy is parameterized in terms of a few bulk parameters,

$\displaystyle S(\rho)=J+Lx+\frac{1}{2}K_{\rm{sym}}x^{2}+...\
,\quad\,x\equiv\frac{\rho\\!-\\!\rho_{0}}{3\rho_{0}},$ (1)

where $J$ is the value, $L$ the slope, and $K_{\rm{sym}}\,$ the curvature of
the symmetry energy at saturation density. In particular, the correlation
between the neutron skin thickness of 208Pb and the radius of a canonical
$1.4M_{\odot}\,$ neutron star indicates that $L$ controls both the thickness
of the neutron skin Brown (2000); Furnstahl (2002); Roca-Maza et al. (2011)
and the radius of low-mass neutron stars Horowitz and Piekarewicz (2001a, b);
Carriere et al. (2003); Piekarewicz and Fattoyev (2019).

With tiny parity-violating asymmetries of the order of parts per million
Adhikari et al. (2021, 2022), both CREX and PREX-2 determined the neutral weak
form factors $F_{\rm wk}$ of 48Ca and 208Pb at a single value of the momentum
transfer. These values are as model independent as the ones extracted decades
ago for the corresponding charge form factors $F_{\rm ch}$ using elastic
electron scattering De Vries et al. (1987). By then exploiting the correlation
between the neutron skin thickness $R_{\mathrm{skin}}\\!=\\!R_{n}\\!-\\!R_{p}$
and the weak skin form factor $F_{\mathrm{Wskin}}\\!=\\!F_{\rm
ch}\\!-\\!F_{\rm wk}$ Thiel et al. (2019) predicted by several relativistic
and non-relativistic models, a constraint on the neutron skin is obtained.
Whereas the theoretical model error is small in the case of 208Pb (about 15%
relative to the experimental error), the theoretical and experimental errors
are comparable for 48Ca. Thus, in comparing against experiment it is better to
use the model independent CREX weak skin form factor rather than the neutron
skin. Although part of the larger theoretical error in 48Ca is associated to
the slightly larger than optimal experimental momentum transfer, 48Ca presents
additional theoretical challenges relative to 208Pb. For example, the tensor
component of the current—often referred to as the spin-orbit
contribution—becomes important in 48Ca because in a mean field approximation
the $f_{7/2}$ orbital is full while its $f_{5/2}$ spin-orbit partner is empty
Horowitz and Piekarewicz (2012). Hence, unlike PREX, the connection between
CREX and the density dependence of the symmetry energy is not as robust CRE .

While some differences between PREX and CREX were anticipated, it came as a
surprise to many when CREX reported a central value for the neutron skin
thickness of 48Ca ($R_{\rm skin}^{48}$) that was significantly smaller than
the corresponding value for 208Pb ($R_{\rm skin}^{208}$). Although the
correlation between $R_{\rm skin}^{48}$ and $R_{\rm skin}^{208}$ may not be as
strong as that observed between 208Pb and other heavier neutron-rich nuclei
Piekarewicz et al. (2012), it is difficult to find a theoretical model than
can reproduce both. Indeed, several theoretical approaches have attempted
unsuccessfully to reconcile both measurements Hu et al. (2022); Reinhard et
al. (2022); Zhang and Chen (2022); Mondal and Gulminelli (2022);
Papakonstantinou (2022); Yüksel and Paar (2023). The common theme that has
emerged from these studies suggests that it is difficult to accommodate within
$1\sigma$ the large value of $R_{\rm skin}^{208}$ within the constraints
imposed from other nuclear observables, particularly the electric dipole
polarizability $\alpha_{D}$. Reference Reinhard et al. (2022) goes as far as
suggesting that the large error bars reported by the PREX collaboration makes
it difficult to use the parity violating asymmetry in 208Pb as a meaningful
constraint on the isovector sector of current energy density functionals.
Finally, a recent ab initio approach using chiral forces predicted a value for
$R_{\rm skin}^{208}$ that is in mild tension with PREX Hu et al. (2022), and
for which the correlation between $R_{\rm skin}^{208}$ and $R_{\rm skin}^{48}$
remains fairly strong. Additionally, the range of parameter space inferred by
ab initio calculations point to the large neutron skin of 208Pb being unlikely
Sammarruca (2024). We note that within the ab initio framework, a prediction
of the neutron skin thickness of 48Ca was made well in advance of the
experiment Hagen et al. (2016) that is entirely consistent with the value
reported by the CREX collaboration. It is the goal of this work to explore
refinements to the isovector sector of covariant energy density functionals in
an attempt to reconcile the PREX-2 and CREX results. We first begin by showing
the correlation of several energy density functional (EDF) predictions of
isovector skins with some of the terms of the symmetry energy in section II.
We then extract constraints on these terms inferred from the results of both
PREX and CREX using a Markov Chain Monte Carlo (MCMC) approach in section III.
Next, in section IV, we generate three new relativistic EDFs which are
consistent with the results extracted above. By using these three new EDFs, we
predict several properties of nuclear and neutron star matter (see section V).
We then conclude in section VI.

##  II Isovector Skin Correlations

In a previous letter, we invoked the strong correlation between $R_{\rm
skin}^{208}$ and $L$ to conclude that the symmetry energy is stiff Reed et al.
(2021). While it may be possible to soften the symmetry energy by
incorporating the CREX results, such an approach seems unnatural given that
several recent studies have concluded that the current picture of density
functional theory is insufficient to explain why 208Pb has a thick skin and
48Ca has a thin skin Reinhard et al. (2022); Papakonstantinou (2022); Yüksel
and Paar (2023). What may hold the key to elucidate the source of the
discrepancy are the higher order terms in the density expansion of the
symmetry energy displayed in Eq.(1), particularly $K_{\rm{sym}}\,$. In what
follows, we actually argue that the CREX results may provide the means to
constrain $K_{\rm{sym}}\,$. This is relevant given that earlier studies have
suggested that the properties of finite nuclei with a relative modest neutron-
proton asymmetry are insensitive to $K_{\rm{sym}}\,$ Tsang et al. (2009,
2012); Danielewicz et al. (2017).

We begin by displaying in Fig.1a the well-known correlation between $L$ and
the neutron skin thickness of 208Pb. For this work, we employ the entire set
of models used in the CREX analysis, which includes a fairly comprehensive set
of covariant energy density functionals (blue dots) and a limited set of non-
relativistic Skyrme models (red dots). Predictions are also included for three
new covariant EDFs labeled as “DINOa-c” (colored stars) that are introduced
below and in Table 6. In both Fig.1 and Fig.2 we display 1$\sigma$ theoretical
error bands obtained by including only the covariant EDFs displayed by the
blue dots. Although the correlation between $L$ and the neutron skin thickness
of 48Ca was anticipated not to be as robust as in the case of 208Pb, Fig.1(b)
suggests that for the entire set of covariant EDFs, a very strong correlation
emerges between $K_{\rm sym}\\!-\\!6L$ and the CREX weak skin form factor
$F_{\rm Wskin}^{48}$ evaluated at the average momentum transfer of the
experiment $q\\!=\\!0.8733\,{\rm fm}^{-1}$.

Figure 1: (a) Model predictions for the neutron skin thickness of 208Pb as a
function of $L$ and (b) for the weak skin form factor of 48Ca as a function of
the linear combination $K_{\rm sym}\\!-\\!6L$. The blue and red dots represent
the entire set of both covariant and non-relativistic EDFs used in the CREX
analysis Adhikari et al. (2022). The new DINOa-c interactions are shown as
stars. The experimental values are shown as horizontal bands and the $1\sigma$
errors bands for the covariant EDFs are displayed as solid lines with blue
bands.

This figure suggests—for the set of covariant EDFs employed in this work—the
CREX form factor may best be correlated with linear combinations of the
symmetry energy parameters. In particular, the quantity $K_{\rm{sym}}\,-6L$
appears in the incompressibility coefficient of asymmetric matter
$K(\alpha)=K_{0}+K_{\tau}\,\alpha^{2}+\ldots$, where

$\displaystyle K_{\tau}=K_{\mathrm{sym}}-6L-\frac{Q_{0}}{K_{0}}L,$ (2)

Here $\alpha\\!=\\!(N\\!-\\!Z)/A$, $K_{0}$ is the incompressibility
coefficient of symmetric nuclear matter, and $Q_{0}$ the corresponding
skewness parameter Piekarewicz and Centelles (2009). Note that although
motivated by the incompressibility of asymmetric matter, we find that $F_{\rm
Wskin}^{48}$ is slightly better correlated to $K_{\rm sym}\\!-6L$ than to
$K_{\tau}$ because of the small model dependency coming from the third term.
We also note that whereas the correlation with the covariant EDFs is very
sharp, the non-relativistic models do not follow the same trend. However, the
collection of non-relativistic EDFs used here is not representative, as it
spans a fairly narrow range in $K_{\rm sym}\\!-\\!6L$. Hence, investigating
whether this correlation develops as one includes a comprehensive set of non-
relativistic models is warranted.

As already mentioned and now displayed in Fig.2, neither the correlation
between $F_{\rm Wskin}^{48}$ and $L$ nor the correlation between $R_{\rm
skin}^{208}$ and $K_{\rm sym}\\!-\\!6L$ is as strong as those displayed in
Fig.1. Nevertheless, we show this behavior to clarify some of the details of
the statistical inference that will be implemented later on.

Figure 2: (a) Model predictions for the weak skin form factor of 48Ca as a
function of $L$ and (b) for the neutron skin thickness of 208Pb as a function
of the linear combination $K_{\rm sym}\\!-\\!6L$. The blue and red dots
represent the entire set of both covariant and non-relativistic EDFs used in
the CREX analysis Adhikari et al. (2022). The new DINOa-c interactions are
shown as stars, see Figure 1. The experimental values are shown as horizontal
bands and our linear fits with $1\sigma$ errors are shown as solid lines with
a blue band.

##  III Bayesian Inference of L and $K_{\rm{sym}}\,$

Although we caution against the universality of the correlation between
$F_{\rm Wskin}^{48}\,$ and $K_{\rm{sym}}\,\\!-6L$, we nevertheless proceed to
constrain both $L$ and $K_{\rm{sym}}\,$ using the PREX-2 and CREX results. To
account for a possible model dependence, we consider two different scenarios.
In Scenario 1 we determine $L$ from PREX-2 and then use the newly-found
correlation between $F_{\rm Wskin}^{48}$ and $K_{\rm{sym}}\,-6L$ to extract
$K_{\rm{sym}}\,$. In this scenario only the set of covariant EDFs is used
together with the correlations identified in Fig.1. In Scenario 2 we now add
the limited set of Skyrme forces together with the correlations displayed in
Fig.2. Under this scenario, which greatly increases the model uncertainty, the
1$\sigma$ theoretical error bands are now expanded to include both non-
relativistic and covariant EDFs. We underscore that such analysis should be
repeated using a larger set of Skyrme EDFs, given that the limited set
included here spans a narrow range of values for $K_{\rm{sym}}\,\\!-6L$.

For our calculation of optimal $L$ and $K_{\rm{sym}}\,$ values, we employ the
use of a $\chi^{2}$ minimization routine. The $\chi^{2}$ takes into account
the experimental values and the linear fits shown in Figure 1 and Figure 2.
For the experimental errors, we use the results of the main PREX-2 Adhikari et
al. (2021) and CREX Adhikari et al. (2022) papers. We construct our $\chi^{2}$
using the standard form:

$\chi^{2}=\sum\frac{\Big{(}y_{i}-\hat{y_{i}}(L,K_{\rm
sym})\Big{)}^{2}}{\sigma_{\rm th}^{2}+\sigma_{\rm exp}^{2}}$ (3)

which we minimize using downhill simplex method Virtanen et al. (2020). Here
$y_{i}$ is a set of experimental observables, whose theoretical predictions
$\hat{y}_{i}(L,K_{\rm sym})$ are sensitive to the symmetry energy parameters
$L$, and $K_{\rm sym}$, whereas $\sigma_{\rm th}$ and $\sigma_{\rm exp}$
indicate $1\sigma$ theoretical and experimental errors, respectively.

The list of models we employ in calculating our linear correlations come
directly from the main CREX paper Adhikari et al. (2022), given in Table 1 and
Table 2. Important for our determination of both parameters is the model
uncertainty $\sigma_{\rm th}$ used in the $\chi^{2}$. To do this, we calculate
the 1$\sigma$ prediction interval (PI) from the linear regression line. At
each point along the regression line, the $1\sigma$ PI is calculated using

$\displaystyle\sigma_{\rm
th}=t_{68\%}\times\sqrt{\frac{\sum_{i=1}^{n}(y_{i}-\hat{y}(x_{i}))^{2}}{n-2}}$
(4)

where $t_{68\%}$ is the critical t-value at 68% confidence and the terms in
the square-root represent the standard deviation of the residuals. The
prediction interval gives us an approximate spread of values around the
regression line which indicates the region where $\sim 68\%$ of all models
should fall. For the Scenario 1 fit, we fit the regression line to the
covariant models only for the relations in Figure 1. For the Scenario 2 fit,
we use the regression line fit to both covariant and the Skyrme models.

Model | $L$ | $K_{\rm{sym}}\,$ | $F_{\rm Wskin}^{48}\,$ | $R_{\rm skin}^{208}\,$
---|---|---|---|---
| (MeV) | (MeV) |  | (fm)
FSUGarnet Chen and Piekarewicz (2015a) | 50.97 | 59.45 | 0.0439 | 0.1611
FSUGold2 Chen and Piekarewicz (2015b) | 112.68 | 25.38 | 0.0568 | 0.2863
FSUGold2 (47) Fattoyev et al. (2018) | 47.0 | 53.98 | 0.0424 | 0.1517
FSUGold2 (50) Fattoyev et al. (2018) | 50.0 | 30.23 | 0.0444 | 0.1688
FSUGold2 (54) Fattoyev et al. (2018) | 54.0 | 3.67 | 0.0464 | 0.1864
FSUGold2 (58) Fattoyev et al. (2018) | 58.0 | -17.8 | 0.0479 | 0.2004
FSUGold2 (69) Fattoyev et al. (2018) | 69.0 | -54.97 | 0.0509 | 0.2284
FSUGold2 (76) Fattoyev et al. (2018) | 76.0 | -64.40 | 0.0523 | 0.2416
FSUGold2 (90) Fattoyev et al. (2018) | 90.0 | -55.47 | 0.0544 | 0.2623
FSUGold2 (100) Fattoyev et al. (2018) | 100.0 | -29.14 | 0.0556 | 0.2740
IUFSU Fattoyev et al. (2010) | 47.21 | 28.53 | 0.044 | 0.1615
NL3 Lalazissis et al. (1997) | 118.2 | 100.88 | 0.0555 | 0.2798
RMF022 Chen and Piekarewicz (2015a) | 63.52 | -28.76 | 0.0496 | 0.2164
RMF028 Chen and Piekarewicz (2015a) | 112.65 | 26.26 | 0.0569 | 0.2851
RMF032 Chen and Piekarewicz (2015a) | 125.63 | 28.68 | 0.0587 | 0.3201
TAMUa Fattoyev and Piekarewicz (2013) | 82.49 | -68.37 | 0.0529 | 0.2502
TAMUb Fattoyev and Piekarewicz (2013) | 122.53 | 45.88 | 0.0569 | 0.3002
TAMUc Fattoyev and Piekarewicz (2013) | 135.25 | 51.64 | 0.0583 | 0.3306
IU$\delta$ (-40) Adhikari et al. (2022) | -40.0 | 318 | 0.0206 | 0.0408
IU$\delta$ (-30) Adhikari et al. (2022) | -30.0 | 266 | 0.0239 | 0.0568
IU$\delta$ (-20) Adhikari et al. (2022) | -20.0 | 217 | 0.0271 | 0.0722
IU$\delta$ (-10) Adhikari et al. (2022) | -10.0 | 175 | 0.0302 | 0.0928
IU$\delta$ (0) Adhikari et al. (2022) | 0.0 | 135 | 0.0333 | 0.1133
IU$\delta$ (10) Adhikari et al. (2022) | 10.0 | 100 | 0.0359 | 0.1247
IU$\delta$ (20) Adhikari et al. (2022) | 20.0 | 65 | 0.0388 | 0.1438
IU$\delta$ (30) Adhikari et al. (2022) | 30.0 | 56 | 0.0405 | 0.1510
IU$\delta$ (40) Adhikari et al. (2022) | 40.0 | 1.0 | 0.0442 | 0.1846
Table 1: List of all relativistic energy density functional (EDF) models used in the CREX analysis. Note there are several FSUGold2 and IU-$\delta$ models with different values of the slope parameter of the symmetry energy $L$. Also listed are $R_{\rm skin}^{208}\,=R_{n}-R_{p}$ for 208Pb and $F_{\rm Wskin}^{48}\,=F_{ch}-F_{wk}$ for 48Ca. Table 2: List of all non-relativistic EDF models used in the CREX analysis, similar to Table 1. Model | $L$ | $K_{\rm{sym}}\,$ | $F_{\rm Wskin}^{48}\,$ | $R_{\rm skin}^{208}\,$
---|---|---|---|---
| (MeV) | (MeV) |  | (fm)
SI Skyrme (1956) | 1.22 | -461.85 | 0.0319 | 0.1138
SIII Vautherin and Brink (1972) | 9.91 | -393.74 | 0.0345 | 0.1246
SKM* Bartel et al. (1982) | 45.78 | -155.94 | 0.0392 | 0.1688
SLy4 Chabanat et al. (1998) | 45.96 | -119.7 | 0.0391 | 0.1596
SLy5 Reinhard et al. (1999) | 48.14 | -112.76 | 0.0403 | 0.1622
SLy7 Reinhard et al. (1999) | 47.22 | -113.32 | 0.0394 | 0.1583
SV-K218 Klupfel et al. (2009) | 34.62 | -206.88 | 0.0381 | 0.1622
SV-K226 Klupfel et al. (2009) | 34.09 | -211.92 | 0.0379 | 0.1599
SV-K241 Klupfel et al. (2009) | 30.95 | -230.77 | 0.0379 | 0.1527
SV-bas Klupfel et al. (2009) | 32.36 | -221.76 | 0.0379 | 0.1559
SV-kap00 Klupfel et al. (2009) | 39.44 | -161.78 | 0.0368 | 0.158
SV-kap02 Klupfel et al. (2009) | 35.54 | -193.2 | 0.0373 | 0.1565
SV-kap06 Klupfel et al. (2009) | 29.33 | -249.76 | 0.0384 | 0.1555
SV-mas07 Klupfel et al. (2009) | 52.15 | -98.77 | 0.0389 | 0.1708
SV-mas08 Klupfel et al. (2009) | 40.15 | -172.39 | 0.0383 | 0.1616
SV-mas10 Klupfel et al. (2009) | 28.03 | -252.51 | 0.0375 | 0.1536
SV-sym28 Klupfel et al. (2009) | 7.21 | -296.51 | 0.033 | 0.1178
SV-sym32 Klupfel et al. (2009) | 57.07 | -148.8 | 0.0421 | 0.1933
SV-sym34 Klupfel et al. (2009) | 80.95 | -79.08 | 0.0455 | 0.2287
SV-min Klupfel et al. (2009) | 44.81 | -156.57 | 0.0384 | 0.1716
TOV-min Erler et al. (2013) | 76.23 | -15.62 | 0.0402 | 0.2064
UNEDF0 Kortelainen et al. (2010a) | 45.08 | -189.68 | 0.0445 | 0.1882
UNEDF1 Kortelainen et al. (2010b) | 40.00 | -179.48 | 0.036 | 0.177

Once optimal values are found from the simplex method, the posterior
distribution is found by maximizing the likelihood function, defined as
$\mathcal{L}=\exp{(-\frac{1}{2}\chi^{2})}$, using a Markov Chain Monte Carlo
approach Foreman-Mackey et al. (2013). We use a uniform prior on this approach
for $L$ in the range of $0<L<200$ MeV and leave $K_{\rm{sym}}\,$
unconstrained. Of importance here is to note that the model error is
interpolated over the x-coordinate by taking half the difference of the upper
and lower PI bound at each point along the regression line. Once the posterior
distribution is sampled, one obtains a distribution of model parameters from
which averages and standard deviations may be computed. Our results for both
scenarios can be found in Table 3 with confidence ellipses shown in Figure 3.

Fit | L (MeV) |  $K_{\rm{sym}}\,$ (MeV)
---|---|---
Scenario 1 | $110\pm 40$ | $970\pm 320$
Scenario 2 | $19\pm 19$ | $-61\pm 280$
Table 3: Averages and standard deviations for the two scenarios discussed in
the text. Figure 3: Confidence ellipses for $L$ and $K_{\rm{sym}}\,$ from
Table 3. We show the central values of each fit as a star with the surrounding
ellipse indicating the 67% confidence interval. Notably, the Scenario 1 and
Scenario 2 fits do not overlap.

Our results paint an interesting picture of the density dependence of the
symmetry energy. First, the value of $L\\!=\\!110\pm 40$ MeV extracted from
Scenario 1 is entirely consistent with our original published value of
$L\\!=\\!106\pm 37$ MeV Reed et al. (2021), that is also displayed in the
figure. The small differences are attributed to the larger set of covariant
EDFs used here. In contrast, Scenario 2 favors a considerably smaller value of
$L\\!\approx\\!19\pm 19$ MeV. Such a smaller value is driven by the small
experimental error in $F_{\rm Wskin}^{48}$ relative to $F_{\rm Wskin}^{208}$.
So even if the correlation to $L$ is weaker for the former than for the
latter, the smaller experimental error generates the significantly smaller
value of $L$.

Whereas our knowledge of $L$ has improved since the culmination of the PREX
campaigns, $K_{\rm{sym}}\,$ has remained largely unconstrained, with an
overwhelming number of theoretical approaches favoring negative values. In
contrast, Scenario 1 favors positive values for $K_{\rm{sym}}\,$ at the
$3\sigma$ level. It is only by including the limited set of non-relativistic
models, in combination with the small values of $L$, that $K_{\rm{sym}}\,$ is
allowed to take negative values. We find this possible conflict particularly
interesting. After all, one of the main motivations behind CREX was the use of
ab initio models to inform and improve the isovector sector of energy density
functionals CRE .

##  IV Newly Created Covariant Model Set

However, before incorporating ab initio predictions, we adopt the Scenario 1
PREX-2–CREX constraints in an attempt to generate three new covariant EDFs
with large and positive values for $K_{\rm{sym}}\,$. These models, named
“DINOa,b,c”, are members of the “FSUGold” class of covariant EDFs Horowitz and
Piekarewicz (2001b); Todd-Rutel and Piekarewicz (2005), but with a notable
addition to the isovector sector. In the past, the isovector sector of this
class of models included only two terms: a Yukawa coupling of the nucleon to
the isovector $\rho$-meson and a non-linear $\omega$-$\rho$ isoscalar-
isovector cross-coupling term. In this work we add a contribution to the
isovector sector arising from the Yukawa coupling of the nucleon to the
scalar-isovector $\delta$-meson Kubis et al. (1998); Hofmann et al. (2001);
Liu et al. (2002); Greco et al. (2003); Bunta and Gmuca (2003); Liu et al.
(2005). The presence of the $\delta$-meson splits the nucleon effective mass,
thereby providing an additional degree of freedom in selecting
$K_{\rm{sym}}\,$. To span the range of reasonable values for $L$, we have
defined the DINOa, DINOb, and DINOc, models by fixing the slope of the
symmetry energy at $L\\!=\\!50,70,90$ MeV, respectively. Our fitting procedure
follows closely that of Ref. Fattoyev et al. (2010), supplemented by the
inclusion of both $R_{\rm skin}^{208}$ and $F_{\rm Wskin}^{48}\,$ into the
calibration.

We begin our construction with the class of FSUGold-like models featuring the
delta meson Todd-Rutel and Piekarewicz (2005); Singh et al. (2014). The
interaction Lagrangian for this class of models is

$\displaystyle\mathcal{L}_{\mathrm{int}}$
$\displaystyle=\bar{\psi}\Big{[}\mathcal{S}(\phi,\bm{\delta})-\mathcal{V}_{\mu}(V_{\mu},\bm{b_{\mu}},A_{\mu})\gamma^{\mu}\Big{]}\psi-U(\phi)$
$\displaystyle+\frac{\zeta}{4!}g_{v}^{4}(V_{\mu}V^{\mu})^{2}+\Lambda_{v}g_{v}^{2}g_{\rho}^{2}V_{\mu}V^{\mu}\bm{b_{\mu}}\cdot\bm{b^{\mu}}$

where

$\displaystyle\mathcal{S}(\phi,\bm{\delta})=g_{s}\phi+\frac{g_{\delta}}{2}\bm{\tau}\cdot\bm{\delta}$
(6)
$\displaystyle\mathcal{V}_{\mu}(V_{\mu},\bm{b_{\mu}},A_{\mu})=g_{v}V_{\mu}+\frac{g_{\rho}}{2}\bm{\tau}\cdot\bm{b_{\mu}}+\frac{e}{2}(1+\tau_{3})A_{\mu}$
$\displaystyle
U(\phi)=\frac{\kappa}{3!}(g_{s}\phi)^{3}+\frac{\lambda}{4!}(g_{s}\phi)^{4}$

This class of models contain the interactions of nucleons ($\psi$) via the
exchange of the scalar $\sigma$ ($\phi$) and vector $\omega$ ($V_{\mu}$)
mesons as well as the vector-isovector $\rho$-meson ($\textbf{b}^{\mu}$), and,
in this work, the additional vector-isoscalar $\delta$-meson ($\bm{\delta}$).
In addition to interactions among the nucleons, this Lagrangian class also
includes many self and cross-coupling terms. These are $\kappa$ and $\lambda$
for the cubic and quartic scalar self-coupling, $\zeta$ for the quartic
$\omega$ self-coupling, and $\Lambda_{v}$ as the $\omega$-$\rho$ cross
coupling.

Motivated by our findings of large positive $K_{\rm{sym}}\,$ value from the
Scenario 1 fit, we choose to generate three models which we expect run the
gamut on parameter space. The generation of new models follows a similar
prescription as described in Fattoyev and Piekarewicz (2010) which we
summarize here briefly. For each model, we fix the value of the symmetry
energy at sub-saturation density $k_{F}=1.15$ fm-1 to be $\tilde{J}=27$ MeV.
This value has been shown to agree well with measurements of nuclear masses.
We then fix $L$ for the three models at 50, 70, and 90 MeV to give a good
range of $L$ values favored from multiple sources including PREX-2. We also
fix the $\zeta$-parameter to $0.015$, which controls the high-density
components of the EOS. This allows our models to predict neutron star masses
consistent with the most recent observational and theoretical constraints on
the maximum mass Demorest et al. (2010); Antoniadis et al. (2013); Rezzolla et
al. (2018); Miller et al. (2021); Riley et al. (2021); Romani et al. (2022).
From here, we take a typical set of bulk nuclear matter properties, and adjust
the isoscalar bulk properties to fit the binding energy and charge radii (with
spin-orbit corrections Horowitz and Kadau (2009)) of 40Ca , 48Ca , and 208Pb .

The inversion of bulk properties to RMF Lagrangian couplings has been
described in Chen and Piekarewicz (2015b). It was found that the Giant
Monopole Resonances (GMR) are sensitive primarily to the incompressibility of
neutron-rich matter, $K(\alpha)$ Piekarewicz (2010). In our previous work
Fattoyev and Piekarewicz (2013) we found that for 208Pb this corresponds to
$K(\alpha)\approx 222$ MeV for the RMF parameterizations that yield consistent
GMRs. Hence we applied this constraint to $K(\alpha)$ in the DINO models. All
that is left now is to pick a value for $K_{\rm{sym}}\,$ which minimizes the
squared difference of $R_{\rm skin}^{208}\,$ and $F_{\rm Wskin}^{48}\,$ from
their experimental values. The difference between the model predictions and
experimental results for each observable are then used to generate a new set
of parameters until a parameter set converges.

It should be noted that a direct inversion of couplings from $\tilde{J}$, $L$,
and $K_{\rm{sym}}\,$ does not presently exist as the presence of the
additional $\delta$-meson introduces much more complexity than has previously
been detailed. As such, we use a downhill simplex routine to determine the
isovector couplings $g_{\rho}^{2}$, $g_{\delta}^{2}$, and $\Lambda_{v}$.
Furthermore, there exist some combinations of bulk nuclear properties that do
not produce a convergent set of parameters which requires careful calculation
schemes.

The coupling constants of the three DINO models can be found in Table 4. We
find that due to the large value of $K_{\rm{sym}}\,$, the isovector Yukawa
couplings $g_{\rho}^{2}$ and $g_{\delta}^{2}$ are very large. Additionally, we
report the bulk parameters of the nuclear EOS in Table 5 as predicted by these
three models. In the next section, using these new parameters set we proceed
to calculate the properties of finite nuclei and neutron stars.

Model | $m_{\rm s}$ | $g_{s}^{2}$ | $g_{\delta}^{2}$ | $g_{v}^{2}$ | $g_{\rho}^{2}$ | $\kappa$ | $\lambda$ | $\zeta$ | $\Lambda_{\rm v}$
---|---|---|---|---|---|---|---|---|---
DINOa | $490.050$ | $93.9422$ | $1115.15$ | $154.436$ | $805.891$ | $4.9860$ | $-0.01370$ | $0.015$ | $0.0016497$
DINOb | $485.795$ | $91.0316$ | $1252.71$ | $150.824$ | $877.121$ | $5.2914$ | $-0.01488$ | $0.015$ | $0.0014014$
DINOc | $484.162$ | $90.6481$ | $1343.25$ | $151.048$ | $922.617$ | $5.3209$ | $-0.01497$ | $0.015$ | $0.0012312$
Table 4: Coupling constants of the DINO set of models. For each of the three models, we fix $m_{\rm v}=782.5$ MeV, $m_{\rho}=763$ MeV, and $m_{\delta}=980$ MeV as well as the quartic vector self-coupling $\zeta=0.015$. Scalar meson mass and nonlinear $\kappa$ values given in MeV. Model |  $\rho_{0}$ (fm-3) |  $\epsilon_{0}$ (MeV) |  $K_{0}$ (MeV) |  $Q_{0}$ (MeV) |  $\tilde{J}$ (MeV) |  $J$ (MeV) |  $L$ (MeV) |  $K_{\rm sym}$ (MeV) |  $K_{\tau}$ (MeV) |  $K(\alpha)$ (MeV)
---|---|---|---|---|---|---|---|---|---|---
DINOa | $0.1522$ | $-16.16$ | $210.0$ | $-361.4$ | $27.00$ | $31.42$ | $50.00$ | $506.0$ | $292.1$ | $223.1$
DINOb | $0.1525$ | $-16.21$ | $207.0$ | $-412.0$ | $27.00$ | $33.07$ | $70.00$ | $610.0$ | $329.3$ | $221.7$
DINOc | $0.1519$ | $-16.22$ | $206.0$ | $-426.1$ | $27.00$ | $34.58$ | $90.00$ | $715.0$ | $361.2$ | $222.2$
Table 5: Bulk properties of nuclear matter for our set of DINO models. Tilde denotes symmetry energy values taken at subsaturation density ($\tilde{k_{F}}=1.15$fm-1). Here $K(\alpha)$ is calculated using the equation 18c of Piekarewicz and Centelles (2009) for $\alpha=(126-82)/208$. Figure 4: Predictions for the weak skin form factor of 208Pb and 48Ca for the entire set of covariant RMF (green squares) and non-relativistic Skyrme (pink circles) EDFs considered in this letter. The blue ellipses represent joint PREX-2 and CREX 67% and 90% probability contours, just as in Figure 2 of Ref.Adhikari et al. (2022). The FSUGold2 prediction is included to illustrate typical statistical uncertainties. Depicted with stars are the prediction of the three new DINO models. Nucleus | Observable | Experiment | IUFSU | FSUGold2 | BigApple | DINOa | DINOb | DINOc
---|---|---|---|---|---|---|---|---
|  $B/A$ (MeV) | 8.667 | 8.549 | 8.621 | 8.531 | 8.670 | 8.667 | 8.667
|  $R_{\rm ch}$ (fm) | $3.477\pm 0.002$ | 3.416 | 3.413 | 3.447 | 3.454 | 3.458 | 3.461
48Ca |  $R_{\rm skin}$ (fm) | $0.121\pm 0.035$ | 0.1731 | 0.2319 | 0.1682 | 0.0996 | 0.1051 | 0.1135
| $F_{\rm Wskin}\,(q_{\mathrm{CREX}})$ | $0.0277\pm 0.0055$ | 0.0439 | 0.0568 | 0.0417 | 0.0333 | 0.0349 | 0.0367
|  $B/A$ (MeV) | 7.867 | 7.896 | 7.872 | 7.872 | 7.867 | 7.867 | 7.867
|  $R_{\rm ch}$ (fm) | $5.501\pm 0.001$ | 5.476 | 5.489 | 5.490 | 5.502 | 5.503 | 5.503
208Pb |  $R_{\rm skin}$ (fm) | $0.283\pm 0.071$ | 0.1615 | 0.2863 | 0.1506 | 0.1748 | 0.1993 | 0.2240
| $F_{\rm Wskin}\,(q_{\mathrm{PREX}})$ | $0.041\pm 0.013$ | 0.0233 | 0.0423 | 0.0215 | 0.0263 | 0.0303 | 0.0344
Table 6: Experimental data for the binding energy per nucleon (in MeV) Wang et
al. (2012), charge radii (in fm) Angeli and Marinova (2013), neutron skins,
form factor skins Adhikari et al. (2021, 2022) for 48Ca and 208Pb nuclei used
in the optimization. Also displayed are the theoretical results obtained with
IUFSU Fattoyev et al. (2010), FSUGold2 Chen and Piekarewicz (2014), BigApple
Fattoyev et al. (2020), and the three new parameterizations.

##  V DINO Model Predictions

###  V.1 Predictions of Nuclei

We start this section by mentioning that the new model parameters, together
with predictions for the bulk properties of infinite nuclear matter are given
Table 4 and Table 5, respectively. It is worth mentioning, however, that in
all three cases, the curvature of the symmetry energy is large and positive,
namely, $K_{\rm sym}\\!=\\!506,610,715$ MeV, for DINOa, DINOb, and DINOc,
respectively. Regardless, we find that the results for the binding energy,
charge radii, and form factors of 48Ca and 208Pb as listed in Table 6 are all
in good agreement with the experiment. In the particular case of $F_{\rm
Wskin}^{208}\,$ and $F_{\rm Wskin}^{48}\,$, predictions for the large set of
models considered here alongside the predictions from the three DINO models
are displayed in Fig.4. This figure highlights the predicament faced by most
energy density functionals: with few exceptions, most models fall outside the
90% confidence ellipse. In this sense, the three DINO models do quite well
compared to the other models considered.

However, the DINO models’ large isovector couplings do begin to introduce
anomalies in the point particle densities. In particular, the proton and
neutron densities in the core of both 48Ca and 208Pb show large density
fluctuations that are non-physical. For example, 208Pb is a heavy closed-shell
nucleus and so we expect the nucleus to saturate at a density value near
$\rho_{0}$ De Vries et al. (1987); Horowitz et al. (2020). For the DINO
models, 208Pb has very large fluctuations and does not saturate. This anomaly
is largely due to the isovector $g_{\delta}^{2}$ coupling splitting the proton
and neutron masses, thereby causing a larger proton density in the core and a
small neutron density. We see similar fluctuations in 48Ca , although the
saturation property is less strict for lighter nuclei.

Figure 5: Point neutron (solid) and proton (dashed) densities in 208Pb for the
DINO family of models. Note that the central proton density for DINOc is shown
off the plot. Figure 6: Point neutron (solid) and proton (dashed) densities in
48Ca for the DINO family of models.

###  V.2 Predictions of Neutron Star Properties

Our calculation of neutron star (NS) structure begins by following the
standard practices of relativistic mean field theory to calculate the equation
of state of the core Glendenning (2000); Fattoyev (2011); Singh et al. (2014).
This includes calculating the pressure and energy density of protons,
neutrons, electrons, and muons in beta equilibrium using the interaction
Lagrangian of section IV. The outer crust is tabulated using the mass table of
Duflo and Zuker Duflo and Zuker (1995) up to neutron drip $\sim 10^{-4}$ fm-3.
The outer crust and core equations of state are then interpolated using a
cubic polynomial ensuring both the density and sound speed are continuous at
each interface as in Piekarewicz and Fattoyev (2019). We should note that the
prescription for the inner crust may account for as much as 0.5km difference
depending on the formalism used. However, the tidal deformability is
insensitive to the crust EOS Piekarewicz and Fattoyev (2019). The inner crust-
core transition density is determined using the RPA-dynamical method Baym and
Flowers (1974); Carriere et al. (2003) with a slight modification to accompany
the presence of the delta meson Reed (2023).

Although the DINO models are consistent with the properties of finite nuclei
(see Table 6) their large values for $K_{\rm{sym}}\,$ induce a dramatic
stiffening of the EOS at the high densities of relevance to neutron stars.
Such a stiffening generates neutron star radii that are too large when
compared against LIGO-Virgo and NICER data. This behavior is illustrated in
Fig.7, with explicit values in Table 7, which displays the so called “Holy
Grail” of neutron star physics: the mass-radius relation. Together with the
predictions from the three DINO models and a few other covariant EDFs, the
figure shows mass-radius determinations by the NICER mission for the two
pulsars PSR J0030+045 Miller et al. (2019); Riley et al. (2019) and PSR
J0740+6620 Miller et al. (2021); Riley et al. (2021).

With radii of about 15 km, the predictions from all three DINO models are
inconsistent with the NICER data for the low mass star. With such large radii,
the models also predict very large tidal deformabilities that are highly
disfavored by GW170817 Abbott et al. (2017); Malik et al. (2018); Most et al.
(2018); Radice and Dai (2019); Tews et al. (2018). In particular, we find that
DINO models predict larger neutron star radii and lower densities for the
direct Urca threshold. To agree with radii observation, one expects that a
phase transition or other strong density dependence develops at intermediate
densities. Moreover, a very small density for the direct Urca threshold
suggest that all neutron stars with normal matter would undergo fast cooling
which does not agree with to the cooling observations Beznogov and Yakovlev
(2015). If these predictions stand true, it raises the intriguing possibility
of a neutron star core permeated with superfluid or superconductive matter
throughout. Or in conjunction with the radii observation, this may even imply
the prospect of an earlier onset of a phase transition to exotic states.
Hence, we must conclude that while the DINO models provide a plausible
solution to the PREX-2–CREX dilemma, they fail to reproduce astrophysical
observations—unless a phase transition or other strong density dependence
develops at intermediate densities. The emergence of a phase
transition—normally accompanied by a softening of the EOS—has been shown to
reduce the radius and deformability of neutron stars Chatziioannou (2022); Han
and Steiner (2018); Xie and Li (2021).

Model | $M_{\rm max}$ | $R_{1.4}$ | $\Lambda_{1.4}$ | $\rho_{\rm t}$ | $\rho_{\rm Urca}$ | $M_{\rm Urca}$
---|---|---|---|---|---|---
| (M⊙) | (km) |  | (fm-3) | (fm-3) | (M⊙)
DINOa | 2.17 | 14.82 | 1050.6 | 0.0914 | 0.1580 | 0.418
DINOb | 2.15 | 15.11 | 1128.2 | 0.0846 | 0.1438 | 0.427
DINOc | 2.14 | 15.41 | 1240.4 | 0.0789 | 0.1373 | 0.458
Table 7: Prediction of various neutron star observables for our models. We
calculate the core-crust transition density as described in Carriere et al.
(2003) with the modification to account for the addition of $\delta$-mesons.
Figure 7: Mass-radius relation for neutron stars. Along with our three DINO
models, predictions are also shown for a few covariant energy density
functionals. The two sets of contours display the NICER mass-radius
constraints of two pulsars, PSR J0030+0451 Miller et al. (2019); Riley et al.
(2019) and PSR J0740+6620 Miller et al. (2021); Riley et al. (2021) at 68% and
95% confidence.

##  VI Conclusions

In conclusion, the recently measured weak skin form factor in 48Ca by the CREX
collaboration paints a peculiar picture for the behavior of neutron rich
matter. When combined with the PREX-2 result that favors a large value for
$L$, the newly-revealed correlation between $F_{\rm Wskin}^{48}$ and the
combination of bulk symmetry energy parameters $K_{\rm sym}\\!-6L$, motivates
the calibration of a new set of covariant EDFs with large and positive values
for $K_{\rm sym}$. This result is in striking contrast to most non-
relativistic EDFs and ab initio models that instead predict negative values
for $K_{\rm sym}$. Nevertheless, the newly created DINO models with such
unconventional values for $K_{\rm sym}$ agree better with PREX and CREX than
many EDFs existent in the literature. However, this accomplishment comes at
the price of stiffening the equation of state at the high densities of
relevance to neutron stars, resulting in a failure to accommodate the recent
NICER constraints. Nonetheless, these models represent a first step in
reconciling PREX and CREX for covariant EDFs.

Moving forward, the Mainz Radius Experiment (MREX) at the future Mainz Energy-
recovery Superconducting Accelerator (MESA) Becker et al. (2018) promises to
increase the precision of $R_{\rm skin}^{208}$ by a factor of two. MREX will
then confirm whether the slope of the symmetry energy is indeed stiff, or if
the PREX-2 measurement represents a large statistical fluctuation. Further,
that most theoretical models disfavor large and positive values for $K_{\rm
sym}$ may suggest that the strong correlation uncovered here may be model
dependent. Should this be the case, then ab initio models could inform how to
improve the isovector sector of energy density functionals CRE . Finally,
although the stiff symmetry energy suggested by the DINO models at high
densities are ruled out by observation, these models could still be fruitful
in describing the properties of atomic nuclei—especially if a phase transition
emerges at intermediate densities. Regardless, the PREX-2–CREX tension
continues to reenergize the theoretical, experimental, and observational
communities in our common quest to understand the behavior of dense, neutron
rich matter.

###### Acknowledgements.

We would like to thank Anna Watts for providing us with the NICER contours.
This material is based upon work supported by the U.S. Department of Energy
Office of Science, National Nuclear Security Administration of U.S. Department
of Energy (Contract No. 89233218CNA000001), the Laboratory Directed Research
and Development program of Los Alamos National Laboratory under project number
20230785PRD1, Office of Nuclear Physics under Awards DE-FG02-87ER40365
(Indiana University), Number DE-FG02-92ER40750 (Florida State University), and
Number DE-SC0008808 (NUCLEI SciDAC Collaboration).

## References

  * Adhikari et al. (2022) D. Adhikari et al. (CREX), Phys. Rev. Lett. 129, 042501 (2022).
  * Abrahamyan et al. (2012) S. Abrahamyan, Z. Ahmed, H. Albataineh, K. Aniol, D. S. Armstrong, et al., Phys. Rev. Lett. 108, 112502 (2012).
  * Horowitz et al. (2012) C. J. Horowitz, Z. Ahmed, C. M. Jen, A. Rakhman, P. A. Souder, et al., Phys. Rev. C85, 032501 (2012).
  * Adhikari et al. (2021) D. Adhikari et al. (PREX), Phys. Rev. Lett. 126, 172502 (2021).
  * Brown (2000) B. A. Brown, Phys. Rev. Lett. 85, 5296 (2000).
  * Furnstahl (2002) R. J. Furnstahl, Nucl. Phys. A706, 85 (2002).
  * Roca-Maza et al. (2011) X. Roca-Maza, M. Centelles, X. Viñas, and M. Warda, Phys. Rev. Lett. 106, 252501 (2011).
  * Horowitz and Piekarewicz (2001a) C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001a).
  * Horowitz and Piekarewicz (2001b) C. J. Horowitz and J. Piekarewicz, Phys. Rev. C64, 062802 (2001b).
  * Carriere et al. (2003) J. Carriere, C. J. Horowitz, and J. Piekarewicz, Astrophys. J. 593, 463 (2003).
  * Piekarewicz and Fattoyev (2019) J. Piekarewicz and F. J. Fattoyev, Physics Today 72, 30 (2019).
  * De Vries et al. (1987) H. De Vries, C. W. De Jager, and C. De Vries, Atom. Data Nucl. Data Tabl. 36, 495 (1987).
  * Thiel et al. (2019) M. Thiel, C. Sfienti, J. Piekarewicz, C. J. Horowitz, and M. Vanderhaeghen, J. Phys. G46, 093003 (2019).
  * Horowitz and Piekarewicz (2012) C. J. Horowitz and J. Piekarewicz, Phys.Rev. C86, 045503 (2012).
  * (15) _CREX: Parity-violating measurement of the weak charge distribution of 48Ca_, http://hallaweb.jlab.org/parity/prex/c-rex/c-rex.pdf.
  * Piekarewicz et al. (2012) J. Piekarewicz, B. Agrawal, G. Colò, W. Nazarewicz, N. Paar, et al., Phys. Rev. C85, 041302(R) (2012).
  * Hu et al. (2022) B. Hu et al., Nature Phys. 18, 1196 (2022).
  * Reinhard et al. (2022) P.-G. Reinhard, X. Roca-Maza, and W. Nazarewicz (2022), eprint 2206.03134.
  * Zhang and Chen (2022) Z. Zhang and L.-W. Chen (2022), eprint 2207.03328.
  * Mondal and Gulminelli (2022) C. Mondal and F. Gulminelli (2022), eprint 2209.05177.
  * Papakonstantinou (2022) P. Papakonstantinou (2022), eprint 2210.02696.
  * Yüksel and Paar (2023) E. Yüksel and N. Paar, Phys. Lett. B 836, 137622 (2023), eprint 2206.06527.
  * Sammarruca (2024) F. Sammarruca, Symmetry 16 (2024), ISSN 2073-8994, URL https://www.mdpi.com/2073-8994/16/1/34.
  * Hagen et al. (2016) G. Hagen et al., Nature Phys. 12, 186 (2016).
  * Reed et al. (2021) B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, Phys. Rev. Lett. 126, 172503 (2021).
  * Papakonstantinou (2022) P. Papakonstantinou, arXiv e-prints arXiv:2210.02696 (2022), eprint 2210.02696.
  * Tsang et al. (2009) M. B. Tsang et al., Phys. Rev. Lett. 102, 122701 (2009).
  * Tsang et al. (2012) M. Tsang, J. Stone, F. Camera, P. Danielewicz, S. Gandolfi, et al., Phys.Rev. C86, 015803 (2012).
  * Danielewicz et al. (2017) P. Danielewicz, P. Singh, and J. Lee, Nucl. Phys. A 958, 147 (2017).
  * Piekarewicz and Centelles (2009) J. Piekarewicz and M. Centelles, Phys. Rev. C79, 054311 (2009).
  * Virtanen et al. (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al., Nature Methods 17, 261 (2020).
  * Chen and Piekarewicz (2015a) W.-C. Chen and J. Piekarewicz, Phys. Rev. Lett. 115, 161101 (2015a).
  * Chen and Piekarewicz (2015b) W.-C. Chen and J. Piekarewicz, Phys. Lett. B748, 284 (2015b).
  * Fattoyev et al. (2018) F. J. Fattoyev, J. Piekarewicz, and C. J. Horowitz, Phys. Rev. Lett. 120, 172702 (2018).
  * Fattoyev et al. (2010) F. J. Fattoyev, C. J. Horowitz, J. Piekarewicz, and G. Shen, Phys. Rev. C82, 055803 (2010).
  * Lalazissis et al. (1997) G. A. Lalazissis, J. Konig, and P. Ring, Phys. Rev. C55, 540 (1997).
  * Fattoyev and Piekarewicz (2013) F. J. Fattoyev and J. Piekarewicz, Phys. Rev. Lett. 111, 162501 (2013).
  * Skyrme (1956) T. H. R. Skyrme, Phil. Mag. 1, 1043 (1956).
  * Vautherin and Brink (1972) D. Vautherin and D. M. Brink, Phys. Rev. C5, 626 (1972).
  * Bartel et al. (1982) J. Bartel, P. Quentin, M. Brack, C. Guet, and H.-B. Håkansson, Nuclear Physics A 386, 79 (1982), ISSN 0375-9474, URL https://www.sciencedirect.com/science/article/pii/0375947482904031.
  * Chabanat et al. (1998) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nuclear Physics A 635, 231 (1998), ISSN 0375-9474, URL https://www.sciencedirect.com/science/article/pii/S0375947498001808.
  * Reinhard et al. (1999) P.-G. Reinhard, D. J. Dean, W. Nazarewicz, J. Dobaczewski, J. A. Maruhn, and M. R. Strayer, Phys. Rev. C 60, 014316 (1999), URL https://link.aps.org/doi/10.1103/PhysRevC.60.014316.
  * Klupfel et al. (2009) P. Klupfel, P.-G. Reinhard, T. Burvenich, and J. Maruhn, Phys.Rev. C79, 034310 (2009), eprint 0804.3385.
  * Erler et al. (2013) J. Erler, C. J. Horowitz, W. Nazarewicz, M. Rafalski, and P.-G. Reinhard, Phys. Rev. C87, 044320 (2013).
  * Kortelainen et al. (2010a) M. Kortelainen, R. Furnstahl, W. Nazarewicz, and M. Stoitsov, Phys. Rev. C82, 011304 (2010a).
  * Kortelainen et al. (2010b) M. Kortelainen, T. Lesinski, J. More, W. Nazarewicz, J. Sarich, et al., Phys. Rev. C82, 024313 (2010b).
  * Foreman-Mackey et al. (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, PASP 125, 306 (2013), eprint 1202.3665.
  * Todd-Rutel and Piekarewicz (2005) B. G. Todd-Rutel and J. Piekarewicz, Phys. Rev. Lett 95, 122501 (2005).
  * Kubis et al. (1998) S. Kubis, M. Kutschera, and S. Stachniewicz, Acta Phys. Polon. B 29, 809 (1998), eprint astro-ph/9802303.
  * Hofmann et al. (2001) F. Hofmann, C. M. Keil, and H. Lenske, Phys. Rev. C 64, 034314 (2001), eprint nucl-th/0007050.
  * Liu et al. (2002) B. Liu, V. Greco, V. Baran, M. Colonna, and M. Di Toro, Phys.Rev. C65, 045201 (2002), eprint nucl-th/0112034.
  * Greco et al. (2003) V. Greco, V. Baran, M. Colonna, M. Di Toro, T. Gaitanos, and H. H. Wolter, Phys. Lett. B 562, 215 (2003), eprint nucl-th/0212102.
  * Bunta and Gmuca (2003) J. K. Bunta and S. Gmuca, Phys. Rev. C 68, 054318 (2003), eprint nucl-th/0309046.
  * Liu et al. (2005) B. Liu, H. Guo, M. Di Toro, and V. Greco, Eur. Phys. J. A 25, 293 (2005), eprint nucl-th/0409014.
  * Singh et al. (2014) S. K. Singh, S. K. Biswal, M. Bhuyan, and S. K. Patra, Phys. Rev. C 89, 044001 (2014).
  * Fattoyev and Piekarewicz (2010) F. J. Fattoyev and J. Piekarewicz, Phys. Rev. C82, 025805 (2010).
  * Demorest et al. (2010) P. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, Nature 467, 1081 (2010).
  * Antoniadis et al. (2013) J. Antoniadis, P. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, et al., Science 340, 6131 (2013).
  * Rezzolla et al. (2018) L. Rezzolla, E. R. Most, and L. R. Weih, Astrophys. J. Lett. 852, L25 (2018), eprint 1711.00314.
  * Miller et al. (2021) M. C. Miller et al., Astrophys. J. Lett. 918, L28 (2021).
  * Riley et al. (2021) T. E. Riley et al., Astrophys. J. Lett. 918, L27 (2021).
  * Romani et al. (2022) R. W. Romani, D. Kandel, A. V. Filippenko, T. G. Brink, and W. Zheng, Astrophys. J. Lett. 934, L18 (2022).
  * Horowitz and Kadau (2009) C. J. Horowitz and K. Kadau, Phys. Rev. Lett. 102, 191102 (2009).
  * Piekarewicz (2010) J. Piekarewicz, J. Phys. G37, 064038 (2010).
  * Wang et al. (2012) M. Wang, G. Audi, A. H. Wapstra, F. G. Kondev, M. MacCormick, X. Xu, and B. Pfeiffer, Chinese Phys. C 36, 1603 (2012).
  * Angeli and Marinova (2013) I. Angeli and K. Marinova, At. Data Nucl. Data Tables 99, 69 (2013).
  * Chen and Piekarewicz (2014) W.-C. Chen and J. Piekarewicz, Phys. Rev. C90, 044305 (2014).
  * Fattoyev et al. (2020) F. J. Fattoyev, C. J. Horowitz, J. Piekarewicz, and B. Reed, Phys. Rev. C 102, 065805 (2020).
  * Horowitz et al. (2020) C. J. Horowitz, J. Piekarewicz, and B. Reed, Phys. Rev. C 102, 044321 (2020).
  * Glendenning (2000) N. K. Glendenning, _Compact Stars_ (Springer-Verlag New York, 2000).
  * Fattoyev (2011) F. Fattoyev, _Sensitivity of neutron star properties to the equation of state_ (Florida State University, Tallahassee, 2011).
  * Duflo and Zuker (1995) J. Duflo and A. Zuker, Phys. Rev. C52, R23 (1995).
  * Baym and Flowers (1974) G. Baym and E. Flowers, Nucl. Phys. A 222, 29 (1974).
  * Reed (2023) B. T. Reed, Ph.D. thesis (2023), URL https://www.proquest.com/dissertations-theses/neutron-stars-skins-crossroads-nuclear/docview/2841238132/se-2.
  * Miller et al. (2019) M. C. Miller et al., Astrophys. J. Lett. 887, L24 (2019).
  * Riley et al. (2019) T. E. Riley et al., Astrophys. J. Lett. 887, L21 (2019).
  * Abbott et al. (2017) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017).
  * Malik et al. (2018) T. Malik, N. Alam, M. Fortin, C. Providencia, B. K. Agrawal, T. K. Jha, B. Kumar, and S. K. Patra, Phys. Rev. C98, 035804 (2018).
  * Most et al. (2018) E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. 120, 261103 (2018).
  * Radice and Dai (2019) D. Radice and L. Dai, Eur. Phys. J. A 55, 50 (2019).
  * Tews et al. (2018) I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C98, 045804 (2018).
  * Beznogov and Yakovlev (2015) M. V. Beznogov and D. G. Yakovlev, Mon. Not. Roy. Astron. Soc. 452, 540 (2015), eprint 1507.04206.
  * Chatziioannou (2022) K. Chatziioannou, Phys. Rev. D 105, 084021 (2022).
  * Han and Steiner (2018) S. Han and A. W. Steiner (2018), eprint 1810.10967.
  * Xie and Li (2021) W.-J. Xie and B.-A. Li, Phys. Rev. C 103, 035802 (2021).
  * Becker et al. (2018) D. Becker et al. (2018), eprint 1802.04759.