# Baryon and meson masses in the Nambu–Jona-Lasinio model: a Bayesian approach

Antoine Pfaff<EMAIL_ADDRESS>Univ. Lyon, Univ. Claude Bernard Lyon 1,
CNRS/IN2P3, Institut de Physique des 2 Infinis de Lyon, UMR 5822, 69622
Villeurbanne, France  SUBATECH, Nantes University, IMT Atlantique, IN2P3/CNRS,
4 rue Alfred Kastler, 44307 Nantes cedex 3, France  Hubert Hansen  Univ. Lyon,
Univ. Claude Bernard Lyon 1, CNRS/IN2P3, Institut de Physique des 2 Infinis de
Lyon, UMR 5822, 69622 Villeurbanne, France  Joerg Aichelin  SUBATECH, Nantes
University, IMT Atlantique, IN2P3/CNRS,
4 rue Alfred Kastler, 44307 Nantes cedex 3, France  Juan M. Torres-Rincon
<EMAIL_ADDRESS>Departament de Física Quàntica i Astrofísica and Institut de
Ciències del Cosmos (ICCUB), Facultat de Física, Universitat de Barcelona,
Martí i Franquès 1, 08028 Barcelona, Spain

###### Abstract

We investigate the capabilities of the Nambu–Jona-Lasinio model to describe
and reproduce fundamental vacuum properties of Quantum Chromodynamics, notably
the hadronic spectrum. Mesons are described as quark-antiquark bound states at
the level of the random phase approximation of the Bethe-Salpeter equation,
while baryons are characterized as quark-diquark bound states within the
static approximation of the Fadeev equation. Within a Bayesian framework, we
constrain the model by phenomenologically known quantities and study the
implications on its parameters and predictions in vacuum, as well as the
correlations between the two. We find that within our framework, the vacuum
masses of mesons and baryons can be reasonably well reproduced. Scalar
diquarks need to be significantly bound in order to correctly reproduce the
masses of the baryon octet, therefore enforcing values of the scalar diquark
coupling larger than what is suggested by the canonical Fierz values. These
findings could have important implications on the phenomenology of strongly-
interacting matter at high temperature and density as well as of compact star
physics.

††preprint: APS/123-QED

##  I Introduction

Modeling of hadrons—and their possible composite nature—is a challenging task
because their internal structure lies out of the realm of perturbative Quantum
Chromodynamics (QCD). In addition, they are confined color neutral objects and
confinement is one of the least understood properties of QCD. Hadrons are,
however, in the center of interest of current research on physical systems
composed by strongly-interacting particles. Besides, these systems can be
exposed to a finite temperature—like the medium formed in relativistic heavy-
ion collisions—or to finite density—such as the nuclear medium, or the
interior of neutron stars. The study of hadrons in matter is therefore in the
center of interest of present nuclear physics.

An attractive possibility is to solve the QCD field equations on a discretized
lattice using a fully numerical approach. This idea has been very successful
to extract hadron properties, both in vacuum [1, 2, 3, 4] and at finite
temperature [5, 6, 7]. The latest lattice-QCD (LQCD) calculations for the
light sector can in most cases approach the “physical point” (with a physical
pion mass ${m_{\pi}\approx 138}~{}$MeV), and the extrapolation to the
continuum limit can address the spectroscopy in a rather faithful way [8, 9,
10]. However, calculations of the baryon spectrum reaching the physical point
have not been implemented yet. In addition, LQCD calculations are
computationally costly, and cannot be straightforwardly extended to finite
baryon densities due to the so-called “numerical sign problem”.

Another option to address hadron spectroscopy is the application of
(truncated) QCD Dyson-Schwinger and Bethe-Salpeter equations which can
accommodate both finite temperature and density (see Ref. [11] and references
therein). Results of this approach on vacuum spectroscopy are presented in the
review [12], while at finite temperature some results for light mesons have
been recently published [13, 14].

In this work we take the approach of effective field theories. In particular,
we would like to describe a substantial set of hadrons—including mesons and
baryons—but incorporating the minimal set of degrees of freedom and parameters
in the effective model. We also would eventually like to extend the
calculations to finite temperature and density (while in this work we restrict
ourselves to the vacuum case). With this goal in mind, we focus on the
Nambu–Jona-Lasinio (NJL) model: an effective approach of QCD at low energies
describing interactions among quarks while gluons are integrated out. The
effective interactions result (in the simplest local version of the model) in
multiquark contact terms with energy-independent strength. This interaction
can be tuned to reproduce the physical value of several observables in the
low-energy hadron sector.

The NJL model has been profusely used in the literature—both in vacuum and in
a thermal/dense medium—due to its technical simplicity and its properties [15,
16, 17, 18, 19, 20]: it describes the spontaneous symmetry breaking property
of QCD and the generation of (pseudo) Goldstone bosons, it can easily be
extended to include in-medium properties, it illustrates the chiral symmetry
restoration at finite temperatures/densities, it can effectively accommodate
the effects of the U(1)A axial anomaly at the meson level
($\eta-\eta^{\prime}$ mass splitting) and, using the Bethe-Salpeter or the
Fadeev equations, it can describe the dynamical generation of colorless mesons
and baryons, which can be identified with physical hadrons.

With respect to the last point, the local NJL model is able to generate the
masses of ground state mesons (scalar, pseudoscalar, vector and axial vectors)
as bound states in the quark-antiquark scattering, as well as baryons (members
of the SU$(3)_{f}$ octet and decuplet) as bound states of a quark and a
diquark (the latter itself being a bound state of two quarks). The original
formalism is based on the relativistic Faddeev equation for three quarks
which, under certain approximations, can be simplified to a two-body equation
[21, 22]. Using the NJL model, some vacuum baryon masses were already obtained
in [22, 23, 24], and more recently, as functions of temperature and/or density
in [25, 26, 27, 28, 29]—notice that at finite temperature/density one can also
use the Polyakov loop-extended NJL (PNJL) model [30, 31], which incorporates
certain gluonic components, allowing to suppress part of the unphysical quark
modes in the partition function at low $T$ and $n$.

In a previous publication [29], some of us studied the thermal evolution of
several baryons with three flavors with the NJL and PNJL models. We were able
to satisfactorily reproduce the vacuum masses of 8 baryon states
($N,\Delta,\Sigma,\Lambda,\Xi,\Sigma^{\star},\Xi^{\star},\Omega$) with only
two additional parameters in the NJL model (scalar and axial-vector couplings
within the diquark sector). This was achieved without spoiling the good
agreement of meson masses, but at the cost of having a sizable light quark
dressed mass of $m_{q}\simeq 480$ MeV. Being not an observable, such a value
for the quark mass is not a particularly troubling result—a typical quark mass
would be of the order of $m_{\Delta}/3\simeq 400$ MeV, as it seems to be
required in previous calculations [22, 23, 24]—but it was unclear whether this
was a necessary condition of the model in order to fit the meson and baryon
spectroscopy, or whether a different parametrization with a smaller quark mass
could also be found.

In this paper we employ the Bayesian formalism to perform a systematic
analysis of the NJL model prediction for the vacuum hadron spectrum. The
Bayesian approach has notably been implemented with the NJL model in the
context of the study of the compact star equation of state, where the models
are conditioned by measurements of the star masses, radii and tidal
deformabilities [32, 33, 34, 35]. However, in the current work the
parametrization of the model is instead constrained by a likelihood informed
with the experimentally measured masses of the full pseudoscalar meson octet,
baryon octet and baryon decuplet as well as the pion and kaon weak decay
constants. Our goals in this work are: (i) Find whether the minimal version of
the NJL model (local, mean field and with the minimal number of couplings) is
able to provide a good matching of these quantities overall in a systematic
approach, (ii) obtain and analyze the actual correlations between parameters,
and between parameters and observables, (iii) address the question whether a
good fit of the hadron spectrum in the NJL model requires to have a relatively
massive quark, as indicated in the recent studies [27, 36, 29], and (iv)
sketch an “optimal” parametrization which can be used in future studies of
hadron spectrum and the phase diagram of the model at finite temperature and
density.

This paper is organized as follows: In Sec. II we introduce our minimal
version of the NJL model and detail the level of approximations employed to
compute both meson and baryon vacuum masses. In Sec. III we describe the
Bayesian method that we apply to explore the full parameter space and confront
the outcome of the model to the experimental results. Finally, in Sec. IV we
analyze the posterior distributions of the parameters and physical observables
(with a particular emphasis on the diquark and baryon sector), and comment the
most interesting correlations found. We present a summary of this work in Sec.
V, and conclude the paper with a set of appendices to which all the technical
details are referred.

##  II Nambu–Jona-Lasinio Model

Our analysis is based on the description of the low-energy hadronic spectrum
of the three-flavor Nambu–Jona-Lasinio (NJL) quark model [15, 16, 17, 18, 19,
20], with the following Lagrangian density,

$\displaystyle\mathcal{L}_{\textrm{NJL}}$
$\displaystyle=\overline{q}(i\not{\partial}-\hat{m}_{0})q+G\
[(\overline{q}\tau^{a}q)^{2}+(\overline{q}i\gamma_{5}\tau^{a}q)^{2}]$
$\displaystyle-K\Big{(}\operatorname{det}\limits_{f}[\overline{q}(1+\gamma_{5})q]+\operatorname{det}\limits_{f}[\overline{q}(1-\gamma_{5})q]\Big{)}$
$\displaystyle+G_{d}\
(\overline{q}i\gamma_{5}C\tau^{A}\lambda^{A^{\prime}}\overline{q}^{T})(q^{T}Ci\gamma_{5}\tau^{A}\lambda^{A^{\prime}}q)$
$\displaystyle-G_{d,V}\
(\overline{q}\gamma^{\mu}C\tau^{S}\lambda^{A^{\prime}}\overline{q}^{T})(q^{T}C\gamma_{\mu}\tau^{S}\lambda^{A^{\prime}}q)\
.$ (1)

In Eq. (1) we denote $\tau^{a}$, ${a=0,1,...,8}$ the $\mathfrak{su}(3)$ Gell-
Mann matrices in flavor space (with
${\tau^{0}=\sqrt{\frac{2}{3}}\mathbbb{1}_{f}}$) and as $\lambda^{a^{\prime}}$,
${a^{\prime}=0,1,...,8}$ the $\mathfrak{su}(3)$ matrices in color space (with
$\lambda^{0}=\sqrt{\frac{2}{3}}\mathbbb{1}_{c}$). $A$ and $A^{\prime}$ are
summed over the indices of antisymmetric Gell-Mann matrices;
$A,A^{^{\prime}}=2,5,7$; while $S$ is summed over the indices of symmetric
Gell-Mann matrices; ${S=0,1,3,4,6,8}$. In the diquark sector—last two lines of
Eq. (1)—the charge conjugation matrix is given by $C=i\gamma_{0}\gamma_{2}$.

The quarks are described by the flavor triplet spinor field ${q=(u,d,s)}$ with
bare masses ${\hat{m}_{0}=\operatorname{diag}(m_{0,u},m_{0,d},m_{0,s})}$.
Quarks interact with each other through four- and six-fermion contact vertices
controlled by four coupling constants: the scalar-pseudoscalar coupling $G$,
the ’t Hooft coupling $K$, the scalar diquark coupling $G_{d}$, and the axial
diquark coupling $G_{d,V}$.

