---
author:
- Dominik Gerstung
- Norbert Kaiser
- Wolfram Weise
title: 'Hyperon-nucleon three-body forces and strangeness in neutron stars'
---

Introduction {#sec:intro}
============

The existence of heavy neutron stars with masses around 2$M_\odot$ [@Demorest2010; @Antoniadis2013; @Fonseca2016; @Cromartie2019] sets strong constraints on the equation of state (EoS) of dense baryonic matter. This EoS must be sufficiently stiff, i.e. the pressure $P(\cal{E})$ at energy densities $\cal{E}\sim$ 1 GeV/fm$^3$ must be large enough to support such massive compact objects against gravitational collapse. The detection of gravitational wave signals from two merging neutron stars [@Abbott2017] adds further important information on the EoS, by providing limits for the tidal deformability and for neutron star radii [@Most2018; @De2018].

The composition and properties of strongly interacting matter at high baryon densities is a topic of continuing interest. Various options are under discussion. A time-honored description of matter in the core of neutron stars uses hadronic degrees of freedom (baryons and mesons) with strong many-body correlations [@APR1998]. A modern version that arrives at a similar EoS, consistent with observations, is based on a chiral nucleon-meson field theory combined with functional renormalization group methods [@DW2015; @DW2017]. In this latter approach neutron star matter can be represented as a relativistic Fermi liquid in the sense of Landau theory [@FW2019]. The density-dependent leading Landau parameters show the characteristic behaviour of a strongly correlated fermionic many-body system, albeit less extreme in comparison with another well-known Fermi system at low temperature, namely liquid $^3$He. Alternative descriptions of neutron star matter involve a (possibly smooth) transition from hadronic matter to some form of quark matter [@Baym2018; @Baym2019; @McLerran2019; @FFM2019].

The required stiffness of the EoS implies strong restrictions on the appearance of hyperons in neutron star matter and for the underlying hyperon-nuclear interactions. A naive introduction of $\Lambda$ hyperons, with only two-body $\Lambda N$ interactions, would appear to be energetically favorable by replacing neutrons at baryon densities around 2-3 $\rho_0$ (in terms of the equilibrium density of normal nuclear matter, $\rho_0 = 0.16$ fm$^{-3}$). However, then the EoS of neutron star matter would become far too soft to satisfy the two-solar-mass constraint [@Djapo2010; @Hell2014; @Lonardoni2015]. A phenomenological analysis [@BS2015], taking into account constraints from both hypernuclear physics and neutron star observations, points out indeed that the interaction between $\Lambda$ hyperons and dense matter has to become repulsive at densities below three times $\rho_0$.

The aim of the present work is to investigate whether and to what extent a microscopic description of hyperon-nucleon two- and three-body forces is capable of providing the necessary repulsion in dense baryonic matter. The key quantities to be calculated are the hyperon ($\Lambda$ and $\Sigma$) single-particle potentials $U_{\Lambda,\Sigma}(\rho)$ in nuclear and neutron matter. The density dependence of the resulting chemical potential $\mu_\Lambda(\rho)= M_\Lambda +U_{\Lambda}(\rho)$ of a $\Lambda$ in neutron matter indicates whether it is favourable to replace neutrons by $\Lambda$ hyperons in the core of neutron stars.

Our starting point is the hyperon-nucleon interaction derived from SU(3) chiral effective field theory (ChEFT) at next-to-leading order (NLO) [@Haidenbauer2013]. This interaction includes exchanges of one and two pseudoscalar octet mesons and four-baryon contact terms with SU(3)-symmetric low-energy constants fitted to the (admittedly scarce) hyperon-nucleon scattering data. Three-body forces enter first at next-to-next-to-leading order (NNLO) in this scheme [@Petschauer2016]. Estimating the strength of contact terms in these forces through the contributions of explicit baryon decuplet resonances in intermediate states, these three-body interactions are promoted from NNLO to NLO [@Petschauer2017] due to the small decuplet-octet baryon mass-splitting. We shall adopt here the same strategy. Previous calculations of hyperon-nuclear potentials using Brueckner theory with these interactions as input have been reported in refs.[@Kohno2018; @Petschauer2016b; @Haidenbauer2017]. The present study builds on these results, with extensions in several directions. Special attention will be paid to the importance of the $\Lambda N\leftrightarrow\Sigma N$ transition potential, a necessary ingredient of any such calculation. In particular, the $\Lambda NN\leftrightarrow\Sigma NN$ coupled channels in the three-body sector, translated into density-dependent effective two-body potentials as in ref.[@Petschauer2017], will be treated explicitly when solving coupled-channel Bethe-Goldstone equations.

Basic baryonic interactions {#sec:int}
===========================

Hyperon-nucleon two-body interactions
-------------------------------------

For the description of the hyperon-nucleon two-body interaction, SU(3) chiral effective field theory is used up to NLO with the Weinberg power counting applied to the potential, as reported in detail in ref.[@Haidenbauer2013]. The leading order potential involves the exchange of a single pseudoscalar octet meson $(\pi, K, \eta)$ and non-derivative four-baryon contact terms. At NLO two-meson exchange diagrams at one-loop level arise together with additional contact terms with explicit momentum dependence of order $p^2$ (see fig.\[fig:pwrcounting\]). The set of SU(3)-symmetric contact terms represents unresolved short-distance dynamics. The corresponding low-energy constants (LECs) are fitted to low-energy hyperon-nucleon scattering data and the hypertriton ${}^3_\Lambda H$ binding energy (and partly to elastic $NN$-scattering phase shifts). SU(3) symmetry breaking is none- theless incorporated through the physical masses of the baryons and of the exchanged pseudoscalar mesons. Further details on input parameters can be found in ref.[@Haidenbauer2013].

![Leading and next-to-leading order diagrams representing the baryon-baryon interaction potential. Solid and dashed lines denote octet baryons ($N,\Lambda,\Sigma$) and mesons ($\pi,K,\eta$), respectively.[]{data-label="fig:pwrcounting"}](FBBcont "fig:"){width="\feynwidth"} ![Leading and next-to-leading order diagrams representing the baryon-baryon interaction potential. Solid and dashed lines denote octet baryons ($N,\Lambda,\Sigma$) and mesons ($\pi,K,\eta$), respectively.[]{data-label="fig:pwrcounting"}](FBBtree "fig:"){width="\feynwidth"}\
![Leading and next-to-leading order diagrams representing the baryon-baryon interaction potential. Solid and dashed lines denote octet baryons ($N,\Lambda,\Sigma$) and mesons ($\pi,K,\eta$), respectively.[]{data-label="fig:pwrcounting"}](FBBcontNLO "fig:"){width="\feynwidth"} ![Leading and next-to-leading order diagrams representing the baryon-baryon interaction potential. Solid and dashed lines denote octet baryons ($N,\Lambda,\Sigma$) and mesons ($\pi,K,\eta$), respectively.[]{data-label="fig:pwrcounting"}](FBBpb "fig:"){width="\feynwidth"}![Leading and next-to-leading order diagrams representing the baryon-baryon interaction potential. Solid and dashed lines denote octet baryons ($N,\Lambda,\Sigma$) and mesons ($\pi,K,\eta$), respectively.[]{data-label="fig:pwrcounting"}](FBBcb "fig:"){width="\feynwidth"}![Leading and next-to-leading order diagrams representing the baryon-baryon interaction potential. Solid and dashed lines denote octet baryons ($N,\Lambda,\Sigma$) and mesons ($\pi,K,\eta$), respectively.[]{data-label="fig:pwrcounting"}](FBBtriR "fig:"){width="\feynwidth"}![Leading and next-to-leading order diagrams representing the baryon-baryon interaction potential. Solid and dashed lines denote octet baryons ($N,\Lambda,\Sigma$) and mesons ($\pi,K,\eta$), respectively.[]{data-label="fig:pwrcounting"}](FBBtriL "fig:"){width="\feynwidth"}![Leading and next-to-leading order diagrams representing the baryon-baryon interaction potential. Solid and dashed lines denote octet baryons ($N,\Lambda,\Sigma$) and mesons ($\pi,K,\eta$), respectively.[]{data-label="fig:pwrcounting"}](FBBfoot "fig:"){width="\feynwidth"}

Solving the coupled-channel Lippmann-Schwinger equations using these chiral potentials involves a regulator function of the form $\exp[-(p'^4+p^4)/\lambda^4]$ depending on the momenta, $p$ and $p'$, of the in- and outgoing baryons. The cutoff $\lambda$ is chosen in a reasonable range 500 ... 650 MeV.

A characteristic feature of the resulting $\Lambda N$ S-wave phase shifts that distingushes the chiral EFT-based potential from all earlier phenomenological interactions, such as Jülich’04 [@HM2005] and Nijmegen NSC97f [@Rijken1999], is the following. As there is no (leading-order) one-pion exchange $\Lambda N$ interaction term, the driving intermediate range attraction in this channel comes from the iterative process $\Lambda N \rightarrow \Sigma N \rightarrow \Lambda N$, involving the exchange of two pions. The $\Lambda\Sigma\pi$ coupling as given by SU(3) symmetry implies strong attraction through this second-order process that is partly balanced by the strong short-distance repulsion. The net attraction is seen in the low-energy $^3S_1$ $\Lambda N$ phase shift. The $^1S_0$ phase shift also shows this attractive behavior at low energies but changes sign around a lab momentum of 600 MeV and becomes repulsive. This momentum-dependent repulsion is generated by the short-distance terms at NLO, as can be clearly seen by selectively turning off the $\Lambda\Sigma\pi$ coupling [@Haidenbauer2017]. The repulsive effects in the $\Lambda N$ interaction at high momenta have their direct impact on the $\Lambda$ single-particle potentials in nuclear and neutron matter which change from attraction to repulsion at densities around or below $2\,\rho_0$, a feature not matched by any of the phenomenological potentials. However, as we shall see, such two-body repulsion is still not strong enough to prevent $\Lambda$ hyperons from occuring at baryon densities encountered in neutron stars and thereby softening the EoS to an unwanted degree.

Three-body interactions of hyperons with nucleons
-------------------------------------------------

Three-nucleon forces are an important element in the quantitative understanding of nuclear few- and many-body systems. Likewise, the $\Lambda NN$ three-body interactions is expected to play a significant role in hypernuclear systems. Our interest in the present work is directed, in particular, to the behaviour of $\Lambda NN$ and $\Sigma NN$ three-body forces (3BF) in dense baryonic matter, where their influence is expected to grow continuously with increasing density.

In the systematic expansion of SU(3) chiral EFT, baryonic three-body forces appear at NNLO. Their detailed derivation is given in ref.[@Petschauer2016] on which we shall build. The leading 3BF diagrams are shown in fig.\[fig:3BF\]. They fall into three classes: six-baryon contact terms, one-meson exchange pieces, and two-meson exchange pieces. The corresponding potentials are denoted by $V^{(0)}, V^{(1)}$ and $V^{(2)}$, respectively. The contact potential $V^{(0)}$ involves all possible combinations of baryon spin operators $\vec \sigma_1,  \vec \sigma_2, \vec \sigma_3$. The one-meson exchange 3BF potential has the following generic form: $$V^{(1)}_i = \frac{1}{2f^2} \frac{\vec\sigma_1\cdot\vec q}{\vec q^{\,2}+m_i^2} (
A \vec\sigma_3+ B \, \mathrm i \vec\sigma_2\times\vec\sigma_3)\cdot\vec q
\,,$$ where the inverse squared pseudoscalar meson decay constant $f^{-2} 4.8\,$fm$^2$ determines the interaction strength. The index $i\in \{\pi^0,\pi^+,\pi^-,K^+,K^-,K^0,\bar K^0,\eta\} $ refers to the exchanged meson with mass $m_i$, carrying the momentum transfer $\vec q$. The two parameters $A$ and $B$ are combinations of low-energy constants.

![ Leading three-baryon interactions: contact term, one-meson exchange and two-meson exchange. \[fig:3BF\] ](FBBBcont "fig:")![ Leading three-baryon interactions: contact term, one-meson exchange and two-meson exchange. \[fig:3BF\] ](FBBB1ME "fig:")![ Leading three-baryon interactions: contact term, one-meson exchange and two-meson exchange. \[fig:3BF\] ](FBBB2ME2 "fig:")

On the other hand, the two-meson exchange potential (depicted in fig.\[fig:3BF\]) has following generic form: $$\begin{aligned}
V^{(2)}_{ij} ={} & -\frac{1}{4f^4} \frac{\vec\sigma_1\cdot\vec q\ ~~\vec\sigma_3\cdot\vec q'}{(\vec q^{\,2}+m_i^2)(\vec q'^{\,2}+m_j^2)} \\
&\qquad\times\Big(A' + B'\vec q\cdot\vec q' +D'\, \mathrm i\vec\sigma_2\cdot(\vec q\times\vec q')\Big) \,. {\addtocounter{equation}{1}\tag{\theequation}}\end{aligned}$$ The parameters $A',B'$ and $D'$ are combinations of low-energy constants determined in the meson-baryon subsector. The transferred momenta $\vec q$ and $\vec q'$ are carried by the mesons with masses $m_i$ and $m_j$, respectively. The complete three-body potential is constructed by summing the contributions from all distinguishable diagrams with all possible exchanged pseudoscalar mesons, including baryon exchange operations where needed.

As they stand, the 3BF terms at NNLO introduce a prohibitively large number of low-energy constants in the contact potentials $V^{(0)}$. The existing data-base is too limited and does not permit a meaningfully constrained determination of these parameters. However, their number can be reduced substantially by employing an approximate scheme referred to as decuplet saturation: contact vertices are resolved by propagating explicit decuplet baryons in intermediate states. Furthermore, it has the welcome feature of promoting three-baryon forces from NNLO to NLO within the chiral hierarchy. In essence, the NNLO diagrams in fig.\[fig:3BF\] are replaced by the NLO diagrams of fig.\[fig:3BFdec\]. In our subsequent application, two of the in- or outgoing baryons are nucleons (proton, neutron), while the third one is a strangeness $S=-1$ hyperon ($\Lambda$ or $\Sigma$) and only the (long-range) pion-exchange is considered. The decuplet intermediate states are then either $\Sigma^*(1385)$ or $\Delta(1232)$.

![Three-baryon forces with explicit decuplet baryons in the intermediate states (represented by double lines) [@Petschauer2017].[]{data-label="fig:3BFdec"}](FBBBdeccont "fig:")![Three-baryon forces with explicit decuplet baryons in the intermediate states (represented by double lines) [@Petschauer2017].[]{data-label="fig:3BFdec"}](FBBBdec1MEu "fig:")![Three-baryon forces with explicit decuplet baryons in the intermediate states (represented by double lines) [@Petschauer2017].[]{data-label="fig:3BFdec"}](FBBBdec1MEd "fig:")![Three-baryon forces with explicit decuplet baryons in the intermediate states (represented by double lines) [@Petschauer2017].[]{data-label="fig:3BFdec"}](FBBBdec2ME "fig:")

In what follows we use the detailed formalism developed in refs.[@Petschauer2016; @Petschauer2017]. Within the decuplet dominance approximation there are only three remaining constants to be determined: the coupling strength $C$ of the transition vertex between an octet baryon and a decuplet baryon with a pseudoscalar meson absorbed or emitted; and two coupling constants, $H_1$ and $H_2$, for the four-point vertices connecting three octet baryon lines and one decuplet baryon [@Petschauer2016; @Petschauer2017]. The coupling $C$ is uniquely determined by the $\pi N\rightarrow\Delta$ transition vertex which is in turn constrained by the decay width $\Gamma(\Delta\rightarrow \pi N)\simeq 115$ MeV. The (large-$N_c$) value $C = 3g_A/4 \simeq 0.95$ with the nucleon axial vector coupling constant, $g_A = 1.26$, is well compatible with this constraint. The remaining two constants, $H_{1,2}$ of dimension $(length)^2$, are still free to choose. We use hypernuclear phenomenology in order to restrict their possible values.

Next, the three-body interactions in the $\Lambda NN\!\rightarrow\!\Lambda NN$, $\Lambda NN\leftrightarrow\Sigma NN$, and $\Sigma NN \rightarrow\Sigma NN$ coupled channels are translated into effective density-dependent hyperon-nucleon potentials, which additionally enter the coupled Bethe-Goldstone equations. This is done as in refs.[@Petschauer2016; @Petschauer2017; @Haidenbauer2017; @Holt2010] by integrating one of the two nucleons over the filled Fermi sea, as illustrated schematically in fig.\[fig:3BFeff\]:

![Density-dependent in-medium potentials generated from 3-body forces. Upper left: contact term; upper right: one-meson exchange term; lower part: two-meson exchange term. Loops with arrows indicate integration over the filled nucleon Fermi sea. []{data-label="fig:3BFeff"}](FBBeff)

$$\label{eq:red}
V_{12}^{\text{eff}}(\rho_p,\rho_n) = \sum_{N=p,n}\int\limits_{|\vec k|\leq k_F^{(N)}}\frac{\mathrm d^3k}{(2\pi)^3}\,{\mathrm tr}_3V_{123} \,,$$

where ${\mathrm tr}_3$ denotes the spin-trace over the third particle and the summation over nucleons (protons and neutrons) in the Fermi sea is performed. The resulting effective two-body potentials depend explicitly on the proton and neutron densities, $$\rho_p = {(k_F^{(p)})^3\over 3\pi^2}~,~~~\rho_n = {(k_F^{(n)})^3\over 3\pi^2}~.$$ As a first example, we present expressions for the $\Lambda N\to \Lambda N$ effective two-body potential in (asymmetric) nuclear matter with density $\rho = \rho_p + \rho_n$. Only the $\Lambda n $ effective potential in the nuclear medium needs to be given since the $\Lambda p$ potential follows by simply interchanging $\rho_n\leftrightarrow\rho_p$. The three parts have the following structure [@Petschauer2017]:\
Contact term: $$V^{(0)}_{\text{eff}}(\Lambda n) = {(H_1+3H_2)^2\over 18 \Delta}\,(\rho_n+2\rho_p)~,$$ with the decuplet-octet baryon mass difference $\Delta$, for which we take an average value of $\Delta = 270$ MeV.\
One-pion exchange part: $$V^{(1)}_{\text{eff}}(\Lambda n) =
{g_A C\over 9f^2 \Delta}(H_1+3H_2)\left(\rho_n+2\rho_p - {m_\pi^2\over 2\pi^2}\Gamma_0\right),$$ where $\Gamma_0(p;k_F^{(N)})$ is a function of the $\Lambda N$ center-of-mass momentum $p$, and the proton or neutron Fermi momentum $ k_F^{(N)}$, given explicitly in ref.[@Petschauer2017].\
Two-pion exchange part: $$V^{(2)}_{\text{eff}}(\Lambda n) = {g_A^2 C^2\over 6f^4 \Delta}\left[\rho_n+2\rho_p+ {\cal F}(p,q;k_F^{(n)},k_F^{(p)})\right]\,,$$ where the function ${\cal F}$ depends additionally on the momentum transfer $q$ in the reduced two-body system. It is a lengthy expression involving the spin-orbit operator $\mathrm i \vec \sigma_2\cdot(\vec q\times \vec p\,)$, that is explicitly given in eq.(46) of ref.[@Petschauer2017].

Next, we present some selected expressions for the effective $\Lambda N \rightarrow \Sigma N$ interaction. As representative examples, consider the direct terms, shown in fig.\[fig:3BFeff\], of the $\Lambda n\rightarrow \Sigma^- p$ effective transition potential. Exchange terms are also included in the actual calculations but are not displayed here for simplicity. In the case of symmetric nuclear matter, the corresponding lengthy formulas can be found in subsection III.D of ref.[@Petschauer2017].\
Contact term: $$\begin{aligned}
 && V^{(0)}_\text{{eff}}(\Lambda n\!\rightarrow\!\Sigma^- p) = {1\over 6\sqrt{6}\Delta}\bigg\{(\rho_p\!+\!\rho_n) \big[4H_2^2\!-\!(H_1\!+\!H_2)^2\big]\nonumber\\ &&+ \vec\sigma_1\!\cdot\!\vec\sigma_2\bigg[\bigg(H_1^2\!+\!H_2^2\!+\!{10\over 3}H_1 H_2\bigg)\rho_n+ (H_1\!+\!3H_2)^2{ \rho_p\over 3}\bigg]\bigg\}\,. \nonumber\\ &&\end{aligned}$$ One-pion exchange part: $$\begin{aligned}
&&V^{(1)}_\text{{eff}}(\Lambda n\!\rightarrow\!\Sigma^- p)={\sqrt{6}g_A C \over 27f^2\Delta}\, {\vec\sigma_1\cdot\vec q\,\vec\sigma_2\cdot\vec q\over m_\pi^2+\vec q^2}\nonumber\\ && \qquad\qquad \times \Big[(2H_1-3H_2)\rho_n-(2H_1+9H_2)\rho_p\Big]\,.\end{aligned}$$ Two-pion exchange part: $$\begin{aligned}
V^{(2)}_\text{{eff}}(\Lambda n\!\rightarrow\!\Sigma^- p)=-{16g_A D C^2\over 9\sqrt{6}f^4\Delta}{\vec\sigma_1\!\cdot\!\vec q\,\vec\sigma_2\!\cdot\!\vec q\over( m_\pi^2+\vec q^2)^2}(\rho_p+\rho_n) \vec q^2 , \nonumber\\ &&\end{aligned}$$ where the SU(3) axial vector coupling constant $D$ stems from the $\Lambda\Sigma\pi$ vertex. For the $\Lambda p\rightarrow \Sigma^+ n$ channel, corresponding expressions hold with $\rho_n\leftrightarrow\rho_p$ interchanged.

Finally, the ladder summations in the coupled hyperon-nucleon channels require also as input the $\Sigma N\leftrightarrow\Sigma N$ effective two-body potentials in different charge combinations. As representative examples we list here expressions for the in-medium $\Sigma^-n\rightarrow\Sigma^- n$ potential:\
Contact term: $$\begin{aligned}
&&V^{(0)}_\text{{eff}}(\Sigma^- n) = 
 {1\over 6\Delta}\,\Big\{\big[(H_1+H_2)^2\rho_p + 4H_1^2\rho_n\big]
\nonumber\\ && \qquad \qquad\qquad\qquad -{1\over 3}\,\vec\sigma_1\cdot\vec\sigma_2\,(H_1+H_2)^2\,\rho_p\Big\}\, .\end{aligned}$$ One-pion exchange part: $$\begin{aligned}
&&V^{(1)}_\text{{eff}}(\Sigma^- n)=
-{4g_A C\over 9f^2\Delta}\,{\vec\sigma_1\!\cdot\!\vec q\,\vec\sigma_2\!\cdot\!\vec q\over m_\pi^2+\vec q^2}\Big[H_1\rho_n \!+\!(H_1\!+\!H_2)\rho_p\Big]\,,\nonumber\\ &&\end{aligned}$$ Two-pion exchange part: $$V^{(2)}_\text{{eff}}(\Sigma^- n)=-{8g_A F C^2\over 9f^4\Delta}\,{\vec\sigma_1\cdot\vec q\,\vec\sigma_2\cdot\vec q\over( m_\pi^2+\vec q^2)^2}\,(\rho_n+\rho_p)\,\vec q^2\,,$$ where now the SU(3) axial vector coupling constant $F$ enters. The $\Sigma^+ p\rightarrow\Sigma^+ p$ effective potential is obtained by interchanging $\rho_n\leftrightarrow\rho_p$. The $\Sigma^0 n\rightarrow\Sigma^0 n$ and $\Sigma^0 p\rightarrow\Sigma^0 p$ effective potentials have similar expressions, except that for their two-pion exchange parts the direct term as shown in fig.\[fig:3BFeff\] vanishes and only an exchange diagram contributes in this sector.

The examples just discussed demonstrate the following. In the diagonal $\Lambda N\rightarrow\Lambda N$ potential the two low-energy constants, $H_1$ and $H_2$, those associated with the three-body contact term using decuplet dominance, appear in the combination $H_1 + 3H_2$. This leaves only one parameter in this channel, as it has been the case in previous calculations of $\Lambda$ single-particle potentials in nuclear and neutron matter performed in refs.[@Kohno2018; @Haidenbauer2017]. In contrast, the explicit inclusion of $\Lambda NN\leftrightarrow\Sigma NN$ coupled channels introduces $H_1$ and $H_2$ in various different combinations so that these parameters have to be constrained independently.

Brueckner-Hartree-Fock approach to hyperon single-particle potentials
=====================================================================

The two- and three-body potentials described in the previous section are used[^1] as input for computing a hyperon-nucleon $G$-matrix within Brueckner theory at first order in the hole-line expansion (the Brueckner-Hartree-Fock approximation [@Day1967]). We focus on the single-particle potentials, i.e. the self-energies of hyperons in nuclear matter. The relevant formalism is briefly summarized below. For more details we refer to refs.[@Schulze1998; @Rijken1999; @Vidana2000; @Kohno2000].

The Brueckner reaction matrix or G-matrix is determined by solving the coupled-channel Bethe-Goldstone equation (in symbolic form) $$G(\omega) = V + V \frac Q{e(\omega)+\mathrm i\epsilon} G(\omega)\,.$$ The energy denominator $e(\omega)$ depends on the starting energy $\omega$. The Pauli blocking operator $Q$ excludes particles in intermediate states from scattering into the filled Fermi sea. The potential $V$ is a matrix (labeled by the outgoing and ingoing two-baryon channels) including two-body and three-body contibutions, with the latter given in their density-dependent effective two-body form.

After angle-averaging (see appendix A for details) the Bethe-Goldstone equation decomposes into partial waves with total angular momentum $J$: $$\begin{aligned}
\label{eq:BGE}
&&G^{J}_{\alpha\beta}(p',p; P,\omega) = V^{J}_{\alpha\beta}(p',p)  +\\
&&\sum_{\nu}\int_0^\infty\frac{\mathrm d k\, k^2}{(2\pi)^3} \, V^{J}_{\alpha\nu}(p',k) 
\frac{\bar Q_{\nu}(P,k)}{\bar e_{\nu}(P,k;\omega)+\mathrm i \epsilon}G^{J}_{\nu\beta}(k,p;P,\omega) \,.\nonumber \end{aligned}$$ Here, $P$ is the total momentum of the two baryons. The indices $\alpha, \beta, \nu$ represent complete sets of channels, including partial wave quantum numbers and pairs of interacting baryons, $(B_1,B_2)$. In our case of interest, $B_1$ in the initial state is a $\Lambda$ hyperon, while $B_2$ is a nucleon within its Fermi sea. The channel coupling turns $B_1$ into a $\Sigma$ which subsequently interacts with the nuclear medium and turns back into a $\Lambda$ in the final state.

The $G$-matrix elements of eq.(\[eq:BGE\]) are calculated at the on-shell starting energy $$\begin{aligned}
\omega &=& E_{B_1}(p_1) +  E_{B_2}(p_2)~,\nonumber\\
\mathrm{with}~~~~~~E_{B_i}(p) &=& M_i + {p^2\over 2M_i} + {\operatorname{Re}}U_{B_i}(p)~,\end{aligned}$$ where $M_i$ is the mass of baryon $B_i$. Hence, the energy denominator $e(\omega)$ (see appendix A) requires the single-particle potentials for both hyperons and nucleons as an input.

In Brueckner-Hartree-Fock approximation the single-particle potential for a baryon $B_1\in\{\Lambda,\Sigma^+,\Sigma^0,\Sigma^-,p,n\}$ interacting with the Fermi sea of nucleons $B_2\in\{p,n\}$ is then given by: $$\begin{aligned}
&U_{B_1}(p_1) = 
\sum_{\alpha}\left[1+\delta_{B_1B_2}(-1)^{L+S}\right]\frac{(1+M_2/M_1)^3}{2}~\times\nonumber \\
& \sum_{J}(2J+1)
\int_{k_\mathrm{min}}^{k_\mathrm{max}}\! \frac{\mathrm dk\, k^2}{(2\pi)^3}\,W(p_1,k) \, G_{\alpha\alpha}^{J}(k,k;\bar P,\omega) \,.
{\addtocounter{equation}{1}\tag{\theequation}}\label{eq:U}\end{aligned}$$ The weight function $W(p_1,k)$ resulting from the angular averaging procedure is specified in appendix A, together with the integration boundaries $k_\mathrm{min}$ and $k_\mathrm{max}$. Evidently, the calculation of the single-particle potentials depends on the single-particle potentials themselves as they appear in the energy denominator. Therefore eqs.  and must be solved self-consistently by iteration until convergence is reached. For the nucleon sector itself the $NN$ input potential is taken from ref.[@Entem2003], employing SU(2) chiral perturbation theory at fourth order (N3LO). Chiral three-nucleon forces at order N2LO are also included.

The computations are carried out using the so-called continuous choice, where single-particle potentials enter the energy denominator of intermediate states for all momenta, below and above the Fermi sea). From previous work [@Schulze1998] it is known that the continuous choice for intermediate states is preferable over the simpler gap choice (with potentials set to zero for momenta above the Fermi sea), because it allows for a reliable determination of the single-particle potentials including their imaginary parts. Some further technical details are given in Appendix B.

Results {#sec:res}
=======

Single-particle potentials of a $\Lambda$ hyperon in nuclear and neutron matter {#subsec:spp}
-------------------------------------------------------------------------------

The next step is now to calculate the single-particle potential, $U_\Lambda(p_1=0,\rho)$, for a $\Lambda$ at rest in symmetric nuclear matter as a function of baryon density $\rho=\rho_p+\rho_n$. Using as input the NLO hyperon-nucleon potentials together with the three-body interactions as specified in Section \[sec:int\], it is still necessary to constrain the two free parameters of the three-body contact terms. We recall that the calculation explicitly includes the $\Lambda NN \leftrightarrow\Sigma NN$ coupled channels, so the respective constants $H_1$ and $H_2$ enter independently.

![image](USNM)

![image](3BFH1H2)

A key constraint comes from hypernuclear phenomenology. Bulk properties and the shell structure of $\Lambda$ hypernuclei are well described by a $\Lambda$-nuclear single-particle potential parametrized in Wood-Saxon form, with an accurately determined potential depth at nuclear central densities [@Gal2016]: $$U_\Lambda(\rho \simeq \rho_0) = - 30~ \mathrm{MeV}~,
\label{eq:USNM}$$ corresponding to about half of the attractive strength of the nucleon single-particle potential at nuclear matter saturation density, $\rho_0 = 0.16$ fm$^{-3}$. This constraint will be satisfied by any of the calculations to be described below.

Fig.\[fig:ULambda1\] shows instructive examples of the calculated\
$U_\Lambda(p_1=0;\rho)$ in nuclear matter at low densities. Notably, the result using only the two-body chiral NLO potential would actually produce too much attraction in the $\Lambda$ potential, causing an overbinding of hypernuclei. This calls for repulsive effects beyond two-body interactions, with three-body forces being the natural extension. The three-body effects are moderate (a few MeV) at densities $\rho\simeq\rho_0$, but their relative importance grows continuously as the density increases. At $\rho \simeq 2\,\rho_0$ the repulsive $\Lambda NN$ 3BF is essential in order to turn the sign of $U_\Lambda$ from attractive to repulsive. The result in fig.\[fig:ULambda1\] is produced using a combination of three-body constants, $(H_1,H_2)$, that evidently satifies the constraint eq.(\[eq:USNM\]). However, this choice of $H_1$ and $H_2$ is not uniquely determined. It turns out that any pair $(H_1,H_2)$ that lies on either of the two solid curves in fig.\[fig:H1H2\] is compatible with the condition, $U_\Lambda(\rho=\rho_0) = - 30$MeV. On these two constrained curves one can identify points $(H_1,H_2)$ for which the $\Lambda$ single-particle potential at a higher density, $U_\Lambda(\rho=3\rho_0)$, becomes maximally repulsive (while avoiding instabilities of the numerical procedure): $(H_1,H_2) = (-2.2,0)$ and $(-2.47,1.2)$ in units of the squared inverse pion decay constant, $f^{-2} =4.8\,$fm$^2$. These are the parameters of choice for extrapolations to higher densities, as we are seeking mechanisms and conditions for avoiding a softening of the equation-of-state that would be in conflict with neutron star observations.

The computations of $U_\Lambda$ based on self-consistent solutions of the Bethe-Goldstone equation (\[eq:BGE\]) together with eq.(\[eq:U\]) involve a momentum cutoff not only in the input hyperon-nucleon potentials but also in the momentum-dependent single-particle potentials $U(k)$ that appear in energy denominators and are integrated over intermediate momenta $k$. We apply a cutoff by a factor $\exp[-(k/\lambda_{spp})^6]$ and choose $\lambda_{spp}$ between 600 and 650 MeV, close to the cutoff in the NLO hyperon-nucleon potential, while not altering the low-energy behavior. The ranges of these cutoffs, together with the two possible choices for the pair of constants ($H_1,H_2$), are reflected in the uncertainty estimates shown in the figures.

The computations of $U_\Lambda$ for both symmetric nuclear matter and pure neutron matter can be carried out safely (i.e. with well converging and numerically stable results) up to densities $\rho\simeq \rho_c=3.5\,\rho_0$. At even higher densities the occurence of numerical instabilities becomes prohibitive. One might argue that, in any case, performing calculations at $\rho\gtrsim 2\rho_0$ is not legitimate given the framework of chiral effective field theory with its limited range of applicability. On the other hand, the non-perturbative ladder summations of the Brueckner calculation lead beyond chiral perturbation theory, and the input interactions in the hyperon-nucleon sector are tested for total momenta $P$ up to and even above 800 MeV, whereas the neutron Fermi momentum even at central neutron star densities around $5\,\rho_0$ does not exceed 600 MeV, the typical cutoff scale in the calculations. We can therefore assume that at least qualitative extrapolations towards such high densities can be performed.

For extrapolations to densities beyond the technically (numerically) accessible $\rho_c \simeq 3.5\,\rho_0$, we argue as follows. Around $\rho_c$, one finds clear evidence in both symmetric nuclear matter and neutron matter that the density dependence of $U_\Lambda$ has a prominent quadratic behaviour, $U_\Lambda(\rho) \propto \rho^2$. This is expected as the three-body terms begin to dominate at high densities and grow as $\rho^2$. The two-body terms at high density increase linearly with $\rho$, and this expected subleading behaviour is also observed in the numerical results already at densities beyond $2\,\rho_0$. Therefore, assuming continuity an extrapolation of $U_\Lambda$ to high denstites is supposed to be well justified using the following ansatz for $\rho > \rho_0$: $$U_\Lambda(\rho) = u_0 + u_1\left({\rho\over\rho_0}-1\right) + u_2\left({\rho\over\rho_0}-1\right)^2\,. 
\label{eq:extrap}$$ The parameters $u_1$ and $u_2$ are determined by fits to the calculated single-particle potentials at densities below $\rho_c$. With $u_0$ fixed at $\rho=\rho_0$, it turns out that the high-density behaviour is indeed governed by the quadratic term, with $u_2 \gg u_1$.

![image](ULambda)

![image](NStarCP)

Results for the high-density extrapolation of $U_\Lambda$ in symmetric nuclear matter and neutron matter are displayed in fig.\[fig:ulambda\]. With two-body $YN$ interactions only, the linearly increasing $U_\Lambda$ at densities $\rho\gtrsim 2\rho_0$ is reminiscent of a Hartree potential. When increasingly repulsive three-body interactions are added, the leading $\rho^2$ dependence of the 3BF takes over at high densities. Note that in neutron matter, the strong repulsion in $U_\Lambda$ exceeds 100 MeV already at $\rho \simeq 3.5 \,\rho_0$ and reaches $U_\Lambda > 300$ MeV at $\rho\simeq 5\,\rho_0$, the density range characteristic of the inner core of neutron stars.

Hyperon and neutron chemical potentials in neutron stars
--------------------------------------------------------

If the chemical potentials of the $\Lambda$ hyperon and neutron fulfill the condition $\mu_\Lambda = \mu_n$ at some density in neutron star matter, it becomes energetically favourable to replace neutrons by $\Lambda$ hyperons via weak interactions. For a baryonic species $j$, the chemical potential $\mu_j$ is determined by the derivative of the energy density ${\cal E}$ with respect to the partial density $\rho_j$ of that species: $$\mu_j = {\partial{\cal E}\over\partial\rho_j}~.$$ Given a realistic neutron star equation-of-state, $P({\cal E})$, the Gibbs-Duhem relation for multicomponet sytems $$P({\cal E}) + {\cal E} = \rho_n\,\mu_n + \rho_p\,\mu_p + \rho_\Lambda\,\mu_\Lambda\,+ \dots~~.$$ determines the chemical potentials. In practice we use the microscopic EoS computed from a chiral nucleon-meson field theory in combination with functional renormalization group methods [@DW2017; @FW2019]. This EoS is consistent with all important nuclear physics constraints and neutron star observations. For simplicity we neglect the small (few percent) proton fraction from beta equilibrium which is part of the neutron star EoS, but its effect is marginal in the present context.

At the possible onset of $\Lambda$ hyperons, their kinetic energy vanishes as there is no $\Lambda$ Fermi sea to start with. Therefore the minimal $\Lambda$ chemical potential at that point is simply given by $$\mu_\Lambda(\rho) = M_\Lambda + U_\Lambda(\rho)~,
\label{eq:mulambda}$$

where $M_\Lambda$ is the $\Lambda$ hyperon mass and $\rho\simeq\rho_n$ is the neutron density. In the following the input potential $U_\Lambda(\rho)$ is the one presented in fig.\[fig:ulambda\] for neutron matter.

The comparison of $\mu_\Lambda$ and $\mu_n$ is shown in fig.\[fig:chempot\]. The uncertainty band of the neutron chemical potential is related primarily to the range of possible values of the nuclear symmetry energy, $E_\text{{sym}} = (32\pm 3)$ MeV. We note that this uncertainty band also includes $\mu_n$ as given in Ref.[@APR1998] for their maximally repulsive interaction (AV18+$\delta v$ + UIX\*).

Fig.\[fig:chempot\] points out that the combined repulsion from two- and three-body hyperon-nuclear interactions is indeed potentially strong enough to avoid the appearance of $\Lambda$ hyperons in neutron stars: one finds $\mu_\Lambda > \mu_n$ throughout the neutron star density range, with a set of three-body parameters that is consistent with basic constraints from hypernuclei. However, two-body $\Lambda N$ interactions alone would not be sufficient to suppress the occurrence of $\Lambda$ hyperons. Their onset would appear at densities around $3\,\rho_0$. The resulting softening of the EoS would be prohibitive.

The Brueckner calculations underlying eq.(\[eq:mulambda\]) are based on the $YN$ potential taken from ref.[@Haidenbauer2013] (referred to as NLO13) which features a strong $\Lambda N\leftrightarrow\Sigma N$ coupling. This channel coupling is a major source of the $\Lambda N$ attraction at intermediate distances. Its partial suppression by the Pauli principle acting on intermediate nucleon states in nuclear or neutron matter is responsible for the turnover from attraction to repulsion in the $\Lambda$-nuclear two-body force at a density well below $3\,\rho_0$ in nuclear matter. From fig.\[fig:ULambda1\] we also recall that the $\Lambda N$ attraction provided by NLO13 is slightly stronger than what is needed for the empirical hypernuclear shell-model potential. This sets the frame for constraining unknown pieces of the repulsive $YNN$ three-body force which in turn governs the behaviour of $\mu_\Lambda$ at high density. One might then ask whether alternative options exist which could imply different constraints and extrapolations for the 3BF sector. Such an alternative $YN$ potential (NLO19) has recently been studied [@Haidenbauer2019]. This potential is equivalent to NLO13 with respect to $\Lambda N$ and $\Sigma N$ scattering data, but it produces a stronger attractive $\Lambda$-nuclear potential than NLO13 at $\rho\simeq\rho_0$. This difference needs to be balanced by a more strongly repulsive three-body $\Lambda NN$ interaction. We have thus performed Brueckner calculations (for continuous choice) using NLO19, together with a 3BF readjusted to meet the hypernuclear $U_\Lambda(\rho=\rho_0) = -30$ MeV constraint. The effect at high density is to increase the $\Lambda$ chemical potential further in comparison with the NLO13 result plotted in fig.\[fig:chempot\], so as to widen the gap between $\mu_\Lambda$ and $\mu_n$ even more.

A further comment concerns $\Sigma$ hyperons in neutron star matter. At high density it is in principle possible to consider the condition $$\mu_{\Sigma^-} = \mu_n+\mu_e = 2\mu_n -\mu_p$$ for the occurrence of $\Sigma^-$ replacing an electron together with a neutron. However, this turns out not to be an option because the $\Sigma^-$ potential in neutron matter is strongly repulsive [@Petschauer2016b], so that $\mu_{\Sigma^-}$ remains separated by at least 50 MeV from $\mu_n + \mu_e$ at all relevant densities.

Summary and outlook {#sec:sum}
===================

Investigations of the properties of hyperons in nuclear and neutron matter have been extended to high baryon densities with the quest of exploring the possible occurrence of hyperons in neutron stars. Previous calculations based entirely on two-body hyperon-nucleon interactions suggested an onset for hyperons already at relatively low densities (between two and three times $\rho_0$). As a consequence the resulting equation-of-state is significantly softened and cannot support two-solar-mass neutron stars. The present work addresses the question whether hyperon-nucleon three-body forces (3BF) can provide the necessary repulsion to stop $\Lambda$ hyperons from replacing neutrons, maintaining the required stiffness of the neutron star EoS.

Calculations have been carried out within the self-consistent Brueckner-Hartree-Fock (BHF) approach using the continuous choice for intermediate spectra. The employed microscopic potentials for the hyperon-nucleon interaction are constructed from SU(3) chiral effective field theory at next-to-leading order. In addition, the pertinent N3LO nucleon-nucleon interaction is taken from chiral EFT.

Our main focus lies on the role of three-body forces, their reduction to density-dependent effective two-body interactions and the explicit treatment of $\Lambda NN\leftrightarrow\Sigma NN$ coupled channels in the BHF matrix equations. Unknown 3BF parameters have been constrained by reproducing the depth, $U_\Lambda(\rho=\rho_0) \simeq -30$ MeV, of the $\Lambda$ single-particle potential in nuclear matter as deduced from the phenomenology of hypernuclei.

The detailed discussion of $\Lambda$ and neutron chemical potentials in neutron matter then indicates that the $\Lambda NN$ three-body forces can indeed potentially develop sufficient repulsion to avoid the occurrence of $\Lambda$ hyperons in neutron stars, without the need to introduce strong extra assumptions in performing the extrapolation to high densities. A comparison using different versions of the chiral EFT-based $YN$ interactions (NLO13 versus NLO19) underlines this expectation further. While we regard these findings as a promising step towards a possible clarification of the so-called hyperon puzzle in neutron stars, we can nevertheless not yet claim that this problem has been “solved”. Detailed properties of the hyperon-nuclear three-body interactions must still be better constrained by data, e.g. in hypernuclear few-body systems. However, one can begin to consider the steadily increasing number of observational facts about neutron stars as part of a developing empirical data base that sets restrictive conditions for mechanisms involving strangeness in dense baryonic matter.\
\
\
\[app:angularaverage\]

Here, we present details of the angular averaging procedure that preceeds the derivation of eq.(\[eq:BGE\]). Total and relative momenta for two baryons $B_1$ and $B_2$ are introduced by $$\vec P = \vec p_1 + \vec p_2\,,\quad
\vec k = \frac{\xi_{12}\vec p_1-\vec p_2}{1+\xi_{12}}\,,\quad
\xi_{12} = \frac{M_{2}}{M_{1}}\,.$$

A standard approximation is applied to replace $Q/e$ by the ratio of angle-averaged quantities, $\bar Q/\bar e$. The Pauli operator, averaged over the angle $\theta$ between $\vec P$ and $\vec k$, is given by $$\begin{aligned}
\bar Q_{\nu}(P,k) & =
\frac12\int_{-1}^1\!\mathrm d\cos\theta\
\Theta\big(|\vec p_1| - k_F^{(1)}\big)\Theta\big(|\vec p_2| - k_F^{(2)}\big) \\
 & = [0|\frac{[-1|z_1|1]+[-1|z_2|1]}2|1] \,,
{\addtocounter{equation}{1}\tag{\theequation}}\end{aligned}$$ with $$\begin{aligned}
z_1 & = \frac{1+\xi_{12}}{2kP}\bigg\{\left(\frac1{1+\xi_{12}}P\right)^2+k^2-(k_F^{(1)})^2\bigg\}\,, \\
z_2 & = \frac{1+1/\xi_{12}}{2kP} \bigg\{\left(\frac{\xi_{12}}{1+\xi_{12}}P\right)^2 + k^2-(k_F^{(2)})^2\bigg\} \,.
{\addtocounter{equation}{1}\tag{\theequation}}\end{aligned}$$ and the notation $[a|b|c]\equiv\operatorname{max}(a,\operatorname{min}(b,c))$ introduced in ref.[@Schulze1998]. The angle-averaged energy denominator is given by $$\begin{aligned}
\bar e_{\nu}(P,k;\omega) = \omega & - \frac{P^2}{2M_\nu}-\frac{k^2}{2\mu_\nu} -M_\nu \\
& - {\operatorname{Re}}U_{B_1}(\bar p_1)- {\operatorname{Re}}U_{B_2}(\bar p_2) \,,
{\addtocounter{equation}{1}\tag{\theequation}}\end{aligned}$$ with $M_{\nu}=M_1+M_2$ and $\mu_{\nu}=M_1 M_2/(M_1+M_2)$. The angle-average is done for the arguments of the single particle potentials $U_{B_i}$ of the intermediate baryons: $$\begin{aligned}
\bar p_1 &= \left(\tfrac1{(1+\xi_{12})^2} P^2 + k^2 + 2\tfrac1{1+\xi_{12}} P k\, \overline{\cos\theta}\right)^{1/2} \,, \\
\bar p_2 &= \left(\tfrac{\xi^2_{12}}{(1+\xi_{12})^2} P^2 + k^2 - 2\tfrac{\xi_{12}}{1+\xi_{12}} P k\, \overline{\cos\theta}\right)^{1/2} \,,
{\addtocounter{equation}{1}\tag{\theequation}}\end{aligned}$$ with $$\begin{aligned}
\overline{\cos\theta} & = \frac{ \int_{-1}^1\!\mathrm d\cos\theta\, \cos\theta\, Q(\vec P,\vec k) }{ \int_{-1}^1\!\mathrm d\cos\theta\, Q(\vec P,\vec k) } \\
& = \frac12\big([-1|z_2|1]-[-1|z_1|1]\big) \,.
{\addtocounter{equation}{1}\tag{\theequation}}\end{aligned}$$ where $Q(\vec P,\vec k)$ is the exact Pauli blocking operator.

It is common practice to introduce a further simplification, replacing the squared momenta $P^2=P^2(\vec p_1, \vec k)$ and $p^2_2=p^2_2(\vec p_1, \vec k)$ entering the Bethe-Goldstone equation by their angle averages: $$\begin{aligned}
\bar{P}^2(p_1,k) & = \frac{\int_{|\vec p_2|\leq k_F^{(2)}}\!\mathrm d\cos\vartheta\, P^2(p_1,k,\cos\vartheta)}{\int_{|\vec p_2|\leq k_F^{(2)}}\!\mathrm d\cos\vartheta} \\
&=
(1+\xi_{12})^2\left[p_1^2+k^2-p_1k(1+[-1|x_0|1])\right] \,, \\
\bar p_2^2(p_1,k)
&=
\frac{\xi_{12}}{1+\xi_{12}}\bar P^2(p_1,k) + (1+\xi_{12})k^2 -\xi_{12}\,p_1^2 \,,
{\addtocounter{equation}{1}\tag{\theequation}}\end{aligned}$$ where $\vartheta$ is the angle between $\vec p_1$ and $\vec k$, and one finds $$x_0 = \frac{\xi_{12}^2p_1^2+(1+\xi_{12})^2k^2-(k_F^{(2)})^2}{2\xi_{12}(1+\xi_{12})p_1k} \,.
{\addtocounter{equation}{1}\tag{\theequation}}$$ Note that baryon $B_2$ in the initial state is within its Fermi sea.

The weight function $W(p,k)$ that appears in the calculation of the single-particle potential $U_B(p)$ via eq.(\[eq:U\]) is given by $$\begin{aligned}
W(p,k) &= \frac1{4\pi}\int_{|\vec p_2|\leq k_F^{(2)}}\!\mathrm d\Omega_k = \frac12(1-[-1|x_0|1]) \,.
{\addtocounter{equation}{1}\tag{\theequation}}\end{aligned}$$ The integration boundaries, $k_\mathrm{min}$ and $k_\mathrm{max}$, of the relative momentum are determined by the condition $W(k_1,k) = 0$, which leads to $$\begin{aligned}
k_\mathrm{min} &= \operatorname{max}\left( 0, \frac{-k_F^{(2)}+\xi_{12} p}{1+\xi_{12}} \right) \,,\
k_\mathrm{max} = \frac{k_F^{(2)}+\xi_{12} p}{1+\xi_{12}} \,.
{\addtocounter{equation}{1}\tag{\theequation}}\end{aligned}$$

\
Eqs. (\[eq:BGE\],\[eq:U\]) are solved numerically by alternately iterating both equations until the potentials $U_B$ converge. For certain kinematical conditions, e.g. involving large momenta in the initial state, the energy denominator in eq.(\[eq:BGE\]) vanishes, giving rise to a pole in the Bethe-Gold- stone equation. To make the integral numerically manageable, the principal value prescription is utilized: $$-\!\!\!\!\!\!\int_0^\infty \!dk \frac{N(k)}{D(k)} = 
    \int_0^\infty \!dk \left(\frac{N(k)}{D(k)}-\frac{N(k_0)}{D^\prime(k_0)} \frac{2k_0}{k^2-k_0^2}\right)\,,$$ for a simple pole at $k_0$, i.e. $D(k_0)=0$ and $D'(k_0)\neq0$. This prescription eventually meets its limits at densities beyond $3.5\,\rho_0$, once second-order poles begin to show up.

For densities $\rho > 2\rho_0$ in pure neutron matter, slow convergence of $U_B$ occurs because of increasingly strong input potentials. Convergence is improved by averaging $U_B(k)$ over subsequent iterations. High-density calculations occasionally require intervention by estimating an appropriate starting point for iterations from lower-density results. However, numerical stability is generally not guaranteed any more for $\rho > \rho_c \simeq 3.5\,\rho_0$. For densities $\rho > \rho_c$ we use the power series extrapolation in eq.(\[eq:extrap\]), fitted to numerically stable results at $\rho \lesssim \rho_c$ as described in the text.

As discussed in ref. [@Kohno2018], the cutoff $\lambda$ in the chiral input potential causes numerical oscillations of $U_B(k)$ for momenta $k > \lambda$. These oscillations tend to slow down the convergence of $U_B$ and are of no physical relevance. They are suppressed by the additional cutoff factor $\exp[-(k/\lambda_{spp})^6]$ attached to the single-particle potentials as described in Section \[subsec:spp\].\

[99]{}

P. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, Nature [**467**]{}, 1081 (2010).

J. Antoniadis [*et al.*]{}, Science [**340**]{}, 6131 (2013).

E. Fonseca [*et al.*]{}, Astroph. J. [**832**]{}, 167 (2016).

H.T. Cromartie [*et al.*]{}, Nat. Astron. [**4**]{}, 72 (2019).

B.P. Abbott [*et al.*]{} (LIGO & Virgo Collaborations), Phys. Rev. Lett. [**119**]{}, 161101 (2017).

E.R. Most, L.R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. [**120**]{}, 261103 (2018).

S. De, D. Finstad, J.M. Lattimer, D.A. Brown, E. Berger, and C.M. Biwer, Phys. Rev Lett. [**121**]{}, 091102 (2018).

A. Akmal, V.R. Pandharipande, and D.G. Ravenhall, Phys. Rev. [**C58**]{}, 1804 (1998).

M. Drews and W. Weise, Phys. Rev. [**C91**]{}, 035802 (2015).

M. Drews, and W. Weise, Prog. Part. Nucl. Phys. [**93**]{}, 69 (2017).

B. Friman and W. Weise, Phys. Rev. [**C 100**]{}, 065807 (2019).

G. Baym, T. Hatsuda, T. Kojo, P.D. Powell, Y. Song, and T. Takatsuka, Rept. Prog. Phys. [**81**]{}, 056902 (2018).

G. Baym, S. Furusawa, T. Hatsuda, T. Kojo, and H. Tagashi, Astrophys. J. [**885**]{}, 42 (2019).

L. McLerran and S. Reddy, Phys. Rev. Lett. [**122**]{}, 122701 (2019).

Y. Fujimoto, K. Fukushima, and K. Murase, arXiv:1903.03400 \[nucl-th\].

H. Djapo, B.J. Schaefer, and J. Wambach, Phys. Rev. [**C 81**]{}, 035803 (2010).

T. Hell and W. Weise, Phys. Rev. [**C 90**]{}, 045801 (2014).

D. Lonardoni, A. Lovato, S. Gandolfi, and F. Pederiva, Phys. Rev. Lett. [**114**]{}, 092301 (2015).

P.F. Bedaque and A.W. Steiner, Phys. Rev. [**C 92**]{}, 025803 (2015).

J. Haidenbauer, S. Petschauer, N. Kaiser, U.-G. Mei[ß]{}ner, A. Nogga, and W. Weise, Nucl. Phys. [**A 915**]{}, 24 (2013).

S. Petschauer, N. Kaiser, J. Haidenbauer, U.-G. Mei[ß]{}ner, and W. Weise, Phys. Rev. [**C 93**]{}, 014001 (2016).

S. Petschauer, J. Haidenbauer, N. Kaiser, U.-G. Mei[ß]{}ner, and W. Weise, Nucl. Phys. [**A 957**]{}, 347 (2017).

M. Kohno, Phys. Rev. [**C 97**]{}, 035206 (2018).

S. Petschauer, J. Haidenbauer, N. Kaiser, U.-G. Mei[ß]{}ner, and W. Weise, Eur. Phys. J. [**A 52**]{}, 15 (2016).

J. Haidenbauer, U.-G. Mei[ß]{}ner, N. Kaiser, and W. Weise, Eur. Phys. J. [**A 53**]{}, 121 (2017).

J. Haidenbauer and U.-G. Mei[ß]{}ner, Phys. Rev. [**C 72**]{}, 044005 (2005).

T.A. Rijken, V.G.J. Stoks, and Y. Yamamoto, Phys. Rev. [**C 59**]{}, 21 (1999).

J.W. Holt, N. Kaiser, and W. Weise, Phys. Rev. [**C 81**]{}, 024002 (2010).

B. Day, Rev. Mod. Phys. [**39**]{}, 719 (1967).

H.-J. Schulze, M. Baldo, U. Lombardo, J. Cugnon, and A. Lejeune, Phys. Rev [**C 57**]{}, 704 (1998).

T.A. Rijken, V.G.J. Stoks, and Y. Yamamoto, Phys. Rev [**C 59**]{}, 21 (1999).

I. Vidana, A. Polls, A. Ramos, M. Horth-Jensen, and V.G.J. Stoks, Phys. Rev [**C 61**]{}, 025802 (2000).

M. Kohno, Y. Fujiwara, T. Fujita, C. Nakamoto, and Y. Suzuki, Nucl. Phys. [**185**]{}, 229 (2000).

D.R. Entem and R. Machleidt, Phys. Rev. [**C 68**]{}, 041001 (2001).

A. Gal, E.V. Hungerford, and D.J. Millener, Rev. Mod. Phys. [**88**]{}, 035004 (2016).

J. Haidenbauer, U.-G. Mei[ß]{}ner, and A. Nogga, arXiv:1906.11681 \[nucl-th\].

[^1]: We remind that due to direct and exchange contributions to the density-dependent effective $YN$ potentials, these have to be weighted with a statistical factor $1/2$ when added to the two-body potentials [@Kohno2018].
