# Characterizing the superconducting instability in a two-orbital $d$-$s$
model:
insights to infinite-layer nickelate superconductors

Mi Jiang  Institute of Theoretical and Applied Physics, Jiangsu Key Laboratory
of Thin Films, School of Physical Science and Technology, Soochow University,
Suzhou 215006, China

###### Abstract

Motivated by the recent realization of the infinite-layer nickelate
superconductivity (SC), we quantitatively investigate a two-orbital $d$-$s$
model by dynamic cluster quantum Monte Carlo calculations. Focusing on the
impact of inter-orbital hybridization on the superconducting properties, our
simulations indicate that the $d$-$s$ hybridization strength has a decisive
role in the suppression of the $d$-wave pairing in the doping regime relevant
to infinite-layer nickelates. Although we confirm on the single-orbital
description at weak hybridization, at relatively large hybridization, there
exists a constructive effect of the non-negligibly finite $s$ orbital
occupancy on inhibiting the suppression of the superconductivity. More
strikingly, there exists a possible SC enhancement at large enough Hubbard
interaction via large $d$-$s$ hybridization. We further provide some insights
on the relevance of $d$-$s$ model to infinite-layer nickelates.

##  I Introduction

The realization of the superconductivity (SC) in Sr-doped NdNiO2 films grown
on a SrTiO3 substrate 2019Nature  is the breakthrough in the decade-long
pursuit of rare-earth nickelates Junjie_review ; Botana_review . Even more
fascinating is the successful synthesis of superconducting doped Pr1-xSrxNiO2
Pr , La1-xCaxNiO2 LaSC , and undoped quintuple-layer Nd6Ni5O12 Nd6Ni5O12 ,
which indicates the advent of Ni-based superconducting family. Undoubtedly,
the most important question is the pairing symmetry and mechanism in this new
family of superconductors and furthermore its similarity and difference from
other transition metal based superconductors like Cu- and Fe-based compounds.
For example, given that the parent compound of infinite-layer NdNiO2 is
discovered not to host the long-range magnetic order Pickett2004 ; Hepting ;
dft1 ; dft9 ; dft16 ; dft18 ; dft26 ; dft25 ; noLRO  in spite of magnetic
correlations mag1 ; mag2 ; mag3 ; mag4 ; mag5 , whether the debated spin-
fluctuation scenario for cuprate superconductors also applies for the
infinite-layer nickelates is an intriguing question.

Physically, the minimal microscopic model to capture the essential physics of
infinite-layer nickelates is a demanding task. Previous numerous electronic
structure calculations based on DFT Pickett2004 ; dft1 ; dft2 ; dft3 ; dft4 ;
dft5 ; dft6 ; dft7 ; dft8 ; dft9 ; dft10 ; dft24 ; dft26  and/or DFT+DMFT
dft101 ; dft102 ; dft11 ; dft12 ; dft13 ; dft14 ; dft15 ; dft16 ; dft17 ;
dft18 ; dft19 ; dft20 ; dft21 ; dft22 ; dft23 ; dft25 ; Held2020 ; Held2022 ;
Hanghuireview  methodologies have provided fruitful understanding of the
electronic structure of infinite-layer nickelates including Rn+1NinO2n+2 with
R = La, Nd, Pr for varying $n$ Botana_review ; Hanghuireview ; Held2022 ,
whose calculated properties have been partly confirmed by recent experiments
for $n=5$ Nd6Ni5O12 . One distinct feature associated with infinite-layer
nickelates lies in the presence of the Nd-derived bands, whose self-doping
effect on the magnetic and superconducting properties is highly debated, which
is manifested by the opposite arguments supporting the single
Ni-$d_{x^{2}-y^{2}}$ band description oneband ; dft101 ; dft4 ; dft24  versus
the multi-orbital scenarios Hepting , e.g. Ni-$e_{g}$ with or without other
Nd-5d orbitals.

The most significant theoretical reasoning underlying the single-orbital
scenario arises from the demonstration of the depletion of the rare-earth
bands, namely the extremely low electron density of rare-earth layer, which
results in a Fermi surface reconstruction that might render the infinite-layer
nickelates essentially akin to cuprates. Despite there have been strong
theoretical Held2022 ; dft101 ; dft4 ; dft24  and experimental oneband
support on this scenario, it is still not completely resolved whether and how
the low occupation of rare-earth bands impact the charge Pr ; Hall1 ; Hall2 ,
magnetic GuangMing ; GuangMing1 , and superconducting properties Mi2020 ;
dswave ; dft2 ; dft3 ; dft101 ; dft102 ; spinfluc1 ; spinfluc2  in the whole
phase diagram.

##  II Model and Methodology

###  II.1 Two-orbital $d$-$s$ model

To provide more insights on the prominent role of the rare-earth bands, we
employ a minimal two-orbital $d$-$s$ model, where $d$-orbital mimics
Ni-$3d_{x^{2}-y^{2}}$. Additionally, we include the rare-earth R-$5d$
orbitals, which can be effectively approximated as an interstitial s-like
orbital located at the missing apical Oxygen above and below Ni in the rare-
earth layers. The Hamiltonian reads as follows

$\displaystyle\hat{H}=$
$\displaystyle\sum_{k\sigma}(E^{d}_{k}n^{d}_{k\sigma}+E^{s}_{k}n^{s}_{k\sigma})+U\sum_{i}n^{d}_{i\uparrow}n^{d}_{i\downarrow}$
$\displaystyle+\sum_{k\sigma}V_{k}(d^{\dagger}_{k\sigma}s^{\phantom{\dagger}}_{k\sigma}+h.c.)$
(1)

with two orbitals’ dispersion and $d$-$s$ hybridization as

$\displaystyle E^{d}_{k}=$ $\displaystyle-2t_{d}(\cos k_{x}+\cos
k_{y})+4t^{\prime}_{d}\cos k_{x}\cos k_{y}-\mu$ $\displaystyle E^{s}_{k}=$
$\displaystyle-2t_{s}(\cos k_{x}+\cos k_{y})+\epsilon_{s}-\mu$ $\displaystyle
V_{k}=$ $\displaystyle-2V(\cos k_{x}+\cos k_{y})$ (2)

where $d^{\dagger}_{k\sigma}(s^{\dagger}_{k\sigma})$ are creation operators in
momentum space for two orbitals. $n^{d(s)}_{i\sigma}$ and $n^{d(s)}_{k\sigma}$
are the associated number operators in the real and momentum spaces
separately. The chemical potential $\mu$ can be tuned for a desired total
electron density. The on-site energy $\epsilon_{s}$ of $s$ orbital is a
tunable parameter to adjust the relative electron density between two orbitals
with the criterion that the $s$ orbital’s occupancy is around 6-10% in an
undoped system (total electron density $n=1$) in order to be consistent with
the experimental finding that there are approximately 8% electrons per formula
unit for NdNiO2 Hall1 . We adopt the convention that the nearest-neighbor
hopping between $d$ orbitals $t_{d}=1$ is set as the energy unit. Following
the recent numerical investigation TPD_ds , we choose
$t^{\prime}_{d}=-0.25,t_{s}=0.5$ as typical hopping integrals. Besides, the
impact of on-site Coulomb repulsion $U$ of $d$ orbital will be studied in the
present work.

The specific dispersion of the $d$-$s$ hybridization $V_{k}$ might
quantitatively affect the detailed evaluation of various physical properties.
Nevertheless, as pointed out by Ref. TPD_ds , although the symmetry analysis
of the real infinite-layer materials reveals that the inter-layer hoppings
between Ni-3d and the rare-earth R-5d orbitals can be negligibly small, the
electrons in the rare-earth layer are highly itinerant so that the detailed
form of $V_{k}$ probably only has minor effects. For simplicity, therefore, we
focus on the model Hamiltonian where $d$-$s$ has only nearest-neighbor
hoppings dft25 .

One important issue we will explore is the magnitude of $V$, which is believed
to be small $V/t_{d}\sim 0.1$ in previous studies TPD_ds . Strikingly,
however, a recent investigation indicated that the dominant hybridization
between Ni-3d orbitals and itinerant electrons of this interstitial s orbital
comes from a large inter-cell hopping dft25 . Therefore, the present work will
not restrict our attention to only small hybridization $V$ but will study the
influence of its magnitude. Specifically, we choose both $V/t_{d}=0.1$ TPD_ds
and $V/t_{d}=0.6$ estimated from Ref. dft25  as two typical values. The focus
in this work is to investigate the relevance of the $d$-$s$ model to the
infinite-layer nickelates, in particular the most attractive superconducting
instabilities in the hole doped systems.

Before proceeding, we remark on some extensively investigated limits of the
$d$-$s$ model. On the one hand, at $V=0$, namely the two-dimensional single-
band Hubbard model, there have been tremendous theoretical effort for its
solutions and more importantly its relevance in strongly correlated quantum
materials like cuprate and nickelate superconductors. Regarding its
superconducting properties, the consensus is that the $d$-wave pairing is
appreciated owing to its retarded nature to minimize the impact of local
repulsion $U$ ScalapinoRMP ; Maier06 . On the other hand, at finite
hybridization $V$ but $t_{d}=t^{\prime}_{d}=0$, the model recovers the
conventional periodic Anderson model (PAM), which is widely treated as the
standard model of describing the heavy fermion physics. The deeper connection
and moreover the transition or crossover between PAM and $d$-$s$ model will be
left for another future independent study.

###  II.2 Dynamical cluster approximation