The NJL Lagrangian shown in Eq. (1) obeys the global
SU(3)${}_{V}\times$SU(3)${}_{A}\times$U(1)B symmetry of QCD in the limit of
vanishing bare quark masses ($\hat{m}_{0}=0$), and substitutes its local color
gauge symmetry with a global invariance under SU$(3)_{c}$ transformations. The
so-called axial anomaly of QCD is accounted for by the presence of the
6-fermion Kobayashi–Maskawa–’t Hooft interaction term, which breaks the U(1)A
symmetry explicitly [37, 38, 39]. The quark-antiquark vector channel is
neglected as it does not have any effect in the vacuum at the mean-field level
and we do not investigate the vector meson spectrum. Although other diquark
channels are in principle needed in order to meet the symmetry requirements,
they are not considered in this work as they do not contribute to the baryons
in our scheme.

The form of the four-fermion part of the NJL Lagrangian can be directly
derived by carrying out a Fierz transformation of a global color current-
current interaction, which can be seen as an approximation of QCD in the limit
of low gluon momentum exchange [20]. Following this procedure, four-fermion
interaction terms in all possible spin, flavor and color channels emerge, and
we only kept in Eq. (1) the ones that give rise to the low-mass hadronic
spectrum of interest. Since the original theory of QCD has only one coupling
constant (the strong coupling $g_{s}$), the Fierz method yields fundamental
relationships between the coupling constants in each interacting channel [29,
15, 20]. For the terms appearing in Eq. (1) and for $N_{c}=3$, the relevant
Fierz relations are $G_{d}=3G/4$ and $G_{d,V}=3G/8$ [20]. In this work
however, we will forego these relations and keep $G_{d}$ and $G_{d,V}$ as free
parameters. As we shall see below, the above relations are not compatible with
the experimental masses of the baryon spectrum in our framework.

From the Lagrangian of Eq. (1), the dressed mass $m_{f}$ of the quark of
flavor $f$ can be obtained by solving the self-consistent Dyson equation [40]
for the quark propagator, which in the Hartree approximation eventually
yields:

$m_{i}=m_{i,0}-4G\langle\overline{q}_{i}q_{i}\rangle+2K\langle\overline{q}_{j}q_{j}\rangle\langle\overline{q}_{k}q_{k}\rangle\
,$ (2)

where ($i$,$j$,$k$) denotes any permutation of ($u$,$d$,$s$). The chiral
condensates $\langle\overline{q}_{f}q_{f}\rangle$ are related themselves to
the quark masses by the single line loop integral:

$\langle\overline{q}_{f}q_{f}\rangle=-i\int\frac{d^{4}k}{(2\pi)^{4}}\operatorname{Tr}S_{f}(k)\
,$ (3)

where the traces are performed in color and spinor spaces, and $S_{f}(p)$ is
the dressed quark propagator, given at zero temperature by

$S_{f}(k)=\frac{\not{k}+m_{f}}{k^{2}-m_{f}^{2}+i\varepsilon}$ (4)

After performing the traces, one simply obtains

$\begin{split}\langle\overline{q}_{f}q_{f}\rangle&=-4im_{f}N_{c}\int^{\Lambda}\frac{d^{4}k}{(2\pi)^{4}}\frac{1}{k^{2}-m_{f}^{2}+i\varepsilon}\\\
&=-4im_{f}N_{c}I_{1}(m_{f},\Lambda)\end{split}\ ,$ (5)

Since the integral $I_{1}$ is divergent, it must be regularized by the
introduction of an ultraviolet cutoff $\Lambda$, whose scale is related to the
asymptotic freedom scale of QCD. In this work, we use a simple 3-momentum
noncovariant regularization scheme in which the cutoff $\Lambda$ is imposed
after carrying out the $k_{0}$ integration. Details on the evaluation of the
integrals are relegated to the App. A.

Mesons and diquarks can be described within the NJL model by solving the
Bethe-Salpeter equation for the scattering of a quark and an antiquark (or two
quarks in the case of diquarks). The existence of a bound state can then be
inferred from the $T$-matrix. In the so-called ring approximation (random
phase approximation) for the $T$-matrix, the infinite resummation of diagrams
is reduced to a geometric sum involving the quartic coupling and the two-
particle loop function $\Pi_{X}$(p),

$\displaystyle T_{X}(p)$
$\displaystyle=\mathcal{G}_{X}+\mathcal{G}_{X}\Pi_{X}(p)\mathcal{G}_{X}+...$
$\displaystyle=\frac{\mathcal{G}_{X}}{1-\mathcal{G}_{X}\Pi_{X}(p)}\ ,$ (6)

where $\mathcal{G}_{X}=2G_{X}$ is, up to a loop factor 2, the four-quark
coupling constant appearing in the Lagrangian, that selects the appropriate
interaction channel associated to the quantum numbers $X$. $T_{X}$ refers to
the $T$-matrix in the channel $X$, which can describe a bound state if a pole
in the physical Riemann sheet of $p$ is generated. The mass of the bound state
can be inferred as the pole position in the real energy axis. Therefore, we
identify the mass $m_{X}$ (at rest) of the meson/diquark $X$ as the solution
of the equation,

$1-\mathcal{G}_{X}\Pi_{X}(p)|_{p_{0}=m_{X},\,\bm{p}=0}=0\ .$ (7)

While this pole can in principle lie outside of the real axis, reflecting a
finite width related to a decay into a quark-antiquark pair, we will consider
cases where $X$ is stable, such that its mass is indeed real. Each bound state
$X$ can be associated to an appropriate quark-antiquark (or quark-quark in the
case of diquarks) interaction channel $\Gamma_{X}$ depending on the quantum
numbers exchanged through the interaction vertex (e.g. the pseudoscalar pion
$\pi^{0}$ is associated to the vertex
${\Gamma_{\pi^{0}}=i\gamma_{5}\tau_{3}\mathbbb{1}_{c}}$). The expression of
the loop function $\Pi_{X}(p)$ can be obtained using the Feynman rules for
fermions, and is given by

$\Pi_{X}(p)=i\int^{\Lambda}\frac{d^{4}k}{(2\pi)^{4}}\operatorname{Tr}\Big{(}\Gamma_{X}S(p+k)\Gamma_{X}^{\dagger}S(k)\Big{)}\
,$ (8)

where ${S(k)=\operatorname{diag}(S_{u}(k),S_{d}(k),S_{s}(k))}$ and the trace
must be performed in color, flavor and spinor spaces. Further details on the
calculation of this loop function in the meson and diquark channels of
interest are given in App. B.

For each bound state, one can additionally compute its coupling strength
$g_{X\rightarrow qq^{\prime}}$ corresponding to the residue of the propagator
of bound state $X$ at its pole

$g_{X\rightarrow qq^{\prime}}^{2}=\Big{(}\frac{\partial\Pi_{X}}{\partial
p^{2}}\Big{|}_{p^{2}=m_{X}^{2}}\Big{)}^{-1}\ .$ (9)

Around its pole, the $T$-matrix element of the bound state $X$ can then be
expressed in the form

$T_{X}(p)\approx-\frac{g_{X\rightarrow qq^{\prime}}^{2}}{p^{2}-m_{X}^{2}}\ ,$
(10)

which can readily be interpreted as the bound state propagator. Note that the
equations above apply strictly speaking only for spin 0 bound states. In the
case of spin 1 (in particular for the axial diquarks), a similar procedure can
be applied in order to take into account the Lorentz structure of the
propagator [29, 17].

Finally, the axial decay constants $f_{X}$ can be obtained from the evaluation
of the transition matrix element between a meson and the vacuum.

$if_{X}p^{\mu}=\langle 0|j^{\mu}_{5,X}|X\rangle\ .$ (11)

In the NJL model, the evaluation of the diagram at one-loop level gives [15]

$if_{X}p^{\mu}=i\frac{g_{X\rightarrow
qq^{\prime}}}{2}\int^{\Lambda}\frac{d^{4}k}{(2\pi)^{4}}\operatorname{Tr}\Big{(}\gamma^{\mu}\Gamma_{X}S(p+k)\Gamma_{X}^{\dagger}S(k)\Big{)}\
,$ (12)

which is computed as explained in App. C. In this work, we will only focus on
the pion and kaon decay constants which are experimentally well-known
quantities.

In an analogous way to the modelling of mesons as quark-antiquark bound
states, baryons can be constructed by solving the so-called Faddeev equation
[41] for the scattering of three quarks. While even in the NJL model this
equation is difficult to solve, the quark-diquark framework provides an
appealing approximation in order to simplify the problem (see [12] and
references therein). In this case, the three body scattering problem can be
simplified to a two body problem with a quark and a diquark, in which case the
Faddeev equation reduces to a Bethe-Salpeter like equation, and the
interaction between a quark and a diquark is brought down to a simple quark
exchange.

To simplify further the calculation, we apply the “static approximation”
consisting in neglecting the momentum carried by the exchanged quarks in the
diagram summation [22]. The calculation is then reduced to a geometric sum
very similar to Eq. (6). Such an approximation is valid as long as the in-
medium mass of the quark remains high and the chiral symmetry remains
spontaneously broken. At high temperature and density, chiral symmetry should
be restored and the static approximation is not justified anymore. However,
since in this work we are only interested in studying chiral properties in the
vacuum, this approximation should remain valid. We refer to Ref. [29] and App.
D for further details on the calculation.

Since diquarks are colored objects, each baryon wavefunction must be
decomposed with a specific superposition of quark and diquark combinations in
order to ensure that the constructed baryons are colorless. Baryons of the
octet (with $J^{\pi}=\frac{1}{2}^{+}$) are assumed to be constituted only of
scalar diquarks, while baryons of the decuplet (with
$J^{\pi}=\frac{3}{2}^{+}$) are assumed to be constituted of axial diquarks.
The corresponding projection for each baryon can be found in Ref. [29, 42],
following the principle that the total baryon wave function should be
antisymmetric in its decomposition.

##  III Methodology

The purpose of this work is to investigate the parameter dependence of the NJL
model predictions, as well as to understand how to choose these parameters in
order to be in agreement with experimentally measured properties of the
hadronic spectrum. To do so, we adopt a Bayesian methodology in order to
constrain Monte Carlo generated models to satisfy the physical constraints. A
large number of NJL parameter sets is generated within a large and uninformed
(uniform) prior. The predictions of each set are then computed and confronted
with experimental data (here observables of the SU$(3)_{f}$ hadronic spectrum)
to select parameter sets that are consistent with real world physics. All
calculations are performed in vacuum.

Our framework based on the NJL model comprises a total of seven different
parameters; two bare masses for the light and strange quarks (we work in the
isospin limit where $m_{0}\equiv m_{0,u}=m_{0,d}$), four coupling constants
($G$, $K$, $G_{d}$ and $G_{d,V}$), and one momentum cutoff $\Lambda$. In other
works using a similar framework [43, 44, 28, 45, 46, 29, 27], these parameters
were partially or totally fine tuned in order to reproduce several
phenomenological quantities of the hadronic spectrum. However, this method
possesses some drawbacks:

  1. (i)

there are more than seven quantities that can be fitted from hadron
phenomenology, meaning that the choice of the fitted quantities can affect the
results significantly and that the resulting parameter set is not necessarily
optimal for reproducing known experimental quantities,

  2. (ii)

this method does not reflect the uncertainties (both experimental and
theoretical) that can blur the model inputs.

  3. (iii)

this method cannot account for the correlations between the model parameters
and the observables, which carry significant information on the relationships
between physical sectors of the model.

For these reasons, we propose a method to consistently take into account the
uncertainties on all available quantities.

###  III.1 Bayesian statistical analysis

The core of a Bayesian statistical analysis lies on Bayes’ theorem in
probability theory [47]. Assuming a set of observables $\\{o_{i}\\}$
constrained by some data, the theorem states

$\mathcal{P}(\\{o_{i}\\}|{\rm data})=\mathcal{N}\ \mathcal{P}({\rm
data}|\\{o_{i}\\})\ \mathcal{P}(\\{o_{i}\\})\ .$ (13)

It relates the posterior probability of a model $\mathcal{P}(\\{o_{i}\\}|{\rm
data})$ to both its prior probability $\mathcal{P}(\\{o_{i}\\})$, and to the
likelihood function $\mathcal{P}({\rm data}|\\{o_{i}\\})$ . The quantity
$\mathcal{N}$ is a normalization factor also known as evidence. Here, each
individual parameter set (which will be also denote as “individual model”, as
it is the accepted term in this framework), is represented by the collection
of predictions it makes on the observables $\\{o_{i}\\}$.

The prior probability represents the initial knowledge that we have on the
observables $\\{o_{i}\\}$. It is entirely determined by the modeling of the
observables and the amount of freedom that is allowed within this modeling. In
order to perform a study that is the most agnostic possible and span the
largest uncertainties, the prior distribution should remain as uninformative
as possible.

The likelihood estimates the compatibility of the model with the data by
attributing a probability to the experimental results (the data) for each
possible underlying model. In this sense, it acts as a filter on the different
models considered, by performing a selection to give a preference to the
models that reproduce best the data. Its effects can be dependent on the
degree of trust one can put on the available experimental results relevant to
the analysis.

The posterior probability distribution illustrates the degree of credibility
attributed to a model given the quality of its reproduction of experimental
data. This distribution can be marginalized over a subset of observables to
obtain marginal probability distribution functions (PDFs) on a single quantity
$o_{j}$

$\mathcal{P}(o_{j}|{\rm data})=\bigg{(}\prod_{i\neq k}\int do_{i}\bigg{)}\
\mathcal{P}(\\{o_{i}\\}|{\rm data})\ .$ (14)

###  III.2 Numerical setup

####  III.2.1 Prior knowledge

The first step consists in generating an uninformative prior on the
parameters. We span the model uncertainty by drawing each of the seven model
parameters random and independently, using a uniform prior distribution. The
interval boundaries considered for each of the model parameters are summarized
in Table 1. These ranges were chosen heuristically from previous SU$(3)_{f}$
parametrizations that have been suggested in the literature using similar
models [19, 29, 36]. Values far outside of these boundaries may typically lead
to unphysical behavior (unbound mesons or baryons, absence of resonances,
unusually large chiral condensate…). One can check a posteriori that they
indeed span a satisfactory uncertainty on the hadron spectrum masses. Note
that we choose as independent variables the quantities ($G\Lambda^{2}$,
$G_{d}/G$, $G_{d,V}/G$) instead of ($G$, $G_{d}$, $G_{d,V}$) for the
convenience of working with dimensionless quantities that are easier to
interpret.

Parameter | $\Lambda$ | $G\Lambda^{2}$ | $K\Lambda^{5}$ | $m_{0}$ | $m_{0,s}$ | $G_{d}/G$ | $G_{d,V}/G$
---|---|---|---|---|---|---|---
| (MeV) |  |  | (MeV) | (MeV) |  |
Min | 550 | 1.75 | 4 | 4 | 95 | 1.0 | 0.75
Max | 625 | 4 | 15 | 7 | 150 | 1.89 | 1.3
Table 1: Minimum and maximum values employed for the prior on each of the
independent parameters of the SU$(3)_{f}$ NJL model.

From the prior distribution, we first generate a large number ($N=5\times
10^{7}$) of random parameter sets. For each one of these parameter sets, the
self-consistent Dyson equations are solved in order to find the vacuum quark
masses, and the masses of the meson pseudoscalar octet as well as baryons from
the octet and decuplet are computed by finding the poles of the corresponding
propagator. All models that cannot account for the full set of the chosen
hadrons are considered as unrealistic; if any of the hadronic states is either
found not to exist or to be unbound (unstable) within the model, the model is
discarded from the analysis. This enforces that our prior model population is
always able to predict at least the existence of the correct spectrum for the
mesons and baryons of interest. This “stability condition” imposes

$0\leq m_{X}\leq m_{f_{1}}+m_{f_{2}}$ (15)

for the mesons and diquarks constituted of a quark and (anti)quark with
flavors $f_{1}$ and $f_{2}$ 111The $\eta^{\prime}$ meson is oftentimes found
to be unstable within this model and violates Eq. (15), as its mass remains
quite large due to the U(1)A anomaly [19]. We will therefore make an exception
for this meson and compute its mass as the real part of the complex pole of
the T-matrix. The associated width is an artifact due to the lack of
confinement in the NJL model.. For the baryons, the analogue condition reads

$0\leq m_{B}\leq m_{d}+m_{q}\ ,$ (16)

where $d$ and $q$ denote the diquark and quark that compose the baryon $B$.
Note that if a baryon state is constructed as a superposition of several
quark-diquark states, we enforce this inequality for all the quark-diquark
pairs involved.

####  III.2.2 Bayesian likelihood from experimental data

In order to select the models that are the most consistent with
phenomenological knowledge, likelihood weights are assigned depending on their
quality to reproduce the low-energy hadronic spectrum properties. The
quantities that we compute and confront with experimental data are

  1. (i)

the masses of the pseudoscalar mesons: $m_{\pi}$, $m_{K}$, $m_{\eta}$,
$m_{\eta^{\prime}}$,

  2. (ii)

the pion decay constant and the ratio between kaon/pion decay constants:
$f_{\pi}$, $f_{K}/f_{\pi}$ ,

  3. (iii)

the light quark chiral condensate:
$\Sigma=-\langle\overline{q}q\rangle^{\frac{1}{3}}$ ,

  4. (iv)

the masses of the baryon octet: $m_{N}$, $m_{\Lambda}$, $m_{\Sigma}$,
$m_{\Xi}$,

  5. (v)

the masses of the baryon decuplet: $m_{\Delta}$, $m_{\Sigma^{\star}}$,
$m_{\Xi^{\star}}$, $m_{\Omega}$.

Note that since we assume isospin symmetry, isospin degenerate states have
been gathered into the same state. In principle, scalar mesons could also be
used to further constrain the parameters, but there is some ambiguity in the
identification of the NJL scalar resonances to the physically observed states.
Their masses are in addition experimentally more difficult to measure [49] and
their rigorous modelling is outside the scope the NJL model [50].

We determine the likelihood by assigning to the models a Gaussian weight for
each quantity $o_{i}$ (except for the chiral condensate; see below),

$w_{i}=\frac{1}{\sqrt{2\pi}\sigma_{i}}\exp\Big{[}-\frac{(o_{i}-\overline{o_{i}})^{2}}{2\sigma_{i}^{2}}\Big{]}\
,$ (17)

where $\overline{o_{i}}$ is the recommended empirical value for the quantity
$o_{i}$, and $\sigma_{i}$ the associated uncertainty. The total likelihood
function is obtained by combining the weights associated to each quantity,

$\mathcal{P}({\rm data}|\\{o_{i}\\})=\prod_{i}w_{i}\ ,$ (18)

where the index $i$ runs over all physical quantities that need to be
constrained. All the mean values $\overline{o_{i}}$ and standard deviations
$\sigma_{i}$ used to implement the Gaussian filters are detailed in Table 2.
We discuss now how these quantities were determined.

For the pion and kaon masses, these values are not necessarily straightforward
to choose since the neutral pion (kaon) and the charged pion (kaon) do not
have the same mass. According to the Particle Data Group (PDG), the
experimental masses are [51]:

$m_{\pi^{0}}=134.9768(5)~{}\text{MeV},\quad
m_{\pi^{\pm}}=139.57039(18)~{}\text{MeV}\ ,$ (19)
$m_{K^{0}}=497.611(13)~{}\text{MeV},\quad
m_{K^{\pm}}=493.677(16)~{}\text{MeV}\ .$ (20)

Since the experimental uncertainties are very small, we can neglect them
compared to systematic uncertainties arising from the various approximations
made within the model. However, since we assume $m_{u,0}=m_{d,0}$ and hence
all our results are isospin independent, there is a problem of identification
of the NJL pion/kaon with the real life hadrons. In the real world, the masses
of $\pi^{0}$ and $\pi^{\pm}$ are the result of the combination of strong,
electromagnetic and weak effects, while in the NJL description only the strong
interaction is taken into account, leading to a degeneracy in our approach.
This problem is discussed in details in Ref. [52] where electromagnetic
corrections were applied to the NJL model to describe the mass splitting of
the pions. This is however beyond the scope of this work. Therefore, we
estimate the mass of the isospin symmetric state to a barycenter of the
physical states, and the systematic uncertainty is simply related to the mass
splitting of the states:

$\overline{m_{\pi}}=\frac{2m_{\pi^{\pm}}+m_{\pi^{0}}}{3}\,,\qquad\sigma_{m_{\pi}}=m_{\pi^{\pm}}-m_{\pi^{0}}$
(21)
$\overline{m_{K}}=\frac{m_{K^{\pm}}+m_{K^{0}}}{2}\,,\qquad\sigma_{m_{K}}=m_{K^{\pm}}-m_{K^{0}}$
(22)

For the $\eta$ and $\eta^{\prime}$ mesons, the PDG reports:

$m_{\eta}=547.862(17)~{}\text{MeV},\quad
m_{\eta^{\prime}}=957.78(6)~{}\text{MeV}$ (23)

While the state identification is unequivocal in this case, the experimental
uncertainties are very small compared to the expected systematics of our
modeling approach. Therefore, we choose to scale the mass uncertainty of these
mesons over their mass similar to the pion ratio $\sigma_{m_{\pi}}/m_{\pi}$:

$\sigma_{m_{\eta}}=\frac{m_{\eta}}{m_{\pi}}\sigma_{m_{\pi}}\,,\qquad\sigma_{m_{\eta^{\prime}}}=\frac{m_{\eta^{\prime}}}{m_{\pi}}\sigma_{m_{\pi}}$
(24)

Concerning the pion decay constant, similar considerations on the isospin
dependence are in order. It has been argued that isospin symmetry breaking
brings only a second order correction in chiral perturbation theory, such that
$f_{\pi^{\pm}}\approx f_{\pi^{0}}$ should be a good approximation [53].
Unfortunately, leptonic decay constants cannot be measured directly from
experiments as only their product with CKM matrix elements appear in their
(measurable) weak decay width. Therefore, their determination must rely on
some additional theoretical framework (typically LQCD together with chiral
perturbation theory). From Ref. [53] the theoretical preferred value for
$f_{\pi}$ is:

$f_{\pi}=\frac{1}{\sqrt{2}}(130.2\pm 1.7)\approx 92.1\pm 1.2~{}\text{MeV}$
(25)

The factor $1/\sqrt{2}$ is due to a different choice of normalization to that
of the NJL literature. To take systematic uncertainties into account, the
recommended uncertainty was inflated by 10%.

Surprisingly, the quantity $f_{K}/f_{\pi}$ is not very well reproduced by the
NJL model [15]. Anticipating our results, the predicted values of the ratio
lie typically between 1 and 1.1, while the experimental value of Ref. [53] is
${f_{K}/f_{\pi}=1.1928(26)}$. Such a discrepancy indicates a large systematic
error coming from the model for this quantity that is not well understood. To
still favor larger values of this ratio, we increase the systematic
uncertainty to $\sigma_{f_{K}/f_{\pi}}=0.075$.