We adopt the dynamical cluster approximation (DCA) Hettler98 ; Maier05 ; code
with a continuous time auxilary field (CT-AUX) quantum Monte Carlo (QMC)
cluster solver GullCTAUX  to numerically solve the $d$-$s$ model Eq. (II.1).
As one of the successful many-body methods, DCA has provided much insight on
the strongly correlated physics especially for the Hubbard-type models
Hettler98 ; Maier05 . In principle, DCA evaluates the physical properties in
the thermodynamic limit via mapping the bulk lattice problem onto a finite
cluster embedded in a mean-field self-consistently, where the effective
cluster problem is exactly solved by various perturbative or non-perturbative
many-body methods Hettler98 ; Maier05 . In practice, DCA is realized by a
self-consistent loop for the convergence between the cluster and coarse-
grained single-particle Green’s function. The truly time-consuming part of the
DCA loop lies in the cluster solver, e.g. CT-AUX adopted here. In other words,
most uncertainties of DCA method arises from the cluster solver method e.g.
quantum Monte Carlo procedure. To simulate a wide range of doping levels and
Hubbard interaction, we stick on the smallest $N_{c}=2\times 2$ DCA cluster to
keep the QMC sign problem manageable GullCTAUX  down to the superconducting
transition temperatures $T_{c}$. In fact, the qualitative trend of our major
results can be verified with more expensive simulations employing larger
clusters such as $N_{c}=8$.

To characterize the superconducting instability, we rely on the leading
eigenvalues of the Bethe-Salpeter equation (BSE) in the eigen-equation form in
the particle-particle channel Maier06 ; Scalapino06

$\displaystyle-\frac{T}{N_{c}}\sum_{K^{\prime}}\Gamma^{pp}(K,K^{\prime})\bar{\chi}_{0}^{pp}(K^{\prime})\phi_{\alpha}(K^{\prime})=\lambda_{\alpha}(T)\phi_{\alpha}(K)$
(3)

where $\Gamma^{pp}(K,K^{\prime})$ denotes the lattice irreducible particle-
particle vertex of the effective cluster problem with combining the cluster
momenta ${\bf K}$ and Matsubara frequencies $\omega_{n}=(2n+1)\pi T$ to
$K=({\bf K},i\omega)$. One of the key DCA assumptions is that the desired
lattice two-particle irreducible vertex $\Gamma$ is approximated by the its
cluster counterpart $\Gamma_{c}$, which can be obtained by

$\displaystyle\chi_{c\sigma\sigma^{\prime}}(q,K,K^{\prime})$
$\displaystyle=\chi^{0}_{c\sigma\sigma^{\prime}}(q,K,K^{\prime})+\chi^{0}_{c\sigma\sigma^{\prime\prime}}(q,K,K^{\prime\prime})$
$\displaystyle\times\Gamma_{c\sigma^{\prime\prime}\sigma^{\prime\prime\prime}}(q,K^{\prime\prime},K^{\prime\prime\prime})\chi_{c\sigma^{\prime\prime\prime}\sigma^{\prime}}(q,K^{\prime\prime\prime},K^{\prime})$
(4)

where the cluster two-particle Green’s function

$\displaystyle\chi_{c\sigma\sigma^{\prime}}(q,K,K^{\prime})$
$\displaystyle=\int^{\beta}_{0}\int^{\beta}_{0}\int^{\beta}_{0}\int^{\beta}_{0}d\tau_{1}d\tau_{2}d\tau_{3}d\tau_{4}$
$\displaystyle\times
e^{i[(\omega_{n}+\nu)\tau_{1}-\omega_{n}\tau_{2}+\omega_{n^{\prime}}\tau_{3}-(\omega_{n^{\prime}}+\nu)\tau_{4}]}$
$\displaystyle\times\langle\mathcal{T}$ $\displaystyle
c^{\dagger}_{K+q,\sigma}(\tau_{1})c^{\phantom{\dagger}}_{K\sigma}(\tau_{2})c^{\dagger}_{K^{\prime}\sigma^{\prime}}(\tau_{3})c^{\phantom{\dagger}}_{K^{\prime}+q,\sigma^{\prime}}(\tau_{4})\rangle$
(5)

with conventional notation $K=(\mathbf{K},i\omega_{n})$,
$K^{\prime}=(\mathbf{K^{\prime}},i\omega_{n^{\prime}})$, $q=(\mathbf{q},i\nu)$
and the time-ordering operator $\mathcal{T}$ can be calculated numerically via
a DCA cluster solver (CT-AUX in our case). Concurrently, the non-interacting
two-particle Green’s function
$\chi^{0}_{c\sigma\sigma^{\prime}}(q,K,K^{\prime})$ is constructed from the
product of a pair of fully dressed single-particle Green’s functions. Note the
usual convention that the summation is to be made for repeated indices. In
this work, we are mostly interested in the even-frequency even-parity (spin
singlet) $d$-wave pairing tendency so that $q=(\mathbf{q},i\nu)=0$ is assumed
Maier06 ; Scalapino06 , despite that there exists odd-frequency pairing
channels at high temperatures that quickly decays with lowering $T$.

In the BSE Eq. (3), the coarse-grained bare particle-particle susceptibility

$\displaystyle\bar{\chi}^{pp}_{0}(K)=\frac{N_{c}}{N}\sum_{k^{\prime}}G(K+k^{\prime})G(-K-k^{\prime})$
(6)

is obtained via the dressed single-particle Green’s function $G(k)\equiv
G({\bf k},i\omega_{n})=[i\omega_{n}+\mu-\varepsilon_{\bf k}-\Sigma({\bf
K},i\omega_{n})]^{-1}$, where $\mathbf{k}$ belongs to the DCA patch
surrounding the cluster momentum $\mathbf{K}$, $\mu$ the chemical potential,
$\varepsilon_{\bf k}$ the dispersion relation, and $\Sigma({\bf
K},i\omega_{n})$ the cluster self-energy.

In practice, we normally choose 16 or more discrete points for both the
positive and negative fermionic Matsubara frequency $\omega_{n}=(2n+1)\pi T$
mesh for measuring the two-particle Green’s functions and irreducible vertices
$\Gamma$. Therefore, the BSE Eq. (3) effectively becomes an eigenvalue problem
of a matrix of size $(32N_{c})\times(32N_{c})$, which can be routinely solved
straightforwardly.

Physically, the magnitude of the eigenvalue $\lambda_{\alpha}(T)$ denotes the
normal state pairing tendency in the superconducting channel. Correspondingly,
the spatial, frequency, and orbital dependence of the eigenvector
$\phi_{\alpha}({\bf K},i\omega_{n})$ can be viewed as the normal state analog
of the superconducting gap function to reflect the structure of the pairing
interaction $\Gamma^{pp}$ Maier06 ; Scalapino06 . The superconducting $T_{c}$
is extracted via the temperature where the leading eigenvalue of Eq. (3)
$\lambda(T_{c})=1$. Throughout this work, it is found that the leading pairing
symmetry occurs for the $d$-wave channel with momentum structure $\cos
K_{x}-\cos K_{y}$ so that we are only concerned in the leading eigenvalues
$\lambda_{d}$ and associated $\phi_{d}({\bf K},i\omega_{n})$. In other words,
the presence of $s$ orbital and its hybridization with $d$ orbital does not
affect the $d$-wave pairing symmetry of the conventional Hubbard model with
solely $d$ orbital. Its implication on the infinite-layer nickelates including
the recent experimental findings will be discussed later on.

##  III Results

Figure 1: (a) Temperature dependence of the leading ($d_{x^{2}-y^{2}}$-wave)
eigenvalue $\lambda_{d}(T)$; (b) Doping dependence of the $d$-wave
superconducting $T_{c}$ extracted from $\lambda_{d}(T_{c})=1$; (c) Orbital
occupancy versus the doping level at $T/t_{d}=0.08$. The parameters are
$U/t_{d}=7,N_{c}=4$.

Fig. 1 summarizes our key results. The panel (a) illustrates the temperature
dependence of the leading $d$-wave eigenvalue $\lambda_{d}(T)$ for different
$V$ at two fixed total electron densities $n=0.9,0.85$, namely the total hole
doping $\delta=0.1,0.15$ respectively. What one sees is that the presence of
the hybridization suppresses $d$-wave pairing regardless of the doping level,
which is expected as the additional hybridization tends to affect the original
hopping integral between $d$-orbital so that the effective $U/t_{d}$ will be
altered. In fact, even small $V/t_{d}=0.1$ strongly suppresses
$\lambda_{d}(T)$ at high temperatures while this suppression is gradually
diminished at lower $T$. Interestingly, at high $T$, smaller $V/t_{d}=0.1$ has
stronger effects than larger $V/t_{d}=0.6$ for both $n=0.9,0.85$. Below a
crossover temperature scale $T/t_{d}\sim 0.2$ this difference is reversed.
Although our focus here is the low temperature superconducting instability, we
remark that this stronger suppression due to small $V$ above $T/t_{d}\sim
0.2$, which might be closely related to the promotion of $s$ orbital on SC
that will be discussed later, might be worthwhile for further investigation.

Fig. 1(b) compares the dependence of $T_{c}$ on the global doping level for
the conventional Hubbard model ($V=0$) and $d$-$s$ model at two typical $V$’s.
Apparently, the general feature is that the larger $V$ induces lower $T_{c}$
for all dopings. Nevertheless, this suppression is stronger at low dopings and
it effectively diminishes at high dopings. The impact of $V$ on the variation
of $T_{c}(\delta)$ is closely related to the evolution of $d$ and $s$
orbital’s occupancies.

As displayed in Fig. 1(c), we start from the undoped system ($\delta=0$) with
self-doped $d$ orbital because of finite but tiny $n_{s}\sim 0.08$ to be
consistent with the experiments Hall1 . At small hybridization $V/t_{d}=0.1$,
the doped hole mainly resides on $s$ orbital as evidenced by the rapid
decrease of the dashed orange curve with $\delta$. Recall that experimentally
the superconductivity dome occurs in the range of $0.12<\delta<0.25$. Hence,
in our relevant doping regime, namely $\delta=0.15,0.2$, the $s$ orbital’s
concentration is as tiny as only $\sim 0.01$. In this sense, as pointed out by
previous work TPD_ds ; Held2020 ; Held2022 , the $d$-$s$ model effectively
reduces to the single $d$ orbital system. We numerically confirmed this
single-band picture in the doping range relevant to the realistic materials.
In this case, the decrease of $T_{c}$ shown in Fig. 1(b) originates from the
decrease of $d$ orbital occupancy.

More interestingly, at stronger $V/t_{d}=0.6$, both decrease of $n_{d}$ and
$n_{s}$ with hole doping is smoother than the case of $V/t_{d}=0.1$. In fact,
the doped hole dominantly goes onto $d$ orbital linearly instead, which
implies that the less depleted $s$ orbital may play an important role. From
the viewpoint of the single-orbital Hubbard model, larger hole doping would
rapidly suppress $T_{c}$. Thus, in the regime of $\delta=0.15,0.2$, it is
striking to see that $T_{c}$ at $V/t_{d}=0.6$ conversely approaches to the
curves for $V/t_{d}=0.1$ in Fig. 1(b) in spite of the increasing reduce of
$n_{d}$ for $V/t_{d}=0.6$ compared to $V/t_{d}=0.1$. The only reason
originates from the plausible “enhancing” effect of the small but finite
$n_{s}$ on inhibiting the suppression of SC. This argument also matches with
the elevated $\lambda_{d}(T)$ for $V/t_{d}=0.6$ at high temperatures, namely
the $s$ orbital’s concentration tends to maintain $d$ orbital’s SC. This
supportive role will be more clear when we discuss the $U$-dependence of
$T_{c}$ in Fig. 3. In fact, the $T_{c}(\delta)$ for three $V$’s seemingly
approaches with one another with further hole doping. Unfortunately, due to
the limitation of the CT-QMC sign problem, we are not able to push to even
higher doping levels and also much larger DCA cluster size to confirm this
numerical phenomena to be truly physical.

Putting another way, assuming that the electronic structure calculation dft25
leading to $V/t_{d}\sim 0.6$ is more realistic than the more common belief of
the smaller hybridization strength $V/t_{d}\sim 0.1$, the presence of the
electronic density in rare-earth layer can play the vital role for supporting
the superconductivity. This constitutes one of our major findings in terms of
the constructive role of $d$-$s$ hybridization, which implies that a minimal
model accounting for the infinite-layer nickelate SC might need consider the
role of the rare-earth layer appropriately in some circumstances.

It is reasonable to argue against the supportive role of $s$ orbital on SC
since the $d$-$s$ model has natural connection to PAM, where the $d$-$s$ Kondo
screening has destructive effects on SC. Here we emphasize that $s$’s support
on SC discussed above might only apply for the systems with tiny $s$ electron
density so that the Kondo screening does not have significant impact. We also
note that this extreme limit of PAM has rarely been studied before so that it
is requisite to explore it in more details in future.