The chiral condensate $\Sigma=-\langle\overline{q}q\rangle^{\frac{1}{3}}$ is
another vacuum parameter that cannot be accessed easily by experiments. We
must therefore rely again on ab initio approaches as LQCD or QCD sum rules
approaches (QCDSR) to obtain estimates. Early estimations from LQCD with two
flavors [54] led to the result:

$\Sigma=231\pm 11~{}\text{MeV}$ (26)

at a renormalization scale of 1 GeV. Estimations for the QCDSR, on the other
hand, suggest the wider range [20, 55]:

$\Sigma=229^{+33}_{-36}~{}\text{MeV}$ (27)

Today, the ranges recommended by the Flavour Lattice Data Group (FLAG) for two
and three flavors are [56]:

$\Sigma=266\pm 10~{}\text{MeV}~{}~{}(\text{2f})\,,\quad\Sigma=272\pm
5~{}\text{MeV}~{}~{}(\text{3f}),$ (28)

which are larger than the previous estimations. These results were however
obtained in the chiral limit (vanishing $m_{u}$ and $m_{d}$) and at the
renormalization scale of 2 GeV (which is conventional in LQCD simulations).

With all these conflicting results, it is difficult to find a proper
restriction on the chiral condensate $\Sigma$. It seems to be reasonable to
exclude models predicting chiral condensates larger than ${\Sigma_{{\rm
max}}\approx 275~{}}$MeV and smaller than $\Sigma_{{\rm min}}\approx$ 210 MeV.
Therefore, in the conservative approach we implemented a Gaussian passband
filter to constrain the chiral condensate in this range:

$w_{\Sigma}=\frac{1}{4}\Big{[}1-\operatorname{erf}\Big{(}\frac{\Sigma-\Sigma_{{\rm
max}}}{\sqrt{2}\sigma_{\Sigma}}\Big{)}\Big{]}\Big{[}1+\operatorname{erf}\Big{(}\frac{\Sigma-\Sigma_{{\rm
min}}}{\sqrt{2}\sigma_{\Sigma}}\Big{)}\Big{]}\ ,$ (29)

with an estimate of the standard deviation $\sigma_{\Sigma}=10$ MeV.

For all baryons of the octet, we have naturally very precise measurements of
their experimental masses from the PDG. Just like for mesons, the isospin
partners are degenerate in our approximation and we determine the mass of the
isospin-averaged state (that we will call nucleon, $\Lambda$, $\Sigma$ and
$\Xi$ baryons) by calculating the barycenter of the masses of the different
physical states. Clearly, the uncertainty coming from the isospin mass
splittings ($\approx 1~{}$MeV for the nucleon, $\lesssim 10~{}$MeV for the
$\Sigma$ and $\Xi$ baryons) are very small. The difficulty lies in the
estimation of the systematic uncertainty of our formalism for the construction
of baryons, that relies on several nested approximations (NJL effective model,
mean field approach, quark-diquark reduction of the 3-body equation and static
approximation). To be on the conservative side we choose to fix the nucleon
mass uncertainty with ${\sigma_{m_{N}}=50\times(m_{n}-m_{p})}\approx
60~{}$MeV, and scale this uncertainty to the other octet baryons with the same
method, as we did for $\eta$ mesons in Eq. (24).

In the case of the $\Delta$ baryon, the PDG database reports a Breit-Wigner
mass $m_{\Delta}=1232\pm 2$ MeV and a full width $\Gamma_{\Delta}=117\pm 3$
MeV. The masses of the individual $\Delta^{++}$, $\Delta^{+}$, $\Delta^{0}$
and $\Delta^{-}$ resonances are not known very precisely. Since the decay of
the $\Delta$ baryon is not predicted within our model, we suggest that the
width of the resonance is a good indicator of the scale of precision that
cannot be accounted for by the NJL model, such that we choose
$\sigma_{m_{\Delta}}=\Gamma_{\Delta}$. For the remaining baryons of the
decuplet ($\Sigma^{*}$, $\Xi^{*}$ and $\Omega$) we determine again the isospin
symmetric mass by averaging the physical masses from the PDG database, and
scale $\sigma_{m_{\Delta}}$ to the mass of the corresponding baryon in a
similar fashion to Eq. (24).

| $m_{\pi}$ | $f_{\pi}$ | $m_{K}$ | $f_{K}/f_{\pi}$ | $m_{\eta}$ | $m_{\eta^{\prime}}$ | $-\langle\overline{q}q\rangle^{1/3}$ |
---|---|---|---|---|---|---|---|---
$\overline{o}$ | 138.04 | 92.07 | 495.64 | 1.1928 | 547.862 | 957.78 | 210 - 275 |
$\sigma$ | 4.6 | 1.3 | 3.93 | 0.075 | 18.23 | 31.87 | 10 |
| $m_{N}$ | $m_{\Lambda}$ | $m_{\Sigma}$ | $m_{\Xi}$ | $m_{\Delta}$ | $m_{\Sigma^{\star}}$ | $m_{\Xi^{\star}}$ | $m_{\Omega}$
$\overline{o}$ | 938.92 | 1115.683 | 1193.15 | 1318.28 | 1232 | 1384.57 | 1533.40 | 1672.45
$\sigma$ | 64.66 | 76.84 | 82.18 | 90.79 | 117 | 131.49 | 145.62 | 158.83
Table 2: Mean and standard deviation used for the Gaussian weights associated
to each quantity. All quantities (except $f_{K}/f_{\pi}$ which is
dimensionless) are given in units of MeV.

##  IV Results and discussion

Figure 1: Posterior distributions of the SU$(3)_{f}$ NJL model parameters with
weights W1 (blue), W2 (red) and W3 (green); see text for details. All
histograms are normalized to unit area.

To compare the effects brought by each weight, they were gathered in three
weighting procedures that are then compared on the following figures:

  * •

W1 (represented in blue): No weight is assigned to the data, i.e.
$\mathcal{P}({\rm data}|\\{o_{i}\\})=1$ according to Eq. (18). The
distribution simply reflects the natural predictions of the model in the
explored prior parameter space. However we will see that the additional
conditions on the stability of the bound states, which may seem minor at
first, already affect the model selection process.

  * •

W2 (represented in red): All weights associated to the properties of the meson
octet (including $f_{\pi}$ and $f_{K}/f_{\pi}$), as well as the chiral
condensate filter are implemented.

  * •

W3 (represented in green): The weights associated to each baryon of the octet
and decuplet are added on top of the ones already in W2.

In this section, we will compare the posterior probability distribution
functions (PDFs) conditioned by each weighting procedure W1, W2 and W3, in
order to reveal the role that is played by the different constraints. As a
summary of our results, we provide in App. E a table collecting the predicted
mean values and uncertainties associated to each quantity examined in this
paper for the weighting W3, that is when all possible physical constraints
have been applied.

###  IV.1 Model parameters

Figure 2: Posterior distributions of the quark mass and chiral condensates for
light and strange flavors with weights W1 (blue), W2 (red) and W3 (green) (see
main text for details). All histograms are normalized to unit area.

In Fig.1, we show the PDFs of the five model parameters (leaving the
discussion of the diquark parameters to section IV.2). As expected, if no
weights are assigned to the data (W1), the probability distributions remain
almost flat, reflecting our use of a flat prior for these parameters. The
deviations from a flat distribution (particularly notable for the quantity
$G\Lambda^{2}$) epitomize the effect of selecting models that are able at
least to produce a stable spectrum of mesons and baryons. Low values of
$G\Lambda^{2}$ ($\lesssim 2$) typically lead to relatively small quark masses,
which may damage the stability of baryons due to Eq. (16). Once we enforce
physical values for the hadron spectrum (W2 and W3), the distributions of all
parameters become peaked around a preferred value, which validates our choice
for the prior. The parameter set employed in a previous work using the same
set of approximations is in good agreement with these results [29]. We should
mention that we also checked with larger priors that there is no statistically
significant increase of these PDFs outside of the parameter ranges used in
this work.

The PDFs of the dressed quark masses and chiral condensates are presented in
Fig.2. The light quark mass is found to lie between 315 MeV and 1015 MeV
within our parameter prior. Once the W2 filter on meson masses is implemented,
this range is narrowed down to $m_{u,d}\lesssim 600~{}$MeV, with a preference
around 430 MeV. The reproduction of baryon masses disfavors smaller values and
narrows down further the light quark mass to $m_{u,d}=515~{}\pm 53$ MeV,
although W2 and W3 remain largely compatible. A similar pattern is observed in
the strange sector. Such an effect can be explained by the requirement of a
finite binding energy for both diquarks and baryons as we shall discuss below.
In addition, it is interesting to notice that values close to the naive
estimate ${m_{u,d}\approx m_{N}/3\approx 313~{}}$MeV are rejected already by
W1, implying that low quark masses are not compatible with the reproduction of
the full hadron spectrum. This is in agreement with the findings of Ref. [22].

The light chiral condensate
$\Sigma=-\langle\overline{q}q\rangle^{\frac{1}{3}}$ is found to lie in the
range $218~{}$MeV$\,<\Sigma<\,$282 MeV in our prior. This is already in good
agreement with previous empirical knowledge on this quantity, implying that
the effects of the weight $w_{\Sigma}$ (29) will be overall very small.
Imposing the matching of the meson spectrum only gives a small range of
compatible values ${\Sigma=243\pm 6~{}}$MeV, which is not affected much by the
additional conditions of W3 on baryons, contrary to the quark mass
distribution. This result is quite small compared to the estimations from LQCD
(see Eq. 28).

Figure 3: Posterior distributions of the meson octet masses and decay
constants with weights W1 (blue), W2 (red) and W3 (green) (see main text for
details). The dashed lines indicate the likelihood weights implemented in W2
and W3 for each quantity. All histograms are normalized to unit area.

We show in Fig. 3 the PDFs of the observables for the pseudo-scalar octet,
i.e. meson masses and decay constants. For each quantity, we also show with a
dotted line the weight distribution of the likelihood illustrating the
constraints implemented in both W2 and W3. In the case of the pion mass, pion
decay constant and kaon mass, the posterior distributions of W2 and W3 follow
very well the likelihood, such that the uncertainties on these quantities are
accurately reproduced. On the other hand, the likelihood on the ratio
$f_{K}/f_{\pi}$ is poorly reproduced, as the model is known to be unable to
produce values larger than $\approx 1.05$ for this ratio [15], unless a very
small quark mass is adopted (see e.g. [57]). For the $\eta$ and
$\eta^{\prime}$ mesons, we see that while the physical values for the masses
of both states are able to be reached by the model, they unfortunately cannot
be simultaneously accounted for. Therefore, when trying to enforce physical
values, one needs to make a compromise with slightly larger
$m_{\eta}^{\prime}$ and smaller $m_{\eta}$. Both of these failures could be
indicators of the limits of the NJL model as well as the mean field and ring
approximation used to treat the problem. It also shows the limits of a
strategy of direct fitting of the model parameters to the data. For instance,
the ’t Hooft coupling $K$ is often used in order to fine tune the
$\eta-\eta^{\prime}$ mass splitting, but it turns out this cannot be done
without damaging the quality of the fit to other quantities. Treating all
experimental quantities on equal footing (factoring the uncertainties), we
guarantee to find the best middle ground the NJL model can offer.

###  IV.2 Predictions for the baryon spectrum

Figure 4: Posterior distributions of the scalar diquark coupling (a) and
axial-vector diquark coupling (b) with weights W1 (blue), W2 (red) and W3
(green) (see main text for details). All histograms are normalized to unit
area.

In Fig. 4, we show the prior and posterior PDFs of the diquark coupling ratios
$G_{d}/G$ and $G_{d,V}/G$. In both cases, we see that relatively large
coupling ratios are needed to sustain the existence of the baryon states. The
meson and baryon constraints altogether (W3) select rather precise values for
the ratios, $G_{d}/G=1.53\pm 0.04$ and $G_{d,V}/G=0.97\pm 0.06$. These results
are surprisingly far from those obtained via the Fierz formulas, which suggest
$G_{d}/G=3/4$ and $G_{d,V}=3/8$, respectively. This confirms that one should
take caution when applying these relationships within a NJL-like effective
model. This value of the diquark coupling $G_{d}/G$ is also typically larger
than the one used in the literature to study the transition to color
superconducting matter at high density [58, 59, 60, 61].

Even though the quantity $G_{d}/G$ affects only the diquark and baryon
sectors, we can also observe in Fig.4a a strong preference towards large
values with weighting W2 (which remains agnostic on the baryon masses). This
is explained by the preference of W2 towards smaller quark masses, as baryons
require stronger binding to remain stable in this case.

Figure 5: Posterior distributions of the scalar diquark binding energy (left)
and axial-vector diquark binding energy (right) with weights W1 (blue), W2
(red) and W3 (green) (see main text for details). All histograms are
normalized to unit area.

In Fig. 5 we show the posterior distributions of the binding energies
$B_{d}=2m_{q}-m_{d}$ for both of the non strange scalar and axial diquarks. We
can observe a striking difference of behavior between the two species; while
the axial diquark is allowed to be very weakly bound within the NJL model,
this is not the case for the scalar diquark, whose binding energy is bounded
to be larger than $B_{(ud)}^{min}\approx$ 500 MeV. This behavior can be
explained by the fact that even with W1, our description of the nucleon
requires a very large diquark coupling to be bound, which prevents the diquark
to be weakly bound. This effect is not present for the axial diquark because
of its larger mass and the lower value of the diquark-vector coupling required
to fit the physical $\Delta$ mass. We observed a similar behavior for the
other scalar and axial diquarks in the strange sector. It should be stressed
that this feature is not a consequence of our optimization of the parameters
to reproduce experimental results, but is inherent to the NJL model (at least
within our approximations).

This property explains why our assumption of the existence of stable baryons
(and the nucleon in particular) requires a sufficiently large quark mass.
Indeed, if we rewrite the stability conditions (Eq.(15) and (16)) in the form:

$3m_{u/d}=m_{N}+B_{(ud)}+B_{N},$ (30)

where ${B_{N}=m_{(ud)}+m_{u/d}-m_{N}\geq 0}$ is the nucleon binding energy, it
becomes apparent that the minimum possible value for the quark mass is
${m_{u/d}=(m_{N}+B_{(ud)}^{min})/3\approx 480~{}}$MeV when $m_{N}$ is fixed to
its physical value.

Figure 6: Posterior distributions of the baryon binding energy
${B_{B}=m_{d}+m_{q}-m_{B}}$ and energy per quark difference ${\Delta
E_{B}=m_{B}/3-m_{d}/2}$ for the nucleon (left column) and the $\Delta$ baryon
(right column). All histograms are normalized to unit area.

Finally, we compare in the top panel of Fig. 6 the posterior distributions of
the binding energy of the nucleon and $\Delta$ baryon. We notably see that the
constraints from W3 restrict the nucleon binding energy to a maximum of about
100 MeV, while values as large as 500 MeV can be obtained for the $\Delta$
baryon. The most probable scenario is that of a very weakly bound nucleon with
$B_{N}<30~{}$MeV. Since our framework describes both the nucleon and its
constituent diquark with a single coupling constant $G_{d}$, we can expect
that the diquark and nucleon binding are very heavily correlated. This lack of
freedom can explain why the only possibility compatible with the physical
value of $m_{N}$ is that of a very unbalanced repartition of the binding
energies with a strongly bound diquark and a weakly bound nucleon. Because of
the larger mass of the $\Delta$, both diquark and baryon can bear reasonable
binding energies. This discrepancy is also related to the larger value of the
scalar coupling compared to the vector one.

On the bottom panel, we also compare the distributions of the quantity $\Delta
E_{B}=m_{B}/3-m_{d}/2$ which represents the energy per quark difference
between the baryon and the diquark. While this difference is always positive
for the nucleon in our model, it is almost always negative for the $\Delta$
baryon. The same pattern is observed for every other members of the baryon
octet and decuplet. This illustrates again the difference of behavior between
the two baryon families in our framework; in the case of the octet, the
diquark is very strongly bound and therefore very light, whereas the binding
energy is well distributed between the axial diquarks and baryons of the
decuplet. The sign of $\Delta E_{B}$ may have strong implications in medium
due to the lack of color confinement of the model [62, 63]. The fact that it
is positive for the baryon octet means that scalar diquarks are energetically
favored over baryons. Since there are no confinement mechanism to prevent the
diquark colored state to remain, they are allowed to survive in medium. At
finite chemical potential and low temperatures, this may lead to the formation
of a diquark Bose-Einstein condensate [64, 65].

###  IV.3 Correlations

To complete our analysis of the posterior distributions, we provide
correlation tables to understand the interplay on one side between the model
parameters and the pseudoscalar sector in Fig. 7, and on the other side
between the model parameters and the baryon sector in Fig. 8. The correlation
between two quantities $o_{i}$ and $o_{j}$ is interpreted through their linear
Pearson coefficient $C_{ij}$ defined by

$C_{ij}=\frac{\operatorname{\mathbb{E}}[o_{i}o_{j}]-\operatorname{\mathbb{E}}[o_{i}]\operatorname{\mathbb{E}}[o_{j}]}{\sqrt{\operatorname{\mathbb{E}}[o_{i}^{2}]-\displaystyle\operatorname{\mathbb{E}}[o_{i}]^{2}}\sqrt{\operatorname{\mathbb{E}}[o_{j}^{2}]-\operatorname{\mathbb{E}}[o_{j}]^{2}}}\,,$
(31)

where $\operatorname{\mathbb{E}}$ is the expectation operator. On each figure,
the coefficients are shown for the probability distributions obtained with
weightings W1 and W3 in the lower left and upper right triangle, respectively.
This allows us to distinguish between the correlations coming from the
modeling itself, and the ones created by the constraints we impose on each
parameter sets. In the following, we highlight the most important conclusions
that we draw from these two figures.

First, we notice that both the light and strange quark masses are very
strongly correlated to the dimensionless scalar coupling $G\Lambda^{2}$. This
is not surprising, as the scalar quartic term is responsible for the mechanism
of spontaneous symmetry breaking in the theory. However, this shows that the
second term of Eq. (2) remains dominant against the bare mass and the ’t Hooft
term for all possible parametrizations. In addition, we conclude that the
quantity $G\Lambda^{2}$ is the most relevant to control the value of both
quark masses, which are in turn heavily correlated with each other. On the
other hand, the chiral condensates (especially the strange one) are better
characterized by the value of the momentum cutoff $\Lambda$. In the absence of
constraints, the light bare mass only affects the pion mass, which can be
easily understood from its origin as a pseudo-Goldstone boson. Once the pion
mass is fixed to its physical value, the value of $m_{0,u}$ is adapted to
satisfy the condition to bring $m_{\pi}$ to its physical value, which brings
minor correlations with the masses and condensates. In a similar fashion, the
quantity $K\Lambda^{5}$ influences the most the masses of the $\eta$ and
$\eta^{\prime}$ mesons due to its role in breaking the U(1)A symmetry, but
inherits significant anticorrelations with the quark masses once the meson
masses are filtered.

Regarding the baryon sector (Fig. 8), we observe, as expected, strong
correlations between the diquark and baryon sectors. In particular, in the
case of W1, the masses of the baryons are strongly determined by the mass of
the corresponding diquark, and to a lesser extent by the quark mass. In
addition, strong a anticorrelation is observed between the diquark couplings
which are naturally proportional to the baryon binding energy. After matching
with the physical meson and baryon masses, the role of the quark mass $m_{u}$
becomes much more significant, the binding energies $B_{(ud)}$, $B_{[ud]}$,
$B_{N}$ and $B_{\Delta}$ having Pearson coefficients close to unity with
$m_{u}$. These correlations are easy to understand from Eq. (30) (and its
equivalent for the $\Delta$ baryon), where the relationship between the total
binding energy and the quark mass appears to be linear assuming the bound
state mass to be fixed. This emphasizes the central role of the quark mass as
an energy scale in the problem, as its value remains not very well constrained
by the physical input (see Fig.2). This correlation with the quark mass is
also very strong for the mass of the axial diquark $m_{[ud]}$, but
surprisingly disappears completely for the scalar diquark mass $m_{(ud)}$ once
the nucleon mass is fixed. We expect that this effect is related to the
different behavior of the scalar diquark observed in section IV.2 as well as
the larger value of the diquark coupling required in scalar sector to match
the mass of the baryon octet, which prevents large deviations of $m_{(ud)}$
from its average value. This difference is also illustrated by the larger
sensibility around the physical point of the nucleon mass towards the mass of
its diquark component ($C_{ij}\approx 1$) than for the $\Delta$ baryon
($C_{ij}\approx 0.4$).

Note that just like in Sec. IV.2, we only discuss the correlation coefficients
related to the nucleon and $\Delta$ baryon (together with their corresponding
diquark), as we observed a very similar correlation pattern for the other
baryons belonging to the same multiplet.

Figure 7: Correlation coefficients between the parameters and observables in
the meson sector. The bottom-left half contains correlations when no weighting
procedure is applied to the data (weight W1; see Sec.IV), while the top right
corner highlights the remaining correlations after all the constraints have
been applied (weight W3). Figure 8: Correlation coefficients between the
parameters and observables in the baryon sector. The bottom-left half contains
correlations when no weighting procedure is applied to the data (weight W1),
while the top right corner highlights the remaining correlations after all the
constraints have been applied (weight W3).

##  V Summary

In the present work, we have performed a Bayesian analysis of the properties
of the meson and baryon spectra as predicted by the SU$(3)_{f}$ Nambu–Jona-
Lasinio model. The Bethe-Salpeter equations in the RPA approximation as well
as the Faddeev equations in a quark-diquark static approximation were solved
to compute the masses of the lowest mesonic and baryonic states for a large
variation within the parameter space. The results were then confronted with
the experimentally measured masses of the physical states to select plausible
models.

As it has been shown by previous studies, we confirm that the NJL model is
able to satisfactorily describe the main traits of QCD phenomenology in the
vacuum, including the spectrum of the pseudoscalar meson octet and of the
baryon octet and decuplet. Nevertheless, some discrepancies were observed for
the $\eta$ and $\eta^{\prime}$ masses as well as the kaon/pion decay constant
ratio $f_{K}/f_{\pi}$. This global view is new in our approach, as in previous
work the parameter set had been chosen to describe a limited and arbitrary set
of experimental observables.

These results highlights the importance of taking into account the maximum
number of constraints to provide a good parametrization of the model.
Neglecting some of the physical input or the uncertainties associated to the
measurement and modeling of the quantities described may lead to
misconceptions regarding the predictive power of the model and overlook the
subsequent uncertainties on the model parameters and predictions. In addition,
our Bayesian framework allows us to inspect the correlations between
parameters and predictions to better understand the interplay between each
quantity.