Figure 2: The Matsubara frequency dependence of the leading $d$-wave
eigenfunction $\phi_{d}({\bf K},i\omega_{n})/\phi_{d}({\bf K},i\pi T)$ with
${\bf K}=(\pi,0)$ and the normalized cluster spin susceptibility
$2\chi_{s}(\mathbf{Q},i\omega_{m})/[\chi_{s}(\mathbf{Q},0)+\chi_{s}(\mathbf{Q},i2\pi
T]$ for ${\bf Q}=(\pi,\pi)$. The parameters are
$U/t_{d}=7,N_{c}=4,T/t_{d}=0.08$ and $n=0.85$.

Now that Fig. 1(b) has established the notion that the impact of $V$ on the
$d$-wave pairing gradually diminishes in the relatively high doping regime of
$\delta=0.15-0.2$, it is natural to further explore the behavior of the
pairing interaction and moreover the applicability of the spin-fluctuation
mediation in $d$-$s$ model. As mentioned earlier, the structure of pairing
interaction $\Gamma^{pp}$ is reflected by the BSE eigenvector
$\phi_{\alpha}({\bf K},i\omega_{n})$, which is the normal state analog of the
superconducting gap function Maier06 ; Scalapino06 ; Mi2022 .

Fig. 2 illustrates the Matsubara frequency dependence of the normalized
$\phi_{d}({\bf K},i\omega_{n})$ with ${\bf K}=(\pi,0)$. In addition, the
normalized cluster spin susceptibility $\chi_{s}(\mathbf{Q},i\omega_{m})$ for
${\bf Q}=(\pi,\pi)$, which is calculated via
$\sum_{K,K^{\prime}}\chi_{c}(q,K,K^{\prime})$ in Eq. (II.2), is given for
comparison to check the role of spin-fluctuation in mediating the $d$-wave
pairing as done in conventional Hubbard model Maier06 ; Scalapino06 . The
decrease of $\phi_{d}$ versus $i\omega_{n}$ has a characteristic energy scale,
which looks increasing with $V$ of a tiny amount. This indicates that the
$d$-wave pairing is weakened slightly at larger $V$, which is consistent with
$T_{c}$’s decrease in Fig. 1(b). Regarding the antiferromagnetic spin
susceptibility, it has been normalized to coincide with
$\phi_{d}(i\omega_{n})$ at $\omega_{n}=\pi T$. Note that the Matsubara
frequency entering $\chi_{s}$ is the bosonic $\omega_{m}=2m\pi T$. It can be
seen that the antiferromagnetic spin-fluctuation spectrum is characterized by
the same energy scale as the dynamics of $\phi_{d}(i\omega_{n})$ reflecting
the pairing interaction for $V/t_{d}=0.1$. Nevertheless,
$\chi_{s}(i\omega_{m})$ almost overlap between two cases of $V/t_{d}=0.1,0.6$
and thereby the spin-fluctuation and pairing interaction have tiny discrepancy
for $V/t_{d}=0.6$, which might be related to the non-negligible effects of $s$
orbital discuss earlier. Fig. 2 implies that the effectively attractive
interaction mediated by the antiferromagnetic fluctuations and the retardation
nature of $d$-wave pairing Scalapino06  are not qualitatively modified in the
presence of $s$ orbital.

Figure 3: $T_{c}$ versus $U$ for $n=0.85,N_{c}=4$. At stronger $V/t_{d}=0.6$,
the “optimal” $U$ for highest $T_{c}$ is shifted to a higher value.

To further characterize the SC properties of $d$-$s$ model, Fig. 3 shows the
evolution of $T_{c}$ versus $U$ for different hybridization $V$ at a fixed
total density $n=0.85$. Apparently, the small $V$ induces a constant downshift
of $T_{c}$ for all $U$’s, indicating that the extremely tiny concentration of
$s$ orbital (shown in Fig. 1(c)) do not qualitatively modify the
superconducting instability. In particular, the existence of an “optimal”
$U/t_{d}\sim 6$ of supporting the highest $T_{c}$ occurs for both single-
orbital Hubbard model $V=0$ and $d$-$s$ model at $V/t_{d}=0.1$. This peak
structure is consistent with the notion that it is important to have strong
short-range antiferromagnetic correlations, which, however, decreases when $U$
becomes sufficiently large compared to the bandwidth. At strong interaction
$U$, the exchange interaction $J\sim 4t^{2}/U$ decays so that $T_{c}$
decreases accordingly.

Distinct from the weak hybridization, at relatively strong $V/t_{d}=0.6$, the
alleviated QMC sign problem allows us to simulate the stronger $U/t_{d}=10$ so
that the similar peak structure of $T_{c}(U)$ can be seen, namely the
“optimal” $U$ is shifted to a higher value compared to the case of
$V/t_{d}=0.1$. Although we are not able to numerically confirm the trend at
even higher $U$ for small $V$, Fig. 3 implies that the additional $s$ orbital
and its relatively strong hybridization $V$ with $d$ orbital can have
potentially supportive role for an elevated $T_{c}$ at strong interaction,
which further backup our previous arguments discussed for Fig. 1. For
instance, at $U/t_{d}=9$, $T_{c}(V/t_{d}=0.6)>T_{c}(V/t_{d}=0.1)$.
Unfortunately, it is quite challenging, if not impossible, to determine
whether $T_{c}$ can even exceed the Hubbard model at larger $U$ with the aid
of the additional $s$ orbital with appropriate hybridization $V$.

##  IV Discussion

Here we rephrase and remark on some insights on the superconducting properties
of the infinite-layer nickelates and more generally other relevant correlated
materials.

$T_{c}$’s magnitude:

Fig. 1(b) reveals the pair-breaking effects of $d$-$s$ hybridization,
especially at relatively strong $V/t_{d}=0.6$. Considering that the cuprate
superconductors’ optimal hole doping occurs around $\delta\sim 12.5\%$ while
the infinite-layer nickelates at higher $\delta\sim 18\%$, $T_{c}$ has a
significant drop regardless of the hybridization strength. This observation is
generally consistent with the experimental fact that $T_{c}$ of infinite-layer
nickelates is much lower than the cuprates.

Single-orbital picture:

The applicability of the single-orbital description depends on the
hybridization $V$ strength. Fig. 1 indicates that if $V$ is indeed
sufficiently small, it only results in a negligible decrease of $T_{c}$
corresponding to the tiny $s$ orbital concentration.
Nevertheless, if the interlayer hybridization between NiO2 and R layers turns
out not to be negligible, a minimal model of the infinite-layer nickelate SC
or other similar compounds might need considering the role of the rare-earth
layer appropriately even though the R layers have only tiny electron
densities.

Undoped system:

At first glance, the curves in Fig. 1(b) might hint towards a finite $T_{c}$
of undoped systems. This obvious discrepancy with the realistic observations
probably arises from the small DCA cluster size. Previous expensive DCA
calculations with larger $N_{c}$ but with smaller $U$ has uncovered the dome-
shaped $T_{c}$ versus hole doping Jarrell2013 . In this sense, our study here
is simply proof-of-principle in terms of the effects of the hybridization $V$.

$V$’s promotion of $T_{c}$:

As discussed in Fig. 1 and Fig. 3, it is striking to observe that the tiny but
finite $s$ orbital concentration has promotion effects on SC so that the
hybridization strength between NiO2 and rare-earth layers might require more
attention, for instance, whether the electronic structure calculation dft25
is more realistic. Even though this picture would be ultimately proved to be
irrelevant to the infinite-layer nickelates, it might provide some insights on
other strongly correlated compounds.

Optimal $U$:

Our numerical calculations uncovered some interesting trend of elevated
$T_{c}$ due to relatively large hybridization $V$. Together with the previous
point, it might provide some hints for the material design for compounds (not
only infinite-layer nickelates) with higher $T_{c}$.

Heavy fermion physics:

Despite that our current focus is only the superconducting properties at large
hole dopings, the potential heavy fermion physics at or close to the undoped
systems deserves further examination. The $d$-$s$ model has natural relation
to the Anderson lattice model except for the extremely low conduction electron
density. It is worthwhile exploring the potential Kondo screening exhaustion
exhaustion2019  and its impact on the magnetic properties of the undoped
infinite-layer nickelates.

Pairing symmetry:

Our present work follows the established framework of $d$-wave symmetries
based pairing mechanism TPD_ds ; Held2022 ; Held2020 ; Millis2022 ;
Xianxin2020 ; Adhikary2020 ; Matianxing2020 . However, the most recent London
penetration measurement strongly challenges this scenario by revealing the
predominantly nodeless pairing nodeless2022 . Specifically, there might exist
a 2D-to-3D superconducting states crossover possibly due to the coupling
between NiO2 plane and rare-earth spacer layer. This similarity with the iron-
based superconductor with nodeless multiband SC is reminiscent of the various
studies on the bilayer Hubbard model at strong inter-layer hybridization,
where the nodeless $s^{\pm}$-wave pairing is appreciated bilayer2011 ;
bilayer2016 ; bilayer2021 . In this regard, therefore, the $d$-$s$ model as
the minimal model seems problematic at least in the regime of relatively small
inter-orbital hybridization $V$. This direction may attractive further
attention and exploration.

##  V Summary

In conclusion, we investigated a two-orbital $d$-$s$ model to mimic the
interplay between Ni-$d_{x^{2}-y^{2}}$ and the effective interstitial s-like
orbital located at the missing apical Oxygen above and below Ni in the rare-
earth layers to account for the rare-earth 5d orbitals. To accurately approach
to the superconducting $T_{c}$, we employed the dynamic cluster quantum Monte
Carlo focussing on the impact of the $d$-$s$ hybridization.

It is found that the $d$-$s$ hybridization suppresses the $d$-wave pairing
regardless of the doping level. The suppression of $T_{c}$ is stronger at low
dopings and diminishes with further dopings, which is closely related to the
evolution of $d$ and $s$ orbital occupancies. In particular, although the
doped hole mainly resides on $s$ orbital at small hybridization $V/t_{d}=0.1$
so that the $d$-$s$ model effectively reduces to the single $d$ orbital
system; at stronger $V/t_{d}=0.6$, $d$ orbital has considerable concentration
even at high doping levels, whose role is manifested by the weakened
suppression trend of $V$ on SC. In other words, $T_{c}$ at $V/t_{d}=0.6$
approaches to that at $V/t_{d}=0.1$ at high dopings. It is thus plausible to
conjecture about a striking promotion effect of the small but finite $s$
orbital occupancy on SC.

To explore how to promote $T_{c}$, we further study the impact of Hubbard $U$.
Although the weak hybridization $V$ does not qualitatively modify the
superconducting properties without $s$ orbital, the relatively large
$V/t_{d}=0.6$ induces that the “optimal” $U$ is shifted to a higher value
compared to the case of $V/t_{d}=0.1$. This might imply a potential strategy
to enhance the SC instability in general strongly correlated systems including
the infinite-layer nickelates via appropriate hybridization with additional
conduction electrons.

In addition, in spite of the additional $s$ orbital, the comparison between
the Matsubara frequency dependence of the $d$-wave eigenvector (reflecting the
pairing interaction) and the antiferromagnetic spin susceptibility supports
the spin fluctuation mediation in $d$-wave pairing of $d$ orbital.

Our presented work provides complemental insights on the quantitative effects
of the hybridization between NiO2 and rare-earth layers on the superconducting
instability relevant to the infinite-layer nickelate superconductors,
especially the potentially enhancing effects of rare-earth layers on SC. As a
final remark, we believe that the necessity of including the rare-earth layer
in a minimal model of infinite-layer nickelates and other relevant strongly
correlated compounds can be crucial in the presence of a considerable
interlayer hybridization.

## Acknowledgments

We acknowledge fruitful discussion with Zhicheng Zhong and Guangming Zhang.
This work was supported by National Natural Science Foundation of China (NSFC)
Grant No. 12174278, the startup fund from Soochow University, and Priority
Academic Program Development (PAPD) of Jiangsu Higher Education Institutions.

## References

  * (1) D. Li, K. Lee, B. Y. Wang, M. Osada, S. Crossley, H. R. Lee, Y. Cui, Y. Hikita, and H. Y. Hwang, Nature (London) 572, 624 (2019).
  * (2) J. Zhang and X. Tao, CrystEngComm,2021,23,3249 (2021).
  * (3) A. S. Botana, F. Bernardini, and A. Cano, arXiv: 2012.02764 (2020).
  * (4) M. Osada, B. Y. Wang, K. Lee, D. Li, and H. Y. Hwang, Phys. Rev. Materials 4, 121801(R) (2020).
  * (5) S. W. Zeng, C. J. Li, L. E. Chow, Y. Cao, Z. T. Zhang, C. S. Tang, X. M. Yin, Z. S. Lim, J. X. Hu, P. Yang, A. Ariando, arXiv:2105.13492 (2021).
  * (6) G. A. Pan, D. F. Segedin, H. LaBollita, Q. Song, E. M. Nica et al, Nat. Mater. (2021). https://doi.org/10.1038/s41563-021-01142-9
  * (7) K.-W. Lee and W. E. Pickett, Phys. Rev. B 70, 165109 (2004).
  * (8) M. Hepting, D. Li, C. J. Jia, H. Lu, E. Paris, Y. Tseng, X. Feng, M. Osada, E. Been, Y. Hikita, Y. D. Chuang, Z. Hussain, K. J. Zhou, A. Nag, M. Garcia-Fernandez, M. Rossi, H. Y. Huang, D. J. Huang, Z. X. Shen, T. Schmitt, H. Y. Hwang, B. Moritz, J. Zaanen, T. P. Devereaux, and W. S. Lee, Nature Mater. 19, 381 (2020).
  * (9) A. S. Botana and M. R. Norman, Phys. Rev. X 10, 011024 (2020).
  * (10) J. Kapeghian and A. S. Botana, Phys. Rev. B 102, 205130 (2020).
  * (11) Z. Liu, Z. Ren, W. Zhu, Z. Wang, and J. Yang, npj Quantum Materials 5, 31 (2020).
  * (12) I. Leonov, S. L. Skornyakov, and S. Y. Savrasov, Phys. Rev. B 101, 241108 (2020).
  * (13) I. Leonov and S. Y. Savrasov, arXiv:2006.05295.
  * (14) Y. Gu, S. Zhu, X. Wang, J. Hu, and H. Chen, Communications Physics 3, 84 (2020).
  * (15) A. Ikeda, Y. Krockenberger, H. Irie, M. Naito, and ?H. Yamamoto, Applied Physics Express 9, 061101 (2016); M. Hayward, M. Green, M. Rosseinsky, and J. Sloan, J. Am. Chem. Soc. 121, 8843 (1999).
  * (16) H. Lu, M. Rossi, A. Nag, M. Osada, D. F. Li, K. Lee, B. Y. Wang, M. Garcia-Fernandez, S. Agrestini, Z. X. Shen, E. M. Been, B. Moritz, T. P. Devereaux1, J. Zaanen, H. Y. Hwang, Ke-Jin Zhou, W. S. Lee, Science 373, 6551 (2021).
  * (17) J. Q. Lin, P. Villar Arribi, G. Fabbris, A. S. Botana, D. Meyers, et al, Phys. Rev. Lett. 126, 087001 (2021).
  * (18) Y. Cui, C. Li, Q. Li, X. Zhu, Z. Hu, Y.-f. Yang, J. S. Zhang, R. Yu, H.-H. Wen, and W. Yu, Chin. Phys. Lett. 38, 067401 (2021).
  * (19) R. A. Ortiz, P. Puphal, M. Klett, F. Hotz, R. K. Kremer, H. Trepka, M. Hemmida, H.-A. Krug von Nidda, M. Isobe, R. Khasanov, H. Luetkens, P. Hansmann, B. Keimer, T. Schafer, and M. Hepting, arXiv:2111.13668 (2021).
  * (20) J. Fowlie, M. Hadjimichael, M. M. Martins, D. Li, M. Osada, B. Y. Wang, K. Lee, Y. Lee, Z. Salman, T. Prokscha, J.-M. Triscone, H. Y. Hwang, A. Suter, arXiv:2201.11943 (2022).
  * (21) H. Sakakibara, H. Usui, K. Suzuki, T. Kotani, H. Aoki, and K. Kuroki, Phys. Rev. Lett. 125, 077003 (2020).
  * (22) X. Wu, D. Di Sante, T. Schwemmer, W. Hanke, H. Y. Hwang, S. Raghu, and R. Thomale, Phys. Rev. B 101, 060504 (2020).
  * (23) Y. Nomura, M. Hirayama, T. Tadano, Y. Yoshimoto, K. Nakamura, and R. Arita, Phys. Rev. B 100, 205138 (2019).
  * (24) J. Gao, S. Peng, Z. Wang, C. Fang, and H. Weng, National Science Review (2020).
  * (25) P. Jiang, L. Si, Z. Liao, and Z. Zhong, Phys. Rev. B 100, 201106 (2019).
  * (26) E. Been, W.-S. Lee, H. Y. Hwang, Y. Cui, J. Zaanen, T. Devereaux, B. Moritz, and C. Jia, arXiv:2002.12300.
  * (27) F. Bernardini, V. Olevano, and A. Cano, Phys. Rev. Research 2, 013219 (2020).
  * (28) J. Krishna, H. LaBollita, A. O. Fumega, V. Pardo, and A. S. Botana, Phys. Rev. B 102, 224506 (2020).
  * (29) H. Zhang, L. Jin, S. Wang, B. Xi, X. Shi, F. Ye, and J.-W. Mei, Phys. Rev. Research 2, 013214 (2020).
  * (30) M. Kitatani, L. Si, O. Janson, R. Arita, Z. Zhong, and K. Held, npj Quantum Materials 5, 59 (2020).
  * (31) K. Held, L. Si, P. Worm, O. Janson, R. Arita, Z. Zhong, J. M. Tomczak, M. Kitatani, arXiv:2201.0122 (2022).
  * (32) H. Chen, A. Hampel, J. Karp, F. Lechermann, and A. Millis, arXiv:2201.02852 (2022).
  * (33) M. Kitatani, L. Si, O. Janson, R. Arita, Z. Zhong, and K. Held, npj Quantum Materials 5, 59 (2020).
  * (34) P. Werner and S. Hoshino, Phys. Rev. B 101, 041104 (2020).
  * (35) V. Olevano, F. Bernardini, X. Blase, and A. Cano, Phys. Rev. B 101, 161102 (2020)
  * (36) F. Petocchi, V. Christiansson, F. Nilsson, F. Aryasetiawan, and P. Werner, Phys. Rev. X 10, 041047 (2020).
  * (37) S. Ryee, H. Yoon, T. J. Kim, M. Y. Jeong, and M. J. Han, Phys. Rev. B 101, 064513 (2020).
  * (38) F. Lechermann, Phys. Rev. B 101, 081110 (2020).
  * (39) J. Karp, A. S. Botana, M. R. Norman, H. Park, M. Zingl, and A. Millis, Phys. Rev. X 10, 021061 (2020).
  * (40) F. Lechermann, Phys. Rev. X 10, 041002 (2020).
  * (41) X. Wan, V. Ivanov, G. Resta, I. Leonov, and S. Y. Savrasov, arXiv:2008.07465.
  * (42) Y. Wang, C.-J. Kang, H. Miao, and G. Kotliar, Phys. Rev. B 102, 161118 (2020).
  * (43) C.-J. Kang and G. Kotliar, arXiv:2007.15383.
  * (44) Z.-J. Lang, R. Jiang, and W. Ku, arXiv:2005.00022.
  * (45) J. Karp, A. Hampel, M. Zingl, A. S. Botana, H. Park, M. R. Norman, and A. J. Millis, Phys. Rev. B 102, 245130 (2020).
  * (46) M. Rossi, H. Lu, A. Nag, D. Li, M. Osada, K. Lee, B. Y. Wang, S. Agrestini, M. Garcia-Fernand ez, Y. D. Chuang, Z. X. Shen, H. Y. Hwang, B. Moritz, K.-J. Zhou, T. P. Devereaux, and W. S. Lee, arXiv:2011.00595.
  * (47) D. Li, B. Y. Wang, K. Lee, S. P. Harvey, M. Osada, B. H. Goodge, L. F. Kourkoutis, and H. Y. Hwang, Phys. Rev. Lett. 125, 027001 (2020).
  * (48) S. Zeng, C. S. Tang, X. Yin, C. Li, M. Li, Z. Huang, J. Hu, W. Liu, G. J. Omar, H. Jani, Z. S. Lim, K. Han, D. Wan, P. Yang, S. J. Pennycook, A. T. S. Wee, and A. Ariando, Phys. Rev. Lett. 125, 147003 (2020).
  * (49) G.-M. Zhang, Y.-F. Yang, and F.-C. Zhang, Phys. Rev. B 101, 020501(R) (2020).
  * (50) Z. Wang, G.-M. Zhang, Y.-F. Yang, and F-C Zhang, Phys. Rev. B102, 220501(R) (2020).
  * (51) Q. Gu, Y. Li, S. Wan, H. Li, W. Guo, H. Yang, Q. Li, X. Zhu, X. Pan, Y. Nie, and H.-H. Wen, Nat. Comm. 11, 6027 (2020).
  * (52) T. Zhou, Y. Gao, and Z. Wang, Science China Physics, Mechanics, and Astronomy 63, 287412 (2020).
  * (53) P. Adhikary, S. Bandyopadhyay, T. Das, I. Dasgupta, and T. Saha-Dasgupta, Phys. Rev. B 102, 100501 (2020).
  * (54) M. Jiang, M. Berciu, and G. A. Sawatzky, Phys. Rev. Lett. 124, 207004 (2020).
  * (55) C. Peng, H.-C. Jiang, B. Moritz, T. P. Devereaux, and C. Jia, arXiv:2110.07593 (2021).
  * (56) D. J. Scalapino, Rev. Mod. Phys. 84 1383 (2012).
  * (57) T. A. Maier, M. S. Jarrell, D. J. Scalapino, Phys. Rev. Lett. 96, 047005 (2006).
  * (58) T. Maier, M. Jarrell, T. Pruschke, and M. Hettler, Rev. Mod. Phys. 77, 1027 (2005).
  * (59) U. R. Hahner, G. Alvarez, T. A. Maier, R. Solca, P. Staar, M. S. Summers, and T. C. Schulthess, Comput. Phys. Commun. 246, 106709 (2020). The DCA++ code used for this project can be obtained at https://github.com/CompFUSE/DCA.
  * (60) M. H. Hettler, A. N. Tahvildar-Zadeh, M. Jarrell, T. Pruschke, and H. R. Krishnamurthy, Phys. Rev. B 58, R7475 (1998).
  * (61) E. Gull, P. Werner, O. Parcollet, M. Troyer, Europhys. Lett. 82, 57003 (2008).
  * (62) D. J. Scalapino, Chapter 13 in the "Handbook of High Temperature Superconductivity", J. R. Schrieffer, editor, Springer (2006).
  * (63) M. Jiang, Phys. Rev. B 105, 024510 (2022).
  * (64) K.-S. Chen, Z. Y. Meng, S.-X. Yang, T. Pruschke, J. Moreno, and M. Jarrell, Phys. Rev. B 88, 245110 (2013).
  * (65) L. Zhang, T. Ma, N. C. Costa, R. R. dos Santos, and R. T. Scalettar, Phys. Rev. B 99, 195147 (2019).
  * (66) J. Karp, A. Hampel, and A. J. Millis, arXiv:2201.10481 (2022).
  * (67) X. Wu, D. D. Sante, T.Schwemmer, W. Hanke, H. Y. Hwang, S. Raghu, and R. Thomale, Phys. Rev. B 101, 060504(R) (2020).
  * (68) P. Adhikary, S. Bandyopadhyay, T. Das, I. Dasgupta, and T. Saha-Dasgupta, Phys. Rev. B 102, 100501(R) (2020)
  * (69) C. Chen, R. Ma, Xuelei Sui, Y. Liang, B. Huang, T. Ma, arXiv:2201.06950 (2022).
  * (70) L. E. Chow, S. K. Sudheesh, P. Nandi, S. W. Zeng, Z. T. Zhang, X. M. Du, Z. S. Lim, Elbert E. M. Chia, A. Ariando, arXiv:2201.10038 (2022).
  * (71) T. A. Maier and D. J. Scalapino, Phys. Rev. B 84, 180513(R) (2011).
  * (72) V. Mishra, D. J. Scalapino, and T. A. Maier, Sci. Rep. 6, 32078 (2016).
  * (73) S. Karakuzu, S. Johnston, and T. A. Maier, arXiv:2107.13996 (2021).


From finite nuclei to neutron stars: the essential role of high-order density
dependence in effective forces ††footnotetext: Supported by the National Key
R&D Program of China (Grant No. 2018YFA0404403), the National Natural Science
Foundation of China (Grant Nos. 11975032, 11835001, 11790325, 11961141003.).
∗Corresponding author. Email<EMAIL_ADDRESS>© 2020 Chinese Physical Society
and IOP Publishing Ltd
C.J. Jiang, Y. Qiang, D.W. Guan, Q.Z. Chai, C.Y. Qiao, and J.C. Pei∗

State Key Laboratory of Nuclear Physics and Technology, School of Physics,
Peking University, Beijing 100871, China
(Received 28 December 2020; accepted xxx; published online )

A unified description of finite nuclei and equation of state of neutron stars
present a major challenge as well as opportunities for understandings of
nuclear interactions. Inspired by the Lee-Huang-Yang formula of hard-sphere
gases, we developed effective nuclear interactions with an additional high-
order density dependent term. The original Skyrme force SLy4 is widely used in
studies of neutron stars but is not satisfied for global descriptions of
finite nuclei. The refitted SLy4′ force can improve descriptions of finite
nuclei but slightly reduces the radius of neutron star of 1.4 solar mass. We
found that the extended SLy4 force with a higher-order density dependence can
properly describe properties of both finite nuclei and GW170817 binary neutron
stars, including the mass-radius relation and the tidal deformability. This
demonstrated the essential role of high-order density dependence at ultrahigh
densities. Our work provides a unified and predictive model for neutron stars,
as well as new insights for the future development of effective interactions.

PACS: 21.65.+f; 26.60.+c; 21.30.-x; 21.10.Dr DOI:
10.1088/0256-307X/37/x/xxxxxx

_Introduction._ — Neutron stars (NSs), which are compact objects born in the
supernova explosion,[1] attract nuclear physicists’ strong interests as they
are a unique laboratory for studying the properties of ultradense nuclear
matter.[2,3] The astrophysical observables of NS’ mass, radius and the mass-
radius relation can provide constraints on nuclear theories extrapolated from
properties of finite nuclei.[4] On 17 August 2017, the gravitational waves
emitted from the merging binary NS GW170817[5] are detected by the Advanced
LIGO and the Advanced Virgo, marking a new era of multi-messenger astronomy.
The GW170817 event not only extends our knowledge of NSs’ mass and radius, but
also provides a new astronomical observable, the quadrupolar tidal
deformability[6], as a further constraint.[7-9]

The theoretical descriptions of NS’s mass, radius and tidal deformability
depend solely on the equation of state (EOS) of nuclear matter.[10-13]
Therefore, to seek a unified description of properties of finite nuclei and
EoS of neutron stars is a major challenge. Because Quantum Chromodynamics
(QCD) is non-perturbative at low energies, it is in principle possible to
obtain realistic model-independent nuclear forces from Lattice QCD
calculations.[14] However, such kind of precise calculations are so formidable
that can be accomplished by today’s supercomputers. The chiral effective field
theory based on QCD can derive accurate realistic forces for light nuclei but
can propagate large uncertainties towards heavy nuclei and dense nuclear
matter, due to the less constrained many-body forces.[15] To this end, various
effective nuclear interactions are extensively explored for EOS of neutron
stars.[16,17] The effective interactions are developed by taking into account
the medium effects, i.e., dependent on surrounding nuclear densities, in
contrast to realistic forces in free spaces. The effective interactions are
rather well constrained around the saturation density with experimental
studies of finite nuclei. However, the extrapolation of effective interactions
is very risky to neutron stars which could be several times denser than the
saturation density. There are experiments to study the correlation between
finite nuclei and neutron stars by measurements of neutron skins of neutron-
rich nuclei, combined with model dependent analysis.[18,19] It is expected
that improved measurements of symmetry energies via heavy ion collisions can
directly study EOS at high densities.[20]

In this Letter, we concern about the applicability of the widely used Skyrme
forces for neutron stars. The Skyrme force is a low-momentum effective forces
which is very successful in descriptions of structures and dynamics of the
whole nuclear chart.[21,22] Among many Skyrme forces, the SLy4
parameterization [23,24] is the very best choice for descriptions of neutron
stars and has been widely used as a reference in astrophysical
calculations.[6] The development of SLy4 has taken into account ab initio
UV14+UVII results of EoS of pure neutron matter at high densities.[23] Many
effective forces with a soft EOS have been excluded by observations of neutron
stars around 2 solar mass.[25] The inclusion of hyperons in NS core is
questionable due to its soft EOS.[16,17] For example, the best descriptions of
binding energies of finite nuclei is obtained by the UNEDF0 parameterization
with a global RMSE of 1.5 MeV, however, which is failed to properly describe
neutron stars.[26] On the other hand, SLy4 force has a global RMSE of 4.37 MeV
which is too large.[27]

In Skyrme forces, it is crucial to include a repulsive density dependent term
to describe the nuclear saturation properties and pressure.[28] The density
dependent term also simulates the effects of 3-body and many-body forces
beyond the low-momentum two-body interactions. It is understood that the
induced 3-body and many-body forces are particularly significant in low-
momentum (or soft) interactions, which would be less significant in realistic
(or hard) interactions with a high momentum cutoff.[29] In this case, the
density dependent term in soft Skyrme forces should be carefully treated. We
speculate that it must be insufficient for Skyrme forces with a simple density
dependent term to describe a variety of nuclear systems from dilute nuclear
halos to dense neutron stars. For systems of hard spheres, its exact energy
density functional was derived by Lee, Huang, Yang in a series of
publications.[30,31] The energy density functional form is written as[32],

$\begin{array}[]{ll}\varepsilon=&\displaystyle\frac{3\hbar^{2}}{10m}(3\pi^{2})^{2/3}\rho^{5/3}+\frac{\hbar^{2}\pi
a}{m}\rho^{2}\vspace{5pt}\\\
&\displaystyle+\frac{2\hbar^{2}a^{2}3^{4/3}\pi^{2/3}}{35m}(11-2\ln
2)\rho^{7/3}\vspace{5pt}\\\
&\displaystyle+0.78\frac{\hbar^{2}a^{3}3^{5/3}\pi^{10/3}}{2m}\rho^{8/3}\\\
\end{array}$ (1)

in which $a$ denotes a finite radius of particles. The higher order density
dependent terms appear naturally as a result of a finite radius of particles.
The higher-order correction certainty impacts the equation of state, which has
been proved precisely by experiments of cold atomic gases.[33] The finite size
effects of nucleons is related to repulsive hard cores in realistic forces but
have never been explicitly considered by effective nuclear forces. The aim of
this work is to study the influences of high-order density dependencies in
extended Skyrme forces on observables of neutron stars.

_Extended Skyrme forces._ — The standard Skyrme force is a low momentum
effective interaction for Skyrme Hartree-Fock calculations and time-dependent
dynamic calculations. It includes two-body terms with the momentum expansion
to the second order and a density dependent term
$\rho^{\gamma}(\vec{r}_{1}-\vec{r}_{2})$. The density dependent term simulates
the three-body and many body correlations. However, the power of density
dependence $\gamma$ is largely uncertain. Inspired by the Lee-Huang-Yang
formula, we extended the Skyrme force with an additional higher-order density
dependent term. The energy per nucleon $\frac{E}{A}(\rho)$ in nuclear matter
in terms of density $\rho$ is written as,[27]

$\displaystyle{\frac{E}{A}(\rho)=}$ | $\displaystyle{{}\frac{3\hbar^{2}}{10m}\left(3\pi^{2}\right)^{2/3}\rho^{2/3}+\frac{1}{4}t_{0}(1-x_{0})\rho}$
---|---
|
$\displaystyle{{}+\frac{3}{40}[t_{1}(1-x_{1})+3t_{2}(1+x_{2})]\left(3\pi^{2}\right)^{2/3}\rho^{5/3}}$
|
$\displaystyle{{}+\frac{1}{24}\sum\limits_{i}t_{3i}(1-x_{3i})\rho^{\gamma_{i}+1}}$
(2)

Note that Eq.(2) is for the pure neutron matter. The energy density functional
of combined protons and neutrons is given in Ref.[23]. In the extended Skyrme
forces, there are two density density dependent terms associated with
$t_{3i}(i=1,2)$ and $\gamma_{i}$. We take $\gamma_{1}$=1/6 and
$\gamma_{2}$=1/2 based on SLy4 forces.[27] This can be seen as a reasonable
and perturbative improvement and it is different from the series of expansions
of $\rho^{n/3}$ terms as in Ref.[34]. The pressure is given by
$P(\rho)=\rho^{2}\frac{\partial\frac{E}{A}(\rho)}{\partial\rho}$. The symmetry
energy is given by the energy difference between pure neuron matter and
symmetric nuclear matter.

For simplicity, we only refit the momentum-independent parameters
$t_{0},t_{3i},x_{0},x_{3i}$. These momentum-independent parameters are
directly correlated in the $s$-wave channel. Other parameters are kept as the
same as SLy4 forces. We obtained different parameter sets by fitting binding
energies and charge radii of groups of light nuclei, heavy nuclei, and global
nuclei. We also refitted the SLy4 force with the group of global nuclei. The
details of the optimizations and lists of parameterizations are given in
Ref.[27]

_Equilibrium of $npe\mu$ matter._— Nuclear matter in neutron star is expected
to be a mixture of baryons (protons and neutrons) and leptons in equilibrium
with respect to weak interactions [13,23]. The relative fractions of $npe\mu$
components are crucial for the rate of neutrino cooling of neutron stars. The
total energy density is a function of particle densities
$\rho_{n},\rho_{p},\rho_{e},\rho_{\mu}$, and is written as:

$\displaystyle{\epsilon(r)}$ | $\displaystyle{{}=\epsilon_{b}(\rho_{n},\rho_{p})+\rho_{n}m_{n}c^{2}+\rho_{p}m_{p}c^{2}+\epsilon_{e}(\rho_{e})}$
---|---
|
$\displaystyle{{}+\epsilon_{\mu}(\rho_{\mu})+\rho_{e}m_{e}c^{2}+\rho_{\mu}m_{\mu}c^{2}}$
(3)

Note that the charge neutrality is kept as $\rho_{p}=\rho_{e}+\rho_{\mu}$. The
equilibrium of the chemical potentials (defined as
$\mu_{j}={\displaystyle{\partial\epsilon\over\partial\rho_{j}}}$) is written
as:

$\displaystyle{\mu_{n}=\mu_{p}+\mu_{e},\quad\mu_{\mu}=\mu_{e}}$
---
(4)

Based on the equilibrium equation, we can calculate the densities of various
components at a given baryon density and thus the equation of state at
$npe\mu$ equilibrium is obtained.

_Tolman-Oppenheimer-Volkov equation._ — The TOV equation describes the static
density distributions of spherical neutron star at gravitational equilibrium,
which is derived from general relativity.[35,36] The TOV differential equation
in terms of pressure $P(r)$, energy density $\epsilon(r)$, and mass $M(r)$ is
written as:

$\displaystyle{\frac{\mathrm{d}P}{\mathrm{d}r}=-\frac{GM\epsilon}{r^{2}c^{2}}(1+\frac{P}{\epsilon})(1+\frac{4\pi
r^{3}P}{Mc^{2}})(1-\frac{2GM}{rc^{2}})^{-1}}$
---
(5)

The mass $M(r)$ is the total mass contained inside radius $r$, while $M(0)$=0.

With calculated EOS $P(\rho)$, the TOV equation can be numerically solved from
the center with an initial density to the surface where the pressure becomes
zero. We obtain the density distribution and total mass, and the associated
radius with certain boundary conditions. It can be seen that EoS is crucial to
determine the mass-radius relation of neutron stars.

_Tidal deformability._ —The dimensionless tidal deformability $\Lambda$ is a
functional of dimensionless compactness $\mathcal{C}=GM/Rc^{2}$.[9] $\Lambda$
becomes an important additional constraints on theoretical EOS, since it was
first inferred from gravitational waves of GW170817.[6] The tidal
deformability $\Lambda$ is obtained by solving equations as follows.[37,38]
First, a differential equation is solved to obtain $y(r)$ . $y_{R}$ is the
value of $y(R)$ when $R$ is the NS radius.

$\frac{\mathrm{d}y(r)}{\mathrm{d}r}=-\frac{1}{r}[{y(r)}^{2}+y(r)F(r)+r^{2}Q(r)]$
(6)

The above functions $F(r)$ and $Q(r)$, related to the total mass $M(r)$ and
pressure $P(r)$, energy density $\epsilon(r)$ at radius $r$, is written as,
[9]

$F(r)=\left\\{1-4\pi r^{2}G[\epsilon-P]/c^{4}\right\\}\left(1-{\displaystyle{2GM\over rc^{2}}}\right)^{-1}$ (7) $\displaystyle{Q(r)=}$ | $\displaystyle{{}\left\\{\frac{4\pi G}{c^{4}}\left[5\epsilon+9P+\frac{\epsilon+P}{\partial P/\partial\epsilon}\right]-\frac{6}{r^{2}}\right\\}}$
---|---
|
$\displaystyle{{}\left(1-\frac{2GM}{rc^{2}}\right)^{-1}-\frac{4G^{2}M^{2}}{r^{4}c^{4}}}$
| $\displaystyle{{}\left(1+\frac{4\pi
r^{3}P}{Mc^{2}}\right)^{2}\left(1-\frac{2GM}{rc^{2}}\right)^{-2}}$
(8)

After obtaining $y_{R}$, the tidal deformability can be calculated as:[37,38]

$\displaystyle{\Lambda=}$ | $\displaystyle{{}\frac{16}{15}(1-2\mathcal{C})^{2}[2+2\mathcal{C}(y_{R}-1)-y_{R}]}$
---|---
| $\displaystyle{{}\\{2\mathcal{C}[6-3y_{R}+3\mathcal{C}(5y_{R}-8)]}$
|
$\displaystyle{{}+4\mathcal{C}^{3}[13-11y_{R}+\mathcal{C}(3y_{R}-2)+2\mathcal{C}^{2}(1+y_{R})]}$
|
$\displaystyle{{}+3(1-2\mathcal{C})^{2}[2-y_{R}+2\mathcal{C}(y_{R}-1)]\ln{(1-2\mathcal{C})}\\}^{-1}}$
(9)

TABLE I. Global calculations of binding energies of 603 even-even nuclei and
charge radii of 339 nuclei are compared with experimental data.[39,40]
Calculations are performed with the original SLy4, the refitted
SLy$4{{}^{\prime}}$, and three extended SLy4E parameter sets. The standard
root-mean-square (rms) deviations of binding energies $\sigma_{E}$ (in MeV)
and charge radii $\sigma_{r}$ (in fm) are listed.

Parameter SLy4 SLy$4{{}^{\prime}}$ global light heavy $\sigma_{E}(Z\leq 28)$
2.22 2.832 2.302 1.426 2.281 $\sigma_{E}(Z\leq 82)$ 3.18 2.582 2.002 2.381
2.373 $\sigma_{E}(all)$ 4.372 2.860 2.303 3.013 2.576 $\sigma_{r}(all)$ 0.028
0.0287 0.0244 0.0203 0.0244

_Results._ — We firstly performed global Skyrme-Hartree-Fock+BCS calculations
of 603 even-even nuclei. The obtained binding energies of 603 nuclei[39] and
charge radii of 339 nuclei[40] are compared with experimental data, shown in
Table I. For SLy4 force, the overall rms of deviations in binding energies is
about 4.38 MeV, which is very large. In particular, SLy4 seriously
underestimates the binding energies of heavy and superheavy nuclei
remarkably.[27] This undermine its descriptions of even larger systems such as
neutron stars. We also see that overall rms of the refitted SLy4′ is about
2.86 MeV, that is improved significantly. We also performed calculations with
extended SLy4 forces with an additional high-order density dependent term. The
SLy4Elight is excellent for light nuclei but its precision reduces towards
heavy nuclei. SLy4Eheavy is acceptable for both light and heavy nuclei.
SLy4Eglobal is further improved compared to SLy4Eheavy and can be used for
global calculations of finite nuclei with a rms of 2.57 MeV. The charge radii
are generally also improved by extended SLy4 forces. The SLy4Elight was
adopted for studies of collective excitations of weakly bound nuclei. [41] The
extended Skyrme forces based on SkM∗ and UNEDF0 are very successful for global
descriptions of finite nuclei and nuclear driplines,[42,43] but are not
successful for descriptions of neutron stars.

Fig. 1. (a) Pressure of pure neutron matter and (b) symmetry energy as a
function of nuclear density, corresponding to different Skyrme parameter sets,
see text for details.

Next we study the influences of high-order density dependent term on nuclear
equation of state. In Fig.1(a), the pressures of neutron matter by different
extended SLy4 forces are shown. For various modified SLy4 forces at the
saturation point ($\rho_{0}$=0.16 fm-3), the pressure ranges from 2.2 to 2.5
MeV. For densities lower than the saturation density, the pressures are almost
the same. In the high density region, we see SLy4′ gives the lowest pressure,
while SLy4 and SLy4Eglobal gives the highest pressure. The pressure is also
related to the sound of speed in dense nuclear matter.[23] We see the high-
order density dependent term can increase the pressure at high densities. This
is understandable that the repulsive high-order term increases the pressure at
high densities, which is particular important for descriptions of neutron
stars. The symmetry energies is related to EoS and has been extensively
studied by experiments of heavy ion collisions.[20] The symmetry energy at the
saturation density is around 32 MeV for SLy4 forces, while it has large
uncertainties at high densities. In our cases, it can be seen that SLy4′,
SLy4Elight, SLy4Eheavy give very soft symmetry energies at high densities. The
symmetry energies of SLy4Eglobal are slightly smaller than that of SLy4
results.

Fig. 2. Mass-radius relation obtained from EOS of different Skyrme parameter
sets. The dashed regions correspond to constrained mass and radius of the two
NSs in GW170817 at $90\%$ confidence level.[6]

We now solve TOV equation to study the properties of neutron stars with
various Skyrme forces. Note that we consider the neutron star matter is at the
$npe\mu$ equilibrium, which leads to a slightly soft EoS compared to the pure
neutron matter. Fig.2 shows the calculated mass-radius relation. The crucial
criteria for testing effective forces is to reproduce the observed maximum
mass of neutron stars, which is about 2 solar mass.[24] The soft EOS can not
support such a large mass. It can be seen that the maximum masses of all
forces are larger than 2 solar mass. The secret of series of successful SLy4
forces is that its $x_{2}$ parameter is $-1.0$.[23] SLy4′ results in a
slightly smaller maximum mass than that of SLy4. The $M-R$ relation from SLy4
and SLy4Eglobal are almost the same. We can see that SLy4′ results in the
smallest radii of 11.42 km for a NS with 1.4$M_{\bigodot}$. On the other hand,
SLy4 and SLy4Eglobal correspond to the largest radii of 11.67 and 11.70 km,
respectively. It can be understood that the radii of 1.4 $M_{\bigodot}$ and
the maximum NS mass are exactly correlated with the pressure at high
densities, rather than the symmetry energy. In contrast to precise
measurements of NS masses, there are large uncertainties in measurements of NS
radii. The Bayesian analysis of GW170817 binary NSs results in a radius of
11.9$\pm$1.4 km for two NSs.[6] We see that our results agree well with the
$M-R$ relation of GW170817 binary neutron stars. Another Bayesian analysis of
GW170817 provides a more stringent constraint of radii of
12.36${}^{+0.52}_{-0.38}$ km and 12.32${}^{+0.66}_{-0.43}$ km,
respectively.[7] The GW170817 measurement excluded very stiff EOS.

Fig. 3. Calculated tidal deformability $\Lambda$ as a function of NS’s mass,
with different Skyrme parameter sets. The inferred tidal deformability with
uncertainties from Bayesian analysis of GW180817 binary NSs are also shown.
The inferred tidal deformability of 1.4 $M_{\bigodot}$ is also shown.

Finally, the calculated tidal deformability $\Lambda$ with the five SLy4
parameter sets are shown in Fig. 3 as a function of NS’s mass. We see that
generally the tidal deformability is very large for NS with small masses. SLy4
and SLy4Eglobal result in the largest tidal deformability, which is related to
the largest $R_{1.4}$. On the other hand, SLy4′ results in smallest tidal
deformability, corresponding to the smallest $R_{1.4}$. The differences in
tidal deformability from various Skyrme forces are correlated with the $M-R$
relation. The Bayesian analysis of the GW170817 event suggests that at $90\%$
confidence level, the tidal deformabilities $\Lambda_{1,2}$ of the binary
stars are $\Lambda_{1}=255^{+416}_{-171}$ and $\Lambda_{2}=661^{+858}_{-375}$,
respectively.[7] The tidal deformability of 1.4 $M_{\bigodot}$ is estimated to
be $\Lambda_{1.4}=190^{+390}_{-120}$ at $90\%$ level.[6,7] We see that SLy4
and SLy4Eglobal results agree better with the observed tidal deformability of
GW170817.

_Conclusion._ — The observations of GW170817 binary neutron stars, including
the tidal deformability from gravitational waves for the first time, provided
a great opportunity for constraining the equation of state of dense nuclear
matter. It is known that the EoS of dense nuclear matter has large
uncertainties which is related to large uncertainties of nuclear interactions
at high densities.[10,44-46] There are also extensive studies to constrain EoS
or symmetry energies at high densities based on correlations with other
nuclear observables.[10,44-46] It is crucial to develop a unified and
consistent description of finite nuclei and neutron stars. We developed a
series of extended SLy4 forces with an additional high-order density dependent
term by fitting the properties of finite nuclei, as inspired by the Lee-Huang-
Yang formula. The original SLy4 force is widely used in studying neutron stars
with the maximum mass above 2.0 solar mass but is not precise for finite
nuclei. The refitted SLy4′ without high-order density dependence can improve
the descriptions of finite nuclei but its EOS becomes soft. The extended
SLy4Eglobal are satisfactory for finite nuclei and neutron stars
simultaneously. Our results agree well with the GW170817 observations. We
demonstrated that the essential role of the high-order density dependent term
in effective nuclear interactions at ultrahigh densities. Our unified
SLy4Eglobal descriptions provide a predictive model for studying neutron
stars, as well as new understandings of effective interactions. The realistic
studies of neutron stars including superfluidity and non-uniform crusts are of
interests in the future. We expect that the more precise observations of
gravitational waves in the future can provide more stringent constraints on
nuclear interactions in ultradense nuclear matter.

## References

[1] Baade W and Zwicky F 1934 National Academy of Sciences 20 259

[2] James M L and Madappa P 2016 Physics Reports 621 127

[3] Oertel M, Hempel M, Klähn T and Typel S 2017 Rev. Mod. Phys. 89 015007

[4] Rikovska S J, Miller J C, Koncewicz R, Stevenson P D and Strayer M R 2003
Phys. Rev. C. 68 034324

[5] Abbott B P $et\ al.$ 2017 Phys. Rev. Lett. 119 161101

[6] Abbott B P $et\ al.$ 2018 Phys. Rev. Lett. 121 161101

[7] Fasano M, Abdelsalhin T, Maselli A and Ferrari V 2019 Phys. Rev. Lett. 123
141101

[8] Abbott B P $et\ al.$ 2019 Phys. Rev. X 9 011001

[9] Hinderer T 2008 Astrophys. J. 677 1216

[10] Burgio G F and Vidaña I 2018 Universe 6 119

[11] Douchin F and Haensel P 2001 Astro. Astrophys. 380 151

[12] James M L and Madappa P 2001 Astrophys. J. 550 426

[13] Stone J R, Stevenson P D, Miller J C and Strayer M R 2002 Phys. Rev. C.
65 064312

[14] Ishii N, Aoki S and Hatsuda T 2007 Phys. Rev. Lett. 99 022001

[15] Tews I, Davoudi Z, Ekstrom A, Holt J D and Lynn J E 2020 J. Phys. G. 47
103001

[16] Dutra M, Lourenço O, SáMartins J S, Delfino, A, Stone J and Stevenson P D
2012 Phys. Rev. C. 85 035201

[17] Li A, Zhu Z Y, Zhou E P, Dong J M, Hu J N and Xia C J 2020 JHEAp 28 19

[18] Tamii, A. $et\ al.$ 2011 Phys. Rev. Lett. 107 062502

[19] Abrahamyan $et\ al.$ 2012 Phys. Rev. Lett. 108 112502

[20] Tsang M B, Zhang Y X, Danielewicz P, Famiano M, Li Z X, Lynch W G and
Steiner A W 2009 Phys. Rev. Lett. 102 122701; Danielewicz P, Lacey R and Lynch
W G 2002 Science 298 1592; Li B A, Chen L W and Che M K 2008 Phys. Rep. 464
113; Xiao Z G, Li B A, Chen L W, Yong G C, and Zhang M 2009 Phys. Rev. Lett.
102 062502

[21] Skyrme T H R 1956 Philos. Mag 1 1043

[22] Vautherin D and Brink D M 1972 Phys. Rev. C. 5 626

[23] Chabanat E, Bonche P, Haensel P, Meyer J and Schaeffer R 1997 Nuclear
Physics A 627 710

[24] Chabanat E, Bonche P, Haensel P, Meyer J and Schaeffer R 1998 Nuclear
Physics A 635 231

[25] Antoniadis John $et\ al.$ 2013 Science 340 1233232

[26] Erler J, Horowitz C J, Nazarewicz W, Rafalski M and Reinhard P G 2013
Phys. Rev. C. 87 044320

[27] Xiong X Y, Pei J C and Chen W J 2016 Phys. Rev. C. 93 024311

[28] Negele J W 1982 Rev. Mod. Phys. 54 913

[29] Grange P, Lejeune A, Martzolff M and Mathiot J F 1989 Phys. Rev. C. 40
1040

[30] Huang K and Yang C N 1957 Phys. Rev. 105 767

[31] Lee T D and Yang C N 1957 Phys. Rev. 105 1119

[32] DeDominicis C and Martin P C 1957 Phys. Rev. 105 1419

[33] Altmeyer A, Riedl S, Kohstall C, Wright M J, Geursen R, Bartenstein M,
Chin C, Denschlag J H and Grimm R 2017 Phys. Rev. Lett. 98 040401

[34] Agrawal B K, Dhiman S K and Kumar R 2006 Phys. Rev. C. 73 034319

[35] Tolman R C 1939 Phys. Rev. 55 364

[36] Oppenheimer J R and Volkoff G M 1939 Phys. Rev. 55 374

[37] Hu J N, Bao S S, Zhang Y, Nakazato K I, Sumiyoshi K and Shen H 2020 Prog.
Theo. Expt. Phys. 2020 10.1093/ptep/ptaa016

[38] Hinderer T, Lackey B D, Lang R N and Read J S 2010 Phys. Rev. D. 81
123016

[39] Wang M, Audi G, Kondev G G, Huang W J, Naimi S and Xu X 2017 Chin. Phys.
C. 41 030003

[40] Angeli I and Marinova K P 2013 Atomic Data and Nuclear Data Tables 99 69

[41] Wang K, Kortelainen M and Pei J C 2017 Phys. Rev. C. 96 031301

[42] Zuo Z W, Pei J C, Xiong X Y and Zhu Y 2018 Chin. Phys. C. 42 064106

[43] Chai Q Z, Pei J C, Fei N and Guan D W 2020 Phys. Rev. C. 102 014312

[44] Malik T, Agrawal B K, De J N, Samaddar S K, Providência C, Mondal C and
Jha T K 2019 Phys. Rev. C. 99 052801

[45] Tsang C Y, Tsang M B, Danielewicz P, Fattoyev F J, Lynch W G 2019 Phys.
Lett. B. 796 1

[46] Zhang Y X, Liu M, Xia C J, Li Z X and Biswal S K 2020 Phys. Rev. C. 101
034303