In order for the predicted baryon masses to be in accordance with the physical
masses, we find that large values of the scalar and vector diquark couplings
must be required (${G_{d}/G=1.53\pm 0.04}$ and ${G_{d,V}/G=0.97\pm 0.06}$). In
the scalar case, we find that the diquarks are systematically strongly bound
(${B_{(ud)}\gtrsim 500}$ MeV), which in turn also imposes the quark mass to be
sufficiently large (${m_{u/d}\gtrsim 480~{}}$MeV). In consequence, the average
binding energy of a quark in the nucleon is always smaller that its average
binding energy in the diquark which is a constituent. These findings could
possibly have significant impact on the phenomenology of dense quark matter in
compact stars, where the strength of the diquark coupling could lead to a
strong effect of quark pairing and color superconductivity. This relies
however on the assumption that the large diquark couplings could be maintained
in-medium, as the density dependences of the couplings are very much
uncertain.

###### Acknowledgements.

This study is part of a project that has received funding from the European
Union’s Horizon 2020 research and innovation program under grant agreement
STRONG – 2020 - No 824093. The work of J.M.T.-R. is supported by Project no.
PID2020-118758GB-I00, financed by the Spanish MCIN/
AEI/10.13039/501100011033/, and Project no. 315477589 - TRR 211 (Strong-
interaction matter under extreme conditions), financed by the Deutsche
Forschungsgemeinschaft.

## Appendix

###  A One-quark loop integral

The loop integral $I_{1}$ appears in the calculation of the chiral condensates
Eq. (5) and it is given by,

$I_{1}(m,\Lambda)=\int^{\Lambda}\frac{d^{4}k}{(2\pi)^{4}}\frac{1}{k^{2}-m^{2}+i\varepsilon}\
.$ (A.1)

In the 3D noncovariant momentum cutoff scheme, one obtains after applying
Cauchy’s integral formula,

$\begin{split}iI_{1}(m,\Lambda)&=\int^{\Lambda}\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2E_{k}}\\\
&=\frac{\Lambda^{2}}{8\pi^{2}}\Big{[}\sqrt{1+\frac{m^{2}}{\Lambda^{2}}}-\frac{m^{2}}{\Lambda^{2}}\log\Big{(}\frac{\Lambda+\sqrt{\Lambda^{2}+m^{2}}}{m}\Big{)}\Big{]}\
,\end{split}$ (A.2)

where $E_{k}=\sqrt{{\bf k}^{2}+m^{2}}$. This expression can be used to solve
the self-consistent gap equation (2) in order to find the vacuum quark masses.

###  B Meson and diquark polarization functions

We start by recalling the general expression of the loop function associated
to the bound state $X$,

$\Pi_{X}(p)=i\int^{\Lambda}\frac{d^{4}k}{(2\pi)^{4}}\operatorname{Tr}\Big{(}\Gamma_{X}S(p+k)\Gamma_{X}^{\dagger}S(k)\Big{)}\
,$ (B.1)

where the dressed fermion propagator $S(k)$ is given by the Eq. (4). All the
interaction vertices $\Gamma_{X}$ in spinor, flavor and color space associated
to each bound state are gathered in Table 3. Pseudo-scalar mesons have
negative parity, yielding a factor $i\gamma_{5}$ in the spinor space. They are
all colorless states, so we select the singlet color representation. For the
pion and kaon, since we work under the isospin approximation one can choose in
flavor space any representant of the pion triplet ($\pi^{+}$, $\pi^{-}$,
$\pi^{0}$) and kaon quadruplet ($K^{+}$, $K^{-}$, $K^{0}$,
$\overline{K^{0}}$); here we choose the $\pi^{0}$ and $K^{0}$ channels for
illustrative purposes. For the $\eta$ and $\eta^{\prime}$ mesons, a
diagonalization in the flavor subspace (0,8) is necessary since the
polarization matrix has non-zero off-diagonal terms. We refer to Ref. [19] for
the technical details on this procedure (see also Ref. [27] for an extension
without the isospin approximation).

In the diquark case, scalar diquarks are needed for the construction of the
baryon octet, while axial diquarks are required to describe the baryon
decuplet. Their spinor channel includes an additional $C$ matrix in order to
select a particle-particle scattering channel. Scalar diquarks belong to the
antisymmetric flavor antitriplet
$\overline{\bm{3}}_{f}\subset\bm{3}_{f}\otimes\bm{3}_{f}$. In the isospin
approximation, we can choose $(us)$ as a representative of the degenerate pair
together with $(ds)$. Axial diquarks belong to the symmetric flavor sextet
${\bm{6}}_{f}\subset{\bm{3}_{f}\otimes\bm{3}_{f}}$. We choose $[ud]$ as the
representative of the isospin triplet ($[ud]$, $[uu]$, $[dd]$) and $[us]$ as
the representative of the degenerate pair together with $[ds]$. In order to
produce colorless baryons, diquarks must belong to the antisymmetric color
antitriplet representation of ${\bm{3}_{c}\otimes\bm{3}_{c}}$. They are
therefore colored objects with antisymmetric color wave functions $|rb-
br\rangle$, $|gb-bg\rangle$ or $|rg-gr\rangle$. The color triplets are of
course degenerate, and are combined together with quarks of the different
colors to form colorless baryons.

$X$ | Spinor | Flavor | Color
---|---|---|---
$\pi$ | $i\gamma^{5}$ | $\tau_{3}$ | $\mathbbb{1}_{c}$
$K$ | $i\gamma^{5}$ | $\frac{1}{\sqrt{2}}(\tau_{6}-i\tau_{7})$ | $\mathbbb{1}_{c}$
$\eta$ | $i\gamma^{5}$ | $(\tau_{0},\tau_{8})$ | $\mathbbb{1}_{c}$
$\eta^{\prime}$ | $i\gamma^{5}$ | $(\tau_{0},\tau_{8})$ | $\mathbbb{1}_{c}$
$(ud)$ | $Ci\gamma^{5}$ | $-\tau_{2}$ |  $-\lambda_{2}$/$\lambda_{5}$/$-\lambda_{7}$
$(us)$ | $Ci\gamma^{5}$ | $\tau_{5}$ |  $-\lambda_{2}$/$\lambda_{5}$/$-\lambda_{7}$
$[ud]$ | $C\gamma^{\mu}$ | $\tau_{1}$ |  $-\lambda_{2}$/$\lambda_{5}$/$-\lambda_{7}$
$[us]$ | $C\gamma^{\mu}$ | $\tau_{4}$ |  $-\lambda_{2}$/$\lambda_{5}$/$-\lambda_{7}$
$[ss]$ | $C\gamma^{\mu}$ | $\frac{1}{\sqrt{3}}\tau_{0}-\sqrt{\frac{2}{3}}\tau_{8}$ |  $-\lambda_{2}$/$\lambda_{5}$/$-\lambda_{7}$
Table 3: Vertex factors in spinor, flavor and color spaces for the mesons and
diquarks considered in this work.

After performing the traces, we obtain the respective polarization functions
for pseudoscalar mesons, scalar diquarks and axial diquarks,

$\Pi_{PS}^{f_{1}f_{2}}(p)=4iN_{c}\Big{[}I_{1}(m_{f_{1}},\Lambda)+I_{1}(m_{f_{2}},\Lambda)\\\
+((m_{f_{1}}-m_{f_{2}})^{2}-p^{2})I_{2}(p,m_{f_{1}},m_{f_{2}},\Lambda)\Big{]}\
,$ (B.2)

$\Pi_{d,S}^{f_{1}f_{2}}(p)=8i\Big{[}I_{1}(m_{f_{1}},\Lambda)+I_{1}(m_{f_{2}},\Lambda)\\\
+((m_{f_{1}}+m_{f_{2}})^{2}-p^{2})I_{2}(p,m_{f_{1}},m_{f_{2}},\Lambda)\Big{]}\
,$ (B.3)

$\Pi_{d,A}^{f_{1}f_{2}}(p)=-\frac{16}{3}i\Big{[}I_{1}(m_{f_{1}},\Lambda)+I_{1}(m_{f_{2}},\Lambda)\\\
+((m_{f_{1}}+m_{f_{2}})^{2}-p^{2})I_{2}(p,m_{f_{1}},m_{f_{2}},\Lambda)\Big{]}\
,$ (B.4)

where $f_{1}$ and $f_{2}$ are the two flavor of the quarks or antiquarks
involved in the bound state. Note that the diquark polarization functions do
not get an additional factor $N_{c}$ since they do not belong to a singlet
representation of ${\bm{3}_{c}\otimes\bm{3}_{c}}$. The two-line loop integral
$I_{2}$ is expressed by,

$I_{2}(p,m_{f_{1}},m_{f_{2}},\Lambda)=\int^{\Lambda}\frac{d^{4}k}{(2\pi)^{4}}\frac{1}{[(p+k)^{2}-m_{f_{1}}^{2}](k^{2}-m_{f_{2}}^{2})}\
.$ (B.5)

In the noncovariant cutoff scheme, and taking the bound state at rest
($\bm{p}=0$), one can reduce this expression to:

$I_{2}(p_{0},m_{f_{1}},m_{f_{2}},\Lambda)=\\\
-\frac{i}{8\pi^{2}p_{0}}\bigg{[}m_{f_{1}}\mathcal{I}\Big{(}\sqrt{1+\frac{\Lambda^{2}}{m_{f_{1}}^{2}}},-\frac{p_{0}^{2}+m_{f_{1}}^{2}-m_{f_{2}}^{2}}{2p_{0}m_{f_{1}}}\Big{)}\\\
-m_{f_{2}}\mathcal{I}\Big{(}\sqrt{1+\frac{\Lambda^{2}}{m_{f_{2}}^{2}}},\frac{p_{0}^{2}+m_{f_{2}}^{2}-m_{f_{1}}^{2}}{2p_{0}m_{f_{2}}}\Big{)}\bigg{]}\
,$ (B.6)

where the remaining integral $\mathcal{I}$ appearing in the previous
expression can be given analytically,

$\mathcal{I}(M,e_{0})=\int_{1}^{M}de\frac{\sqrt{e^{2}-1}}{e-e_{0}}=\Theta(M-1)\Bigg{(}\sqrt{M^{2}-1}+e_{0}\log(M+\sqrt{M^{2}-1})\\\
+\begin{cases}\sqrt{e_{0}^{2}-1}\log\Big{(}\Big{|}\frac{M-e_{0}}{-1+Me_{0}+\sqrt{(M^{2}-1)(e_{0}^{2}-1)}}\Big{|}\Big{)}\quad&\text{if}~{}|e_{0}|\geq
1\\\
-\sqrt{1-e_{0}^{2}}\arccos\Big{(}\frac{1-Me_{0}}{M-e_{0}}\Big{)}\quad&\text{if}~{}|e_{0}|<1\\\
\end{cases}+i\pi\sqrt{e_{0}^{2}-1}\ \Theta[(e_{0}-1)(M-e_{0})]\Bigg{)}\ ,$
(B.7)

with $\Theta$ being the Heaviside step function.

###  C Pion and kaon decay constants

The pion and kaon decay constants are computed applying Eq. (12) to the
corresponding channels (see the first two lines of Table 3). After performing
the traces, one obtains,

$f_{\pi/K}(p)=-2iN_{c}g_{\pi/K\rightarrow
ij}(m_{i}+m_{j})I_{2}(p,m_{i},m_{j},\Lambda)\ ,$ (C.1)

where the $i$ and $j$ subscripts denote any two light flavors in the pion
case, and one light flavor and one strange flavor for the kaon. The vacuum
values of the decay constants are then computed assuming the meson at rest and
on shell ($p_{0}=m_{X},{\bf p}=0$).

###  D Quark-diquark loop functions for baryons

The evaluation of the baryon propagator in the quark-diquark approximation is
made through the quark-diquark polarization function:

$\Pi_{B}^{f,d}(p)=i\int\frac{d^{4}k}{(2\pi)^{4}}S_{f}(p+k)\Delta_{d}(k)\ ,$
(D.1)

where $\Delta_{d}$ is the propagator of the diquark $d$, given by $T_{d}(p)$
in Eq. (10). Assuming the baryon at rest (${\bm{p}=0}$), we follow the
decomposition suggested in Appendix F of Ref. [29] and obtain,

$\Pi_{B}^{f,d}(p_{0})=-\frac{g_{d\rightarrow
qq}^{2}}{8\pi^{2}}\left[J_{F}^{+}(p_{0})+J_{F}^{-}(p_{0})+J_{B}^{+}(p_{0})+J_{B}^{-}(p_{0})\right]\
,$ (D.2)

where the functions $J_{F/B}^{\pm}$ can be expressed analytically,

$J_{F}^{\pm}(p_{0})=\pm\frac{1}{2p_{0}}m_{f}^{2}\Bigg{(}\mathcal{I}\bigg{(}\sqrt{1+\frac{\Lambda^{2}}{m_{f}^{2}}},\pm\frac{m_{d}^{2}-m_{f}^{2}-p_{0}^{2}}{2p_{0}m_{f}}\bigg{)}\mp\gamma_{0}\mathcal{I}^{\prime}\bigg{(}\sqrt{1+\frac{\Lambda^{2}}{m_{f}^{2}}},\pm\frac{m_{d}^{2}-m_{f}^{2}-p_{0}^{2}}{2p_{0}m_{f}}\bigg{)}\Bigg{)}\
,$ (D.3)
$J_{B}^{\pm}(p_{0})=\frac{1}{2p_{0}}m_{d}\Bigg{(}\pm(m_{f}+\gamma_{0}p_{0})\mathcal{I}\bigg{(}\sqrt{1+\frac{\Lambda^{2}}{m_{d}^{2}}},\pm\frac{m_{f}^{2}-m_{d}^{2}-p_{0}^{2}}{2p_{0}m_{d}}\bigg{)}+\gamma_{0}m_{d}\mathcal{I}^{\prime}\bigg{(}\sqrt{1+\frac{\Lambda^{2}}{m_{d}^{2}}},\pm\frac{m_{f}^{2}-m_{d}^{2}-p_{0}^{2}}{2p_{0}m_{d}}\bigg{)}\Bigg{)}\
,$ (D.4)

where the function $\mathcal{I}$ is shown in Eq. (B.7), while the function
$\mathcal{I}^{\prime}$ is given by,

$\mathcal{I}^{\prime}(M,e_{0})=\int_{1}^{M}de\,e\frac{\sqrt{e^{2}-1}}{e-e_{0}}=\frac{1}{2}\Theta(M-1)\Bigg{(}(2e_{0}^{2}+M)\sqrt{M^{2}-1}+(2e_{0}^{2}-1)\log(M+\sqrt{M^{2}-1})\\\
+\begin{cases}2e_{0}\sqrt{e_{0}^{2}-1}\log\Big{(}\Big{|}\frac{M-e_{0}}{-1+Me_{0}+\sqrt{(M^{2}-1)(e_{0}^{2}-1)}}\Big{|}\Big{)}\quad&\text{if}~{}|e_{0}|\geq
1\\\
-2e_{0}\sqrt{1-e_{0}^{2}}\arccos\Big{(}\frac{1-Me_{0}}{M-e_{0}}\Big{)}\quad&\text{if}~{}|e_{0}|<1\\\
\end{cases}+i\pi e_{0}\sqrt{e_{0}^{2}-1}\Theta((e_{0}-1)(M-e_{0}))\Bigg{)}\ .$
(D.5)

###  E Summary of the results

To summarize our results, we provide in Table. 4 the first two moments of the
posterior distributions with the weight W3 for the parameters and observables
investigated in this work. The former can be interpreted as a good estimate of
an optimal parametrization of the NJL model accounting for the constraints
imposed by W3. We indicate in addition the minimum and maximal values obtained
within our prior exploration for each quantity.

| Unit | Mean | $\sigma$ | Min | Max
---|---|---|---|---|---
$\Lambda$ | MeV | 568 | 8 | 550 | 625
$G\Lambda^{2}$ |  | 2.44 | 0.23 | 1.76 | 4.00
$K\Lambda^{5}$ |  | 10.05 | 1.35 | 4.00 | 15.0
$m_{0}$ | MeV | 5.69 | 0.42 | 4.00 | 7.00
$m_{0,s}$ | MeV | 134 | 5 | 95 | 150
$m_{u,d}$ | MeV | 515 | 53 | 315 | 1045
$m_{s}$ | MeV | 678 | 46 | 476 | 1160
$-\langle\overline{q}q\rangle^{\frac{1}{3}}$ | MeV | 243 | 4 | 219 | 282
$-\langle\overline{s}s\rangle^{\frac{1}{3}}$ | MeV | 250 | 4 | 235 | 284
$m_{\pi}$ | MeV | 138.2 | 4.5 | 110.3 | 187.4
$f_{\pi}$ | MeV | 92.1 | 1.2 | 82.6 | 101.8
$m_{K}$ | MeV | 497 | 4 | 437 | 554
$f_{K}/f_{\pi}$ |  | 1.01 | 0.01 | 0.98 | 1.15
$m_{\eta}$ | MeV | 510 | 8 | 357 | 604
$m_{\eta^{\prime}}$ | MeV | 987 | 27 | 680 | 1428
$G_{d}/G$ |  | 1.53 | 0.04 | 1.00 | 1.89
$G_{d,V}/G$ |  | 0.97 | 0.06 | 0.75 | 1.30
$m_{(ud)}$ | MeV | 489 | 29 | 0 | 1442
$m_{(us)}$ | MeV | 682 | 22 | 368 | 1556
$m_{[ud]}$ | MeV | 918 | 54 | 586 | 1721
$m_{[us]}$ | MeV | 1060 | 49 | 728 | 1832
$m_{[ss]}$ | MeV | 1197 | 46 | 852 | 1944
$m_{N}$ | MeV | 952 | 41 | 0 | 2454
$m_{\Lambda}$ | MeV | 1099 | 40 | 292 | 2560
$m_{\Sigma}$ | MeV | 1183 | 37 | 691 | 2577
$m_{\Xi}$ | MeV | 1294 | 35 | 768 | 2647
$m_{\Delta}$ | MeV | 1240 | 73 | 0 | 2376
$m_{\Sigma^{\star}}$ | MeV | 1389 | 70 | 181 | 2500
$m_{\Xi^{\star}}$ | MeV | 1538 | 67 | 361 | 2624
$m_{\Omega}$ | MeV | 1688 | 64 | 539 | 2752
$B_{(ud)}$ | MeV | 541 | 114 | 407 | 1984
$B_{[ud]}$ | MeV | 112 | 59 | 10 | 808
$B_{N}$ | MeV | 53 | 51 | 2 | 1100
$B_{\Delta}$ | MeV | 193 | 120 | 0 | 1951
Table 4: Summary of the posterior mean, standard deviation, and minimal and
maximal values for the model parameters and predicted hadronic observables
with the full weights W3.

## References

  * Durr  _et al._ [2008] S. Durr _et al._ , Ab-Initio Determination of Light Hadron Masses, Science 322, 1224 (2008), arXiv:0906.3599 [hep-lat] .
  * Fodor and Hoelbling [2012] Z. Fodor and C. Hoelbling, Light Hadron Masses from Lattice QCD, Rev. Mod. Phys. 84, 449 (2012), arXiv:1203.4789 [hep-lat] .
  * Dudek  _et al._ [2013] J. J. Dudek, R. G. Edwards, P. Guo, and C. E. Thomas (Hadron Spectrum), Toward the excited isoscalar meson spectrum from lattice QCD, Phys. Rev. D 88, 094505 (2013), arXiv:1309.2608 [hep-lat] .
  * Detmold  _et al._ [2019] W. Detmold, R. G. Edwards, J. J. Dudek, M. Engelhardt, H.-W. Lin, S. Meinel, K. Orginos, and P. Shanahan (USQCD), Hadrons and Nuclei, Eur. Phys. J. A 55, 193 (2019), arXiv:1904.09512 [hep-lat] .
  * Aarts  _et al._ [2017] G. Aarts, C. Allton, D. De Boni, S. Hands, B. Jäger, C. Praki, and J.-I. Skullerud, Light baryons below and above the deconfinement transition: medium effects and parity doubling, JHEP 06, 034, arXiv:1703.09246 [hep-lat] .
  * Bazavov  _et al._ [2019] A. Bazavov _et al._ , Meson screening masses in (2+1)-flavor QCD, Phys. Rev. D 100, 094510 (2019), arXiv:1908.09552 [hep-lat] .
  * García-Mascaraque  _et al._ [2022] S. C. García-Mascaraque, G. Aarts, C. Allton, T. Burns, S. Hands, and B. Jäger, Meson thermal masses at different temperatures, PoS LATTICE2021, 523 (2022), arXiv:2111.00784 [hep-lat] .
  * Blum  _et al._ [2016] T. Blum _et al._ (RBC, UKQCD), Domain wall QCD with physical quark masses, Phys. Rev. D 93, 074505 (2016), arXiv:1411.7017 [hep-lat] .
  * Fischer  _et al._ [2021] M. Fischer, B. Kostrzewa, M. Mai, M. Petschlies, F. Pittler, M. Ueding, C. Urbach, and M. Werner (Extended Twisted Mass, ETM), The $\rho$-resonance from $N_{f}$ = 2 lattice QCD including the physical pion mass, Phys. Lett. B 819, 136449 (2021), arXiv:2006.13805 [hep-lat] .
  * Paul  _et al._ [2022] S. Paul, A. D. Hanlon, B. Hörz, D. Mohler, C. Morningstar, and H. Wittig, I=1 $\pi$-$\pi$ scattering at the physical point, PoS LATTICE2021, 551 (2022), arXiv:2112.07385 [hep-lat] .
  * Roberts and Schmidt [2000] C. D. Roberts and S. M. Schmidt, Dyson-Schwinger equations: Density, temperature and continuum strong QCD, Prog. Part. Nucl. Phys. 45, S1 (2000), arXiv:nucl-th/0005064 .
  * Eichmann  _et al._ [2016] G. Eichmann, H. Sanchis-Alepuz, R. Williams, R. Alkofer, and C. S. Fischer, Baryons as relativistic three-quark bound states, Prog. Part. Nucl. Phys. 91, 1 (2016), arXiv:1606.09602 [hep-ph] .
  * Dorkin  _et al._ [2019] S. M. Dorkin, L. P. Kaptari, and B. Kämpfer, Pseudo-Scalar $\mathbf{q}\bar{\mathbf{q}}$ Bound States at Finite Temperatures, Few Body Syst. 60, 20 (2019), arXiv:1807.10075 [nucl-th] .
  * Gao and Ding [2020] F. Gao and M. Ding, Thermal properties of $\pi$ and $\rho$ meson, Eur. Phys. J. C 80, 1171 (2020), arXiv:2006.05909 [hep-ph] .
  * Klevansky [1992] S. P. Klevansky, The Nambu-Jona-Lasinio model of quantum chromodynamics, Rev. Mod. Phys. 64, 649 (1992).
  * Hatsuda and Kunihiro [1994] T. Hatsuda and T. Kunihiro, QCD phenomenology based on a chiral effective lagrangian, Physics Reports 247, 221–367 (1994).
  * Klimt  _et al._ [1990] S. Klimt, M. Lutz, U. Vogl, and W. Weise, Generalized SU(3) Nambu-Jona-Lasinio model: (I). Mesonic modes, Nuclear Physics A 516, 429 (1990).
  * Vogl and Weise [1991] U. Vogl and W. Weise, The nambu and jona-lasinio model: Its implications for hadrons and nuclei, Progress in Particle and Nuclear Physics 27, 195 (1991).
  * Rehberg  _et al._ [1996] P. Rehberg, S. P. Klevansky, and J. Hüfner, Hadronization in the SU(3) Nambu–Jona-Lasinio model, Physical Review C 53, 410–429 (1996).
  * Buballa [2005] M. Buballa, NJL-model analysis of dense quark matter, Physics Reports 407, 205–376 (2005).
  * Reinhardt [1990] H. Reinhardt, Hadronization of Quark Flavor Dynamics, Phys. Lett. B 244, 316 (1990).
  * Buck  _et al._ [1992] A. Buck, R. Alkofer, and H. Reinhardt, Baryons as bound states of diquarks and quarks in the Nambu-Jona-Lasinio model, Phys. Lett. B 286, 29 (1992).
  * Ebert  _et al._ [1994] D. Ebert, H. Reinhardt, and M. K. Volkov, Effective hadron theory of QCD, Prog. Part. Nucl. Phys. 33, 1 (1994).
  * Ishii  _et al._ [1995] N. Ishii, W. Bentz, and K. Yazaki, Baryons in the NJL model as solutions of the relativistic Faddeev equation, Nucl. Phys. A 587, 617 (1995).
  * Gastineau and Aichelin [2002] F. Gastineau and J. Aichelin, Strange baryons in a hot and dense medium within the Nambu-Jona-Lasinio model, J. Phys. G 28, 2017 (2002), arXiv:nucl-th/0201063 .
  * Wang  _et al._ [2011] J.-C. Wang, Q. Wang, and D. H. Rischke, Baryon formation and dissociation in dense hadronic and quark matter, Phys. Lett. B 704, 347 (2011), arXiv:1008.4029 [nucl-th] .
  * Blanquier [2011] E. Blanquier, Standard particles in the SU(3) Nambu-Jona-Lasinio model and the Polyakov model, Journal of Physics G: Nuclear and Particle Physics 38, 10.1088/0954-3899/38/10/105003 (2011).
  * Mu  _et al._ [2012] C.-f. Mu, Y. Jiang, P.-f. Zhuang, and Y.-x. Liu, Nucleons at finite temperature and low density in Faddeev equation approach, Phys. Rev. D 85, 014033 (2012).
  * Torres-Rincon  _et al._ [2015] J.-M. Torres-Rincon, B. Sintes, and J. Aichelin, Flavor dependence of baryon melting temperature in effective models of QCD, Physical Review C 91, 10.1103/PhysRevC.91.065206 (2015).
  * Fukushima [2004] K. Fukushima, Chiral effective model with the polyakov loop, Physics Letters B 591, 277–284 (2004).
  * Ratti  _et al._ [2006] C. Ratti, M. A. Thaler, and W. Weise, Phases of QCD: Lattice thermodynamics and a field theoretical model, Physical Review D 73, 10.1103/physrevd.73.014019 (2006).
  * Alvarez-Castillo  _et al._ [2016] D. Alvarez-Castillo, A. Ayriyan, S. Benic, D. Blaschke, H. Grigorian, and S. Typel, New class of hybrid EoS and Bayesian M - R data analysis, The European Physical Journal A 52, 10.1140/epja/i2016-16069-2 (2016).
  * Ayriyan  _et al._ [2021a] A. Ayriyan, D. Blaschke, A. G. Grunfeld, D. Alvarez-Castillo, H. Grigorian, and V. Abgaryan, Bayesian analysis of multimessenger M-R data with interpolated hybrid EoS, Eur. Phys. J. A 57, 318 (2021a), arXiv:2102.13485 [astro-ph.HE] .
  * Shahrbaf  _et al._ [2021] M. Shahrbaf, S. Antić, A. Ayriyan, D. Blaschke, and A. G. Grunfeld, Constraining free parameters of a color superconducting non-local Nambu-Jona-Lasinio model using Bayesian analysis of neutron stars mass and radius measurements (2021), arXiv:2105.00029 [nucl-th] .
  * Pfaff  _et al._ [2022] A. Pfaff, H. Hansen, and F. Gulminelli, Bayesian analysis of the properties of hybrid stars with the Nambu–Jona-Lasinio model, Phys. Rev. C 105, 035802 (2022).
  * Blanquier [2013] E. Blanquier,  _Le modèle de Polyakov, Nambu et Jona-Lasinio et ses applications pour décrire les particules sub-nucléaires_ , Ph.D. thesis, Toulouse III U. (2013).
  * Kobayashi  _et al._ [1971] M. Kobayashi, H. Kondo, and T. Maskawa, Symmetry Breaking of the Chiral U(3)xU(3) and the Quark Model, Progress of Theoretical Physics 45, 1955 (1971).
  * ’t Hooft [1976a] G. ’t Hooft, Computation of the quantum effects due to a four-dimensional pseudoparticle, Phys. Rev. D 14, 3432 (1976a).
  * ’t Hooft [1976b] G. ’t Hooft, Symmetry breaking through Bell-Jackiw anomalies, Phys. Rev. Lett. 37, 8 (1976b).
  * Dyson [1949] F. J. Dyson, The $s$ matrix in quantum electrodynamics, Phys. Rev. 75, 1736 (1949).
  * Faddeev [1960] L. D. Faddeev, Scattering theory for a three particle system, Zh. Eksp. Teor. Fiz. 39, 1459 (1960).
  * Hanhart and Krewald [1995] C. Hanhart and S. Krewald, Faddeev approach to the octet and decuplet baryons, Physics Letters B 344, 55 (1995).
  * Gastineau  _et al._ [2002] F. Gastineau, R. Nebauer, and J. Aichelin, Thermodynamics of the three flavor NJL model: Chiral symmetry breaking and color superconductivity, Phys. Rev. C 65, 045204 (2002), arXiv:hep-ph/0101289 .
  * Osipov  _et al._ [2004] A. A. Osipov, H. Hansen, and B. Hiller, Long distance expansion for the NJL model with SU(3) and U(1)A breaking, Nucl. Phys. A 745, 81 (2004), arXiv:hep-ph/0406112 .
  * Ebert  _et al._ [2005] D. Ebert, K. G. Klimenko, and V. L. Yudichev, Mesons and diquarks in the color neutral 2SC phase of dense cold quark matter, Phys. Rev. D 72, 056007 (2005), arXiv:hep-ph/0504218 .
  * Lawley  _et al._ [2006] S. Lawley, W. Bentz, and A. Thomas, The phases of isospin-asymmetric matter in the two-flavor NJL model, Physics Letters B 632, 495 (2006).
  * Sivia and Skilling [2006] D. Sivia and J. Skilling,  _Data Analysis: A Bayesian Tutorial_, Oxford science publications (OUP Oxford, 2006).
  * Note [1] The $\eta^{\prime}$ meson is oftentimes found to be unstable within this model and violates Eq. (15), as its mass remains quite large due to the U(1)A anomaly [19]. We will therefore make an exception for this meson and compute its mass as the real part of the complex pole of the T-matrix. The associated width is an artifact due to the lack of confinement in the NJL model.
  * Peláez [2016] J. R. Peláez, From controversy to precision on the sigma meson: A review on the status of the non-ordinary $f_{0}(500)$ resonance, Physics Reports 658, 1 (2016).
  * Celenza  _et al._ [1995] L. S. Celenza, C. M. Shakin, W.-D. Sun, J. Szweda, and X.-Q. Zhu, Quark model calculations of current correlators in the nonperturbative domain, Annals Phys. 241, 1 (1995).
  * Zyla  _et al._ [2020] P. Zyla _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2020, 083C01 (2020).
  * Dmitrasinovic  _et al._ [1992] V. Dmitrasinovic, R. Lemmer, and R. Tegen, The mass difference between the charged and neutral pions in the Nambu-Jona-Lasinio model, Physics Letters B 284, 201 (1992).
  * Rosner  _et al._ [2019] J. L. Rosner, S. Stone, and R. S. V. de Water, Leptonic decays of charged pseudoscalar mesons - 2019 (2019), arXiv:1509.02220 [hep-ph] .
  * Giusti  _et al._ [1999] L. Giusti, F. Rapuano, M. Talevi, and A. Vladikas, The QCD chiral condensate from the lattice, Nuclear Physics B 538, 249–277 (1999).
  * Dosch and Narison [1998] H. Dosch and S. Narison, Direct extraction of the chiral quark condensate and bounds on the light quark masses, Physics Letters B 417, 173–176 (1998).
  * Aoki  _et al._ [2020] S. Aoki _et al._ , FLAG review 2019, The European Physical Journal C 80, 10.1140/epjc/s10052-019-7354-7 (2020).
  * Ninomiya  _et al._ [2015] Y. Ninomiya, W. Bentz, and I. C. Cloët, Dressed Quark Mass Dependence of Pion and Kaon Form Factors, Phys. Rev. C 91, 025202 (2015), arXiv:1406.7212 [nucl-th] .
  * Bonanno and Sedrakian [2012] L. Bonanno and A. Sedrakian, Composition and stability of hybrid stars with hyperons and quark color-superconductivity, A& A 539, A16 (2012).
  * Klähn  _et al._ [2013] T. Klähn, R. Łastowiecki, and D. Blaschke, Implications of the measurement of pulsars with two solar masses for quark matter in compact stars and heavy-ion collisions: A Nambu–Jona-Lasinio model case study, Physical Review D 88, 10.1103/physrevd.88.085001 (2013).
  * Ayriyan  _et al._ [2021b] A. Ayriyan, D. Blaschke, A. G. Grunfeld, D. Alvarez-Castillo, H. Grigorian, and V. Abgaryan, Bayesian analysis of multimessenger M-R data with interpolated hybrid eos, The European Physical Journal A 57, 10.1140/epja/s10050-021-00619-0 (2021b).
  * Contrera  _et al._ [2022] G. A. Contrera, D. Blaschke, J. P. Carlomagno, A. G. Grunfeld, and S. Liebing, Quark-nuclear hybrid equation of state for neutron stars under modern observational constraints, Phys. Rev. C 105, 045808 (2022), arXiv:2201.00477 [nucl-th] .
  * Birse  _et al._ [1999] M. C. Birse, J. A. McGovern, S. Pepin, and N. R. Walet, Relativistic Faddeev approach to the NJL model at finite density (1999).
  * Pepin  _et al._ [2000] S. Pepin, M. C. Birse, J. A. McGovern, and N. R. Walet, Nucleons or diquarks: Competition between clustering and color superconductivity in quark matter, Physical Review C 61, 10.1103/physrevc.61.055209 (2000).
  * Hands  _et al._ [2006] S. Hands, S. Kim, and J.-I. Skullerud, Deconfinement in dense 2-color QCD, Eur. Phys. J. C 48, 193 (2006), arXiv:hep-lat/0604004 .
  * Abuki  _et al._ [2010] H. Abuki, G. Baym, T. Hatsuda, and N. Yamamoto, The NJL model of dense three-flavor matter with axial anomaly: the low temperature critical point and BEC-BCS diquark crossover, Phys. Rev. D 81, 125010 (2010), arXiv:1003.0408 [hep-ph] .