# Constraining the high-density behavior of nuclear symmetry energy with
direct Urca processes

Olfa Boukari<EMAIL_ADDRESS>ISEPBG-Soukra, University of
Carthage, Avenue de la République BP 77-1054 Amilcar, Tunisia  Tuhin Malik
<EMAIL_ADDRESS>CFisUC, Department of Physics, University of Coimbra, 3004-516
Coimbra, Portugal  Aziz Rabhi<EMAIL_ADDRESS>CFisUC, Department of Physics,
University of Coimbra, 3004-516 Coimbra, Portugal  LR11ES23, University of
Carthage, Avenue de la République BP 77-1054 Amilcar, Tunisia  Constança
Providência<EMAIL_ADDRESS>CFisUC, Department of Physics, University of Coimbra,
3004-516 Coimbra, Portugal

###### Abstract

The density dependence of the symmetry energy in relativistic mean-field
models with density dependent couplings is discussed in terms of the possible
opening of nucleonic direct Urca processes inside neutron stars, which induce
a very rapid cooling of the star. The modification of the parametrization of
the isospin channel of two models, DD2 and DDMEX, keeping the same isoscalar
properties is considered and the implications are discussed. Within the models
discussed it is not possible the onset of nucleonic direct Urca processes in
stars with a mass below $\sim 1.6-1.8\,M_{\odot}$. The lowest masses that
allow direct Urca processes are associated to a slope of the symmetry energy
$\sim 60$ MeV and a symmetry energy incompressibility close to zero. It is
shown that the parametrization of the isospin channel proposed destroys the
correlation between symmetry energy slope and incompressibility previously
identified in several works.

##  I Introduction

In the past decade, multi-messenger astronomy has made significant progress,
particularly with the recent detection of gravitational waves (GW) by the LIGO
and Virgo Collaboration (LVC). Notably, they detected GW originating from two
distinct binary neutron star (BNS) mergers. The first event, known as GW170817
[1, 2, 3], involved a BNS merger with a total mass of 2.7M⊙ and mass
components ranging between 1.17M⊙ and 1.6M⊙. Similarly, the subsequent
detection of the GW190425 merger, as detailed in reference [4], further
expanded our knowledge. This event featured a total mass of
$3.4^{+0.3}_{-0.1}\mathrm{M}_{\odot}$, with mass components spanning from
1.12M⊙ to 2.52M⊙. These observations provide valuable insights into high-
density stellar matter [5, 6, 7], and offer new constraints for understanding
the Equation of State (EoS) of stellar matter [8, 9, 10, 11].

The ambitious goal of the astrophysics community is to combine these multi-
messenger results to constrain the EoS of neutron stars (NSs)[5, 12, 13, 14,
15, 16, 6, 7], and infer the NS composition. The nuclear symmetry energy
constitutes an essential ingredient to explore the macroscopic properties of
NSs, and should satisfy theoretical constraints such as those obtained from
chiral effective field theory (chEFT) calculations for pure neutron matter
[17]. It should also be consistent with experimental and observational data.

Thanks to the considerable efforts deployed in astrophysics and nuclear
physics over the last two decades, significant progress has been made in
determining the symmetry energy $E_{\textrm{sym}}(\rho_{B})$, in particular
around and below the saturation density of nuclear matter $\rho_{0}$. The
symmetry energy $E_{\rm{sym}}(\rho_{B})$ above saturation density and a
possible hadron-quark phase transition are among the most uncertain parts of
the EoS of dense neutron rich matter [18, 19, 20, 21]. A precise determination
of the behavior of $E_{\rm{sym}}(\rho_{B})$ at densities beyond saturation is
essential. To probe the symmetry energy in this density range, data obtained
from observations of NSs have revealed to be particularly useful compared with
terrestrial experiments. Since the detection of GW170817, astrophysical data
have stimulated many interesting studies of the symmetry energy. Recently Bao-
Jun Cai et al. [22] have studied the behaviour of the nuclear symmetry energy
at $2-3\rho_{0}$ on the basis of its slope $L$, its curvature
$K_{\textbf{sym}}$ and its skewness $Q_{\textbf{sym}}$ at saturation
$\rho_{0}$. This investigation was achieved by developing the function
$E_{sym}(\rho_{B})$ in terms of a suitable auxiliary function, allowing an
accurate prediction at high densities of the symmetry energy.

Several attempts have been made to obtain the EOS of supernuclear matter in
NSs from different models. In particular, by adopting density-dependent DD
models such as DDRMF, DD-MEX [23], DD-LZ1 [24], massive NSs with masses around
2.33 - 2.48 M⊙ could be generated. The DD models do not include mesonic
nonlinear mixing or self-interaction terms. Instead, a density dependence of
the nucleon-meson couplings is introduced, which takes into account the
effects of the nuclear medium on the couplings, as predicted, for example, by
relativistic Dirac-Brueckner-Hartree-Fock calculations [25]. However, DD
models are generally characterized by an isovector coupling $\Gamma_{\rho}$
which tends to zero at higher densities. This has the consequence of softening
the symmetry energy at high densities, allowing large neutron densities, and
therefore these models do not predict nucleonic neutrino emission processes
inside the NS that are responsible for the fast cooling of the NS, i.e. the
nucleon direct Urca (DU) processes [26, 27]. In order to induce the DU
processes in 1.6 to 1.8 $M_{\odot}$ stars, as predicted in [28], the authors
in [29, 27] have included hyperons in the star. This results in a softer EoS,
which makes predicting stars with masses around or above 2$M_{\odot}$
difficult.

Bayesian inference studies have supported the possible existence of NSs with
$2.50$-$2.60\mathrm{M}_{\odot}$ under the constraints of NS properties of
$1.4\mathrm{M}_{\odot}$ [30]. This methodology based on microscopic models
offers as a major advantage, the possibility, once the inference is complete,
of discussing the composition of nuclear matter properties through the EOS and
the symmetry energy at high densities [31, 32, 33]. The knowledge of the
symmetric nuclear matter EOS has proved to be indispensable for studying the
$\beta$ equilibrium condition [34, 35, 36, 32]. In this context, the authors
of [33] have recently developed a Bayesian inference approach, in the
framework of a density dependent model, in order to determine how the GW and
NICER data constrain the high-density symmetry energy values. They have
confirmed that these type of models do not allow the nucleon DU processes
inside NSs. This is a characteristic of all density dependent models that
describe the nucleon-$\rho$ meson coupling as an exponential decreasing
function with density. However, NS cooling curves seem to indicate that these
kind of processes should be possible inside NSs with a mass $\gtrsim
1.6\,M_{\odot}$ [28]. In order to adapt these DD models at high densities,
modifications have recently been introduced by Malik et al. (described in Ref.
[37]) which predict nucleon DU processes inside NS. The proposed modified DD
models allow to control the slope $L$ over a large range of densities both at
high and low densities. For a different parametrization of the $\rho$-meson
coupling see also [38].

In the following we propose a set of DD models, that satisfy the set of
constraints on the symmetry energy proposed in several studies [39, 37], as
well as constraints from the neutron matter calculation within a chEFT
approach [17], which predict the onset of DU cooling processes inside the NS.
In the following, we will refer as $M_{DU}$ the minimum NS mass where nucleon
DU processes are allowed to occur. We will also discuss how some constraints
on the properties of symmetry energy, as delineated in [39, 40, 37], are
satisfied. In particular, considering the analysis of recent NS observations,
the authors of [39] have obtained at saturation density for the symmetry
energy slope $L\approx 57.7\pm 19$ MeV, the curvature $K_{sym}\approx-107\pm
88$ MeV, and at twice the saturation density the symmetry energy
$E_{sym}(2\rho_{0})\approx 51\pm 13$ MeV, all determined at a confidence level
of $68\%$. These values are consistent with previous values obtained from
other analyses based on experimental data, ab-initio calculations and NS
observations, such as [41, 42, 40]. We also consider the constraints found in
[37], based on the strong correlation between the symmetry energy slope at
$\rho=2.5\rho_{0}$ and the DU mass $M_{DU}$. From this correlation it was
found that the symmetry energy slope range at $2.5\rho_{0}$ should be in the
interval $L(2.5\rho_{0})\approx 54-48$ MeV to predict DU processes occurring
inside stars with a mass $M_{DU}\gtrsim 1.6-1.8\,M_{\odot}$.

The present study is organized as follows. In Sec. II we review the formalism
used in the analysis. In Sec. III we discuss the properties of a set of DD
models, based on the DD2 [43] and DDMEX [44] models, which predict the
occurrence of nucleon DU processes within NS by modifying the parameterization
of the isovector channel, keeping the same isoscalar description and the
symmetry energy at saturation. In Sec. IV, we generalize our study and apply a
Bayesian inference approach to study the complete isovector channel, also
modifying the symmetry energy at saturation and identifying possible
correlations. Finally, some conclusions are drawn in Sec. V.

##  II Models

In this section, we review the description of nuclear matter EoS with
relativistic mean-field (RMF) models with density dependent couplings, and the
characterization of the density dependence of the symmetry energy through its
expansion at the saturation density.

###  II.1 Field Theoretical Models with DDRMF Lagrangian

In order to study the symmetry energy at high density, we will consider RMF
models with density dependent baryon-meson couplings that avoid self-
interacting and mixed terms between mesons and which we designate by DD
models. Within this approach, we start from a Lorentz-covariant Lagrangian
density which describes baryons interacting with mesons. It is assumed the
minimal coupling between the baryons and the mesons.

In the DD model, the nucleons interact with each other by exchanging scalar-
isoscalar ($\sigma$), vector-isoscalar ($\omega^{\mu}$), and vector-isovector
($\vec{\rho}^{\mu}$) mesons. The DDRH Lagrangian density can be written as:

$\displaystyle\mathcal{L}_{DD}=$
$\displaystyle\sum_{i=p,~{}n}\overline{\psi}_{i}\left[\gamma^{\mu}\left(i\partial_{\mu}-\Gamma_{\omega}(\rho_{B})\omega_{\mu}-\frac{\Gamma_{\rho}(\rho_{B})}{2}\vec{\tau}\cdot\vec{\rho}_{\mu}\right)\right.$
$\displaystyle\left.-(M-\Gamma_{\sigma}(\rho_{B})\sigma)\right]\psi_{i}+\frac{1}{2}\left(\partial^{\mu}\sigma\partial_{\mu}\sigma-
m_{\sigma}^{2}\sigma^{2}\right)$
$\displaystyle+\frac{1}{2}m_{\omega}^{2}\omega_{\mu}\omega^{\mu}-\frac{1}{4}\Omega^{\mu\nu}\Omega_{\mu\nu}$
$\displaystyle+\frac{1}{2}m_{\rho}^{2}\vec{\rho}_{\mu}\cdot\vec{\rho}^{\mu}-\frac{1}{4}\mathbf{B}_{\mu\nu}\cdot\mathbf{B}^{\mu\nu},$
(1)

where $\psi_{i}$ represents the Dirac spinor of the nucleons, $M$ is the
nucleon mass. $\Omega_{\mu\nu}$ and $B_{\mu\nu}$ are the vector meson field
tensors, $\gamma^{\mu}$ and $\vec{\tau}$ are the Dirac and the Pauli matrices,
respectively. The $\Gamma_{\sigma}(\rho_{B})$, $\Gamma_{\omega}(\rho_{B})$ and
$\Gamma_{\rho}(\rho_{B})$ are the coupling constants of the nucleons to the
meson fields $\sigma$, $\omega$, and $\rho$ respectively, with a corresponding
meson masses are $m_{\sigma},\,m_{\omega}$ and $m_{\rho}$.

We focus on the parameterizations of DDRMF that depend on the vector density.
The nucleon-meson density-dependent coupling parameters are written in the
form:

$\Gamma_{j}(\rho_{B})=\Gamma_{j}(\rho_{0})h_{M}(x)~{},\quad
x=\rho_{B}/\rho_{0}~{},\quad j=\sigma,\ \omega,\ \rho$ (2)

where $\rho_{0}$ is the saturation density of symmetric nuclear matter. In the
present study, we consider the parameterizations DD2 [45] and DDMEX [23]. For
these two parameterizations, the coupling constants of the $\sigma$ and
$\omega$ mesons to the nucleons are given in terms of the functions $h_{j}$
[45],

$h_{j}(x)=a_{j}\frac{1+b_{j}(x+d_{j})^{2}}{1+c_{j}(x+d_{j})^{2}},$ (3)

with $x=\rho_{B}/\rho_{0}$ and $\rho_{B}$ the baryonic density. The parameters
$a_{j},b_{j},c_{j},$ and $d_{j}$ are defined in [45] for DD2 and [23] for
DDMEX. In the DD2 and DDMEX models, the nucleon DU processes do not occur in
the NS interior because the $\rho_{DU}$ density is greater than the central
density of the most massive star, $\rho_{c}$. In order to overcome this
property in [37] a generalization of the isovector channel $\rho$-meson
coupling was proposed, including a new parameter $y$ which controls the high-
density behaviour. This parametrization will be considered in the present
study:

$h_{\rho}(x)=y~{}\exp[-a_{\rho}(x-1)]+(1-y)~{},\quad 0<y\leq 1~{}.$ (4)

From the Lagrangian density, we obtain the following meson field equations in
the mean-field approximation

$\displaystyle
m_{\sigma}^{2}\sigma=\Gamma_{\sigma}(\rho_{B})\rho_{s}=\Gamma_{\sigma}(\rho_{B})\left\langle\overline{\psi}\psi\right\rangle,$
$\displaystyle
m_{\omega}^{2}\omega=\Gamma_{\omega}(\rho_{B})\rho_{B}=\Gamma_{\omega}(\rho_{B})\left\langle\psi^{{\dagger}}\psi\right\rangle,$
$\displaystyle
m_{\rho}^{2}\rho=\frac{\Gamma_{\rho}(\rho_{B})}{2}\rho_{3}=\frac{\Gamma_{\rho}(\rho_{B})}{2}\left\langle\psi^{{\dagger}}\tau_{3}\psi\right\rangle,$
(5)

where $\rho_{s}$ is the scalar density, and $\rho_{3}=\rho_{p}-\rho_{n}$ with
$\rho_{p}$ and $\rho_{n}$, respectively, the proton and the neutron densities.
The Dirac equations for the nucleons are given by

$\displaystyle\left[i\gamma^{\mu}\partial_{\mu}-\gamma^{0}\left(\Gamma_{\omega}(\rho_{B})\omega+\frac{\Gamma_{\rho}(\rho_{B})}{2}\rho\tau_{3}+\Sigma^{R}_{0}(\rho_{B})\right)\right.$
$\displaystyle\left.-M_{B}^{*}\right]\psi_{B}=0,$ (6)

where $\omega$ and $\rho$ define the time components of $\omega^{\mu}$ and of
the third component of $\vec{\rho}^{\mu}$, the nucleon isospin third
components take the values $\tau_{3}=1$ and $\tau_{3}=-1$ for protons and
neutrons, respectively, the effective mass of nucleons is given in terms of
the scalar meson $\sigma$

$\displaystyle M_{B}^{*}$ $\displaystyle=M-\Gamma_{\sigma}(\rho_{B})\sigma.$
(7)

and the rearrangement term $\Sigma_{0}^{R}$, due to the density dependence of
coupling constants [45, 46, 47], is given by

$\displaystyle\Sigma^{R}_{0}(\rho_{B})=\frac{\partial\Gamma_{\omega}(\rho_{B})}{\partial\rho_{B}}\omega\rho_{B}+\frac{1}{2}\frac{\partial\Gamma_{\rho}(\rho_{B})}{\partial\rho_{B}}\rho\rho_{3}-\frac{\partial\Gamma_{\sigma}(\rho_{B})}{\partial\rho_{B}}\sigma\rho_{s}.$
(8)

From the energy-momentum tensor in a uniform system, the energy density,
$\mathcal{E}$, and the pressure, $P$, of infinite nuclear matter can be
obtained, respectively, as

$\displaystyle\mathcal{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2}m_{\sigma}^{2}\sigma^{2}+\frac{1}{2}m_{\omega}^{2}\omega^{2}+\frac{1}{2}m_{\rho}^{2}\rho^{2}$
(9) $\displaystyle+$
$\displaystyle\sum_{i=p,n}\frac{\gamma}{2\pi^{2}}\int_{0}^{k_{Fi}}k^{2}\sqrt{k^{2}+{M_{i}^{*}}^{2}}dk,$
(10) $\displaystyle P$ $\displaystyle=$
$\displaystyle\rho_{B}\Sigma_{R}(\rho_{B})-\frac{1}{2}m_{\sigma}^{2}\sigma^{2}+\frac{1}{2}m_{\omega}^{2}\omega^{2}+\frac{1}{2}m_{\rho}^{2}\rho^{2}$
(11) $\displaystyle+$
$\displaystyle\sum_{i=p,n}\frac{\gamma}{6\pi^{2}}\int_{0}^{k_{Fi}}\frac{k^{4}dk}{\sqrt{k^{2}+{M_{i}^{*}}^{2}}},$
(12)

where $\gamma=2$ is the spin degeneracy factor, and $k_{Fi}$ is the Fermi
momentum of nucleon $i$. The binding energy per nucleon is defined by

$\frac{E_{b}}{A}=\frac{\mathcal{E}}{\rho_{B}}-M.$ (13)

###  II.2 Nuclear symmetry energy

For nuclear matter consisting of protons and neutrons at the density
$\rho_{B}$ with the isospin asymmetry parameter
$t=(\rho_{n}-\rho_{p})/\rho_{B}$, the energy per baryon(nucleon)
$E=\displaystyle\frac{\mathcal{E}}{\rho_{B}}$, can be expanded with respect to
$t$ and $\rho_{B}$:

${E}(\rho_{B},\delta)\simeq E_{SNM}(\rho_{B})+E_{sym}(\rho_{B})t^{2}+{\cal
O}(t^{4})\ ,$ (14)

where $E_{SNM}(\rho_{B})$ is the energy per particle of symmetric nuclear
matter and the symmetry energy $E_{sym}(\rho_{B})$ quantifies the energy
needed to make nuclear matter more neutron rich. Assuming the charge symmetry
of the nuclear forces the symmetry energy can be identified as a quadratic
term in $t$ of the energy per particle:

$E_{sym}(\rho_{B})=\frac{1}{2}\frac{\partial^{2}E}{\partial
t^{2}}\Big{|}_{t=0}.$ (15)

To study the density dependence in a range of densities below the quark-hadron
phase transition, it is common to characterize the density dependence of the
SNM EOS $E_{0}(\rho_{B})$ and the symmetry energy $E_{\rm{sym}}(\rho_{B})$
around the saturation density $\rho_{0}$ in terms of the parameters of the
Taylor series. In the following, we consider terms until third order on the
variable $\displaystyle x=\frac{\rho_{B}-\rho_{0}}{3\rho_{0}}$,

$\displaystyle E_{0}(\rho_{B})$ $\displaystyle=$ $\displaystyle
E_{0}(\rho_{0})+\frac{K_{0}}{2}x^{2}+\frac{Q_{0}}{6}x^{3},$ (16)
$\displaystyle E_{\rm{sym}}(\rho_{B})$ $\displaystyle=$ $\displaystyle
E_{\rm{sym}}(\rho_{0})+L_{0}\,x+\frac{K_{\rm{sym},0}}{2}\,x^{2}$ (17)
$\displaystyle+$ $\displaystyle\frac{Q_{\rm{sym},0}}{6}x^{3},$ (18)

$E_{0}(\rho_{B})\equiv E_{\rm{SNM}}(\rho_{B})$ is the energy per nucleon for
symmetric nuclear matter (SNM). The incompreensibility coefficient $K_{0}$ and
the the skewness $Q_{0}$ are, respectively, defined by:

$\displaystyle K_{0}$ $\displaystyle=$ $\displaystyle
9\rho_{0}^{2}\dfrac{\partial^{2}E}{\partial\rho_{B}^{2}}\Bigg{|}_{\rho_{B}=\rho_{0}}=9\left[\frac{dP_{0}}{d\rho_{B}}\right]_{\rho_{B}=\rho_{0}}$
$\displaystyle Q_{0}$ $\displaystyle=$ $\displaystyle
27\rho_{0}^{3}\dfrac{\partial^{3}E}{\partial\rho_{B}^{3}}\Bigg{|}_{\rho_{B}=\rho_{0}}$

with $P_{0}$ the SNM pressure. $E_{\rm{sym}}(\rho_{0})$ is the symmetry energy
at saturation and the quantities $L_{0}$, $K_{sym}$ and $Q_{sym}$ are its
slope, curvature and skewness respectively, at saturation:

$\begin{array}[]{c}\displaystyle{L_{0}=3\rho_{0}\frac{\partial
E_{sym}(\rho_{B})}{\partial\rho_{B}}\Big{|}_{\rho_{B}=\rho_{0}}}\ ,\\\
\displaystyle{K_{sym,0}=9\rho_{0}^{2}\frac{\partial^{2}E_{sym}(\rho_{B})}{\partial\rho_{B}^{2}}\Big{|}_{\rho_{B}=\rho_{0}}}\\\
\ \
\displaystyle{Q_{sym,0}=27\rho_{0}^{3}\frac{\partial^{3}E_{sym}(\rho_{B})}{\partial\rho_{B}^{3}}\Big{|}_{\rho_{B}=\rho_{0}}}.\end{array}$
(20)

This is considered to be a good approximation even for small proton fraction
$y_{p}$ [48]. These parameters characterize the density dependence of nuclear
symmetry energy around normal nuclear matter density and thus provide
important information on the behavior of nuclear symmetry energy at both high
and low densities. Also, the curvature parameter $K_{sym,0}$ would be a
significant measurement which distinguishes the different parametrizations.
The shift of the incompressibility with asymmetry, is given by

$K_{asy,0}=K_{sym,0}-6L_{0}-\frac{Q_{0}}{K_{0}}L_{0}\approx K_{sym,0}-6L_{0}.$
(21)

This value can be correlated to experimental observations of the giant
monopole resonance (GMR) of neutron-rich nuclei [49] such as on even-$A$ Sn
isotopes with a value of $K_{asy,0}=-550\pm 100$ MeV, see discussion [50, 51,
52].

##  III Results and discussion

The parameterisation of the $\rho$ meson coupling proposed in Eq. (4) allows
us to obtain a set of modified DD models consistent with both low and high
density constraints. We obtain for a given symmetry energy at saturation
density a range of different allowed isovector properties, in particular the
symmetry energy slope $L_{0}$ and its curvature $K_{sym,0}$, for the same
isoscalar properties i.e incompressibility, binding energy and effective mass
of the baryons as shown in Table 1.

The proton fraction and other NS characteristics which are sensitive to the
high density behavior of the symmetry energy are used to constrain the new
parameter, $y$, in Eq. (4). In particular, the onset of the nucleon DU
processes as discussed in [28] is an observation that constrains $y$, see
[37].

We build a set of EOS based on the modified DD2 and DDMEX models that
essentially satisfy neutron matter chiral effective field theory (chEFT) [17,
53] constraints within $2\sigma$, as shown in Fig. 1 top panels, and verify
the constraints on the symmetry energy at saturation density referred in the
Introduction [39, 41, 42, 37]: a slope of nuclear symmetry energy at
saturation density $\rho_{0}$ satisfying $L_{0}\approx 57.7\pm 19$ MeV at 68%
confidence level and a symmetry energy at twice $\rho_{0}$,
$E_{\rm{sym}}(2\rho_{0})\approx 51\pm 13$ MeV at 68% confidence level [39];

the symmetry energy and respective slope satisfying
$E_{sym}(\rho_{0})=(31.7\pm 3.2)$ MeV and $L_{0}=(58.7\pm 28.1)$ MeV [40]; the
slope $L$ for a density $\approx 2.5\rho_{0}$ taking a value in the range
$48-54$ MeV for the nucleon DU processes to occur inside neutron stars [37].

In order to constraint the high density behavior of the symmetry energy, the
parameters $(y,a_{\rho})$ were chosen such that a given $M_{DU}$ is obtained,
i.e. 1.6, 1.8, 2.0 $M_{\odot}$. Acceptable models should also satisfy the low
density constraints imposed by chEFT neutron matter calculations. In Table 2,
some symmetric nuclear matter properties calculated at saturation density and
at $2.5\rho_{0}$ are given for all proposed EOS. As referred before, by
construction the two families of models DD2 and DDMEX have the symmetry energy
at saturation $E_{sym,0}$ of the original DD2 and DDMEX, given in Table 1.

By analyzing Table 2, it is seen that $L_{0}$ is correlated (anti-correlated)
with $y$ ($a_{\rho}$) while $L(2.5\rho_{0})$ are anti-correlated (correlated)
with $y$ ($a_{\rho}$): The combined effect of the two parameters is to
increase the slope of $L$ at high densities, preventing the symmetry energy
from softening. On the other hand, $K_{sym,0}$ is correlated with $a_{\rho}$
and anti-correlated with $y$: it increases when $a_{\rho}$ increases and
remains negative for DD2, while it can reach positive values for DDMEX.
Recently, positive values for $K_{sym,0}$ have been obtained in [54] in an
attempt to simultaneously describe the results of PREX2 [55] and CREX [56]
within a model that included the $\delta$-meson. However, previous analysis
seem to prefer negative values: in [39], the authors predict
$K_{sym,0}\equiv-107\pm 88$ MeV at 68% confidence level, considering a large
set observational data, a value consistent with other analyses that also
considered experimental data -100$\pm 100$ MeV [57] and -112$\pm 71$ MeV [58].

In general, the effect of including the parameter $y$ is to increase the
symmetry energy slope $L_{0}$ at saturation. However, an adequate choice of
both parameters $(y,\,a_{\rho})$ also allows to fix $L_{0}$ with a value
smaller than the one of DD2, but with a much larger value of $K_{sym,0}$. For
the DDMEX family we have only obtained EOS with $L_{0}>49$ MeV, the value that
characterizes DDMEX.

In the Fig. 1 (top panels), the pure neutron matter (PNM) chEFT EOS from [17]
is included considering $1\sigma$ (dark gray) and $2\sigma$ (light gray), as
well as the PNM EOS of the DD modified models for $M_{DU}$ equal to 1.6, 1.8
and 2.0 $M_{\odot}$, the DD2 family in the left panel and the DDMEX in the
right one. The constraints derived for PNM exclude DU processes inside stars
with a mass $\lesssim$ 1.6 $M_{\odot}$ at $2\sigma$. According to [28], the NS
cooling curves seem to indicate that MDU $\sim 1.6-1.8$ M⊙. Within the DDMEX
family we are not able to go below M${}_{\rm DU}=1.8M_{\odot}$ and
simultaneously satisfy chEFT constraints. It is seen that within the two
families and the parametrization defined in Eq. (4) for the $\rho$-meson
coupling, the chEFT constraints of [34] are only satisfied in the complete
density range within the 2$\sigma$ interval. Introducing the $y$ contribution
makes the symmetry energy harder above saturation density, as expected since
this term does not allow the $\rho$-meson coupling to converge to zero, see
the middle panels. The bottom panels show the symmetry energy slope as a
function of the density and the effect of the $y$ parameter is easily
identified giving rise to a kind of plateau before a soft decrease occurs at
the larger densities.

Table 1:  The properties of symmetric nuclear matter at saturation density for the models under study: the nuclear saturation density $\rho_{0}$, the binding energy per particle $B/A$, the incompressibility $K_{0}$, the skewness coefficient $Q_{0}$, the symmetry energy $E_{sym,0}$, the slope of the symmetry energy $L_{0}$, and the effective nucleon mass $M^{*}$. All quantities are in MeV, except $\rho_{0}$ which is in fm-3, and the effective nucleon mass is normalized to the nucleon mass. Model | $\rho_{0}$ | $B/A$ | $K_{0}$ | $E_{sym,0}$ | $M^{*}/M$ | $Q_{0}$
---|---|---|---|---|---|---
DD2 Family | 0.15 | -16.03 | 243 | 30.8 | 0.55 | 169
DDMEX Family | 0.15 | -16.09 | 267 | 32.3 | 0.55 | 877
Table 2: The symmetric nuclear matter properties at saturation density for the models under study: The slope of the symmetry energy $L_{0}$ and the curvature $K_{sym,0}$ and $K_{asy,0}$ at saturation density, the slope of the symmetry energy $L$ at $2.5$ saturation density for the models under study for a set of ($a_{\rho}$,y). Model | $a_{\rho}$ | y | $L_{0}$ | $L(2.5\rho_{0})$ | $K_{sym,0}$ | $K_{asy,0}$
---|---|---|---|---|---|---
|  |  | MeV
DD2 family: |  |  |  |  |
DD2 | 0.5 | 1.0 | 55.12 | 30.27 | -93.21 | -415.93
DD2-2a | 0.8 | 0.51 | 63.85 | 43.16 | -52.18 | -429.84
DD2-2b | 1.0 | 0.47 | 59.38 | 45.53 | -34.93 | -386.78
DD2-2c | 1.2 | 0.45 | 53.59 | 46.60 | -9.62 | -328.576
DD2-1.8a | 1.05 | 0.4 | 62.18 | 48.60 | -28.76 | -398.01
DD2-1.8b | 1.2 | 0.38 | 59.55 | 50.05 | -13.45 | -367.96
DD2-1.6a | 1.25 | 0.29 | 66.33 | 55.45 | -9.636 | -405.16
DD2-1.6b | 1.3 | 0.3 | 64.50 | 54.95 | -5.078 | -406.86
DDMEX family: |  |  |  |  |
DDMEX | 0.6 | 1.0 | 49.66 | 29.87 | -71.57 | -369.17
DDMEX-2a | 1.0 | 0.42 | 65.41 | 48.27 | -8.08 | -410
DDMEX-2b | 1.1 | 0.4 | 63.85 | 49.60 | -8.09 | -390.84
DDMEX-2c | 1.2 | 0.41 | 59.75 | 49.34 | 3.92 | -354.25
DDMEX-1.8a | 1.2 | 0.33 | 67.31 | 54.18 | 1.92 | -401.68
DDMEX-1.8b | 1.3 | 0.3 | 67.80 | 56.38 | 11.1 | -395.4
DDMEX-1.6 | 1.3 | 0.27 | 70.86 | 58.53 | 10.78 | -411.88
|
---|---
|
|
Figure 1: (Color online) Upper panels show the pure neutron matter pressure as
a function of the baryonic density for DD2 (left) and DDMEX (right) families
under consideration varying ($a_{\rho},y$) to obtain a given MDU i.e. 1.6,
1.8, 2.0 M⊙ and including the 1$\sigma$ and 2$\sigma$ bands from chiral
effective field theoretical calculations [17]. Middle panels reveal the
symmetry energy and in the bottom its slope as function of the baryon density
for DD2 (left) and DDMEX (right) families.

By varying $(y,a_{\rho})$ it is possible to increase the slope of the symmetry
energy at $2.5\rho_{0}$ with respect to the original models, DD2 and DDMEX:
$L(2.5\rho_{0}$) varies in the range [43,55] MeV for DD2 family and [48,58]
MeV for DDMEX family. This intervals are reasonably compatible with
$L(2.5\rho_{0})\equiv 48-54$ MeV, taking into account that they have been
deduced from the strong correlation existing between $L$ and MDU [37],

$\frac{L}{\rm MeV}=\frac{-31.224}{\rm MeV}\frac{M_{\rm
DU}}{M_{\odot}}+\frac{104.339}{\rm MeV}$ (22)

with a Pearson coefficient of 0.9.

Figure 2: $K_{sym,0}$ versus $L_{0}$. Original DD2 (black solid line) and
DDMEX (black dashed line) varying $y$, modified DD2 (colored solid line) and
DDMEX (colored dashed line) using different values of the parameter $a_{\rho}$
and varying $y$.

DDMEX model presents a softer EoS with a lower slope L${}_{0}=49.7$MeV that
satisfies the 1$\sigma$ constraint as shown in Fig.1. We notice, the higher
the slope the stiffer is the pressure. The DDMEX-1.6 model has
L${}_{0}=70.9$MeV and misses the chiral EFT above $0.13$fm-3.

In Fig. 2, $K_{\rm sym,0}$ is plotted versus $L_{0}$ changing the parameter
$y$ between 0 and 1, and considering several values for the $a_{\rho}$
parameter. A negative correlation between $L_{0}$ and $K_{\rm sym,0}$ is
obtained when $a_{\rho}$ takes the original values of models DD2[43] and DDMEX
[44]. A similar correlation was obtained in [59, 60, 61], where it was shown
that the coefficients $L_{0}$ and $K_{\rm sym,0}$ taken at saturation density
for a large set of RMF and Skyrme force models present a linear anti-
correlation between them. However, for the modified EOS with a larger
parameter $a_{\rho}$, the $L_{0}$-$K_{\rm sym,0}$ variation is not monotonic,
showing that the $L_{0}$-$K_{\rm sym,0}$ correlation is washed out for
sufficiently large values of the parameter $a_{\rho}$.

###  III.1 Direct Urca process with modified DD models

In this section, we investigate the opening density of the direct DU process
for the modified density-dependent model families. The nucleonic DU processes
[62] by neutrino emission are described by the equations :

$n\rightarrow p+e^{-}+\bar{\nu}_{e},\qquad\textrm{and}\qquad
p+e^{-}\rightarrow n+\nu_{e}.$ (23)

This process is considered to be the most efficient NS cooling mechanism. It
is only triggered if there is conservation of momentum, i.e. $p_{{\rm F}n}\leq
p_{{\rm F}p}+p_{{\rm F}e}$, where $p_{\rm Fi}$ is the Fermi momentum of the
particle species $i$.

Therefore, the proton fraction needs to equal or exceed a minimum proton
fraction $y_{p}^{\rm min}$ [63], to allow the DU process to operate:

$y_{p}^{\rm min}=\frac{1}{1+\left(1+x_{e}^{1/3}\right)^{3}},$ (24)

with $x_{e}=\rho_{e}/\left(\rho_{e}+\rho_{\mu}\right)$, and $\rho_{e}$ and
$\rho_{\mu}$ the electron and muon densities. In the following, $\rho_{\rm
DU}$ denotes the baryon density at which the DU process starts to operate and
$M_{\rm DU}$ is the mass of the star with a central density equal to
$\rho_{\rm DU}$.

Figure 3: (Color online) The DU onset for a given MDU, i.e 1.6, 1.8, 2.0
$M_{\odot}$ (blue, green, red), for the DD2 family (left) and the DDMEX family
(right).

In previous studies [64, 59], the authors investigated how the symmetry energy
in various non-linear RMF models affects the DU process, emphasizing its
connection to the density dependence of symmetry energy. They discovered that
this dependence notably influences this process. In cases like DD2 and DDMEX
models, nucleonic DU processes are prohibited, as in most conventional DD
models. However, fast cooling can still occur with hyperon onset. In such
cases, hyperonic DU processes, involving hyperon decay into other hyperons or
nucleons with neutrino emission, begin just above hyperon onset densities
[26].

To explore the impact of the density dependence of the symmetry energy slope
$L$ on the DU processes, we depict in Fig. 3 the onset density of nuclear DU
processes

the DD modified models discussed in Table II to illustrate the relationship
between baryon density and the ratio of proton fractions, $y_{p}^{\rm
min}/y_{p}$, where $y_{p}$ is the $\beta-$equilibrium proton fraction. The
onset of the DU process occurs when this ratio, $y_{p}^{\rm min}/y_{p}$, is
less or equal to 1. By analyzing the onset of the DU process at $\rho_{dU}$,
we can classify three sets of EOS based on similar densities: $\sim
2.5\rho_{0}$, $\sim 2.75\rho_{0}$ and $\sim 3.1\rho_{0}$ for $M_{DU}$ values
of 1.6, 1.8, and 2.0$M_{\odot}$, respectively. As expected, the highest values
of $M_{DU}$ correspond to the greatest values of $\rho_{DU}$, as shown in
Table 3. Since the symmetry energy at saturation remains constant for both
families, it is the parameter $L$ that governs the onset of DU processes. A
larger value of $L$ results in a lower onset density, as a larger $L$
facilitates faster growth of symmetry energy beyond the saturation density.

###  III.2 Symmetry energy and Mass-Radius

The mass-radius relationship of a static NS can be determined by solving the
Tolman-Oppenheimer-Volkoff (TOV) equation [65, 66], which incorporates the
equation of state (EOS) of the NS, representing pressure as a function of
energy density. We have included the Baym-Pethick-Sutherland outer crust [67],
and the inner crust was treated as in [33]: the outer crust is matched to the
outer core EoS at $0.04$fm-3 using a polytrope. The matching density was
chosen in order to reduce the error introduced, which was shown to be below
100 m [68]. Table 3 provides details on the maximum gravitational mass, its
baryonic mass and radius, the mass $M_{\rm DU}$ and its central density
$\rho_{\rm DU}$, as well as the radius and central density of stars with
masses of 1.6, 1.8, and 2.0 $M_{\odot}$. While the properties of the maximum
mass star remain largely unchanged with the modifications introduced in the
isovector channel, there is a noticeable impact on the radius, central
density, DU mass and radius of stars with masses less than or equal to
2.0$M_{\odot}$. Specifically, a larger slope $L$ results in smaller $M_{\rm
DU}$ and density $\rho_{DU}$, leading to reduced NS radii and central
densities.

In Fig.4, we plot the mass-radius relation ($M-R$) for the modified DD2 and
DDMEX models. Also included are some recent astrophysical observations (see
the caption for details). As referred before, it is clearly seen that the
maximum masses of NSs calculated using DD2 and DDMEX models are not very
sensitive to the isovector channel, suffering a small increase not larger than
0.02$M_{\odot}$. Notably, the DDMEX models, also the generalized
parametrizations, demonstrate the capability to predict NSs exceeding
2.5$M_{\odot}$, consistent with the observed mass range of the secondary
compact object in GW190814, reported as $2.50-2.67M_{\odot}$ (see [44]). It is
clearly seen that the effect of $y$ and new values of $a_{\rho}$ that allow DU
processes to occur inside NS is not negligible: the radius of medium and high
mass stars is 150-300 m larger due the much stiffer symmetry energy. In this
figure we have also included the results of a Bayesian inference calculation
that will be discussed below.

Table 3: The NS properties are given: the maximum mass $M_{\rm max}$ and the DU mass $M_{DU}$, respective baryon mass $M_{\rm B,max}$ and radius $R_{\rm max}$, the radius of 1.4, 1.6, 1.8 and 2.0$M_{\odot}$ stars, the DU onset baryon density $\rho_{DU}$ and the central baryon density $\rho_{c}$. Model | $M_{\rm max}$ | $M_{DU}$ | $M_{\rm B,max}$ | $R_{\rm max}$ | $R_{1.6}$ | $R_{1.8}$ | $R_{2.0}$ | $\rho_{\rm DU}$ | $\rho_{c}$
---|---|---|---|---|---|---|---|---|---
$M_{\odot}$ | km | fm-3
DD2 | 2.42 | – | 2.92 | 11.87 | 13.22 | 13.21 | 13.20 | – | 0.852
DD2-1.6a | 2.44 | 1.6 | 2.94 | 12.07 | 13.54 | 13.51 | 13.42 | $0.366$ | 0.827
DD2-1.6b | 2.44 | 1.6 | 2.94 | 12.06 | 13.50 | 13.48 | 13.40 | $0.37$ | 0.828
DD2-1.8a | 2.44 | 1.8 | 2.94 | 12.02 | 13.44 | 13.43 | 13.33 | $0.417$ | 0.833
DD2-1.8b | 2.44 | 1.8 | 2.94 | 12.02 | 13.40 | 13.39 | 13.31 | $0.415$ | 0.832
DD2-2a | 2.43 | 2.0 | 2.93 | 11.99 | 13.42 | 13.40 | 13.31 | $0.466$ | 0.838
DD2-2b | 2.44 | 2.0 | 2.94 | 11.98 | 13.36 | 13.35 | 13.27 | $0.472$ | 0.836
DD2-2c | 2.44 | 2.0 | 2.94 | 11.97 | 13.28 | 13.28 | 13.22 | $0.464$ | 0.834
DDMEX | 2.56 | – | 3.11 | 12.35 | 13.46 | 13.51 | 13.53 | – | 0.776
DDMEX-1.6 | 2.58 | 1.6 | 3.13 | 12.56 | 13.87 | 13.90 | 13.87 | $0.343$ | 0.756
DDMEX-1.8a | 2.58 | 1.8 | 3.13 | 12.53 | 13.80 | 13.84 | 13.83 | $0.366$ | 0.758
DDMEX-1.8b | 2.58 | 1.8 | 3.13 | 12.54 | 13.81 | 13.86 | 13.83 | $0.357$ | 0.758
DDMEX-2a | 2.57 | 2.0 | 3.12 | 12.49 | 13.76 | 13.80 | 13.78 | $0.40$ | 0.763
DDMEX-2b | 2.57 | 2.0 | 3.12 | 12.49 | 13.74 | 13.78 | 13.77 | $0.397$ | 0.763
DDMEX-2c | 2.57 | 2.0 | 3.13 | 12.47 | 13.67 | 13.73 | 13.72 | $0.403$ | 0.763

In addition to the constraints provided by observations of massive NSs such as
PSR J1614-2230, PSR J034+0432, and PSR J0740+6620, recent measurements from
NICER have simultaneously determined the mass and radius of a NS in the
intermediate mass region, PSR J0030+0451. According to [69] its mass is
1.44${}_{-0.14}^{+0.15}$M⊙ with a radius of $13.02_{-1.06}^{+1.24}$km, while
[70] reports a mass of $1.34_{-0.16}^{+0.15}~{}M_{\odot}$ with a radius of
$12.71_{-1.19}^{+1.14}$ km. Across various models, as verified in previous
studies [71, 64, 72], it is noted that a larger slope corresponds to larger
radii for NSs with masses less than or equal to 2.0M⊙. Specifically, for 1.4M⊙
stars, the radius increases from 13 km to 13.5 km.

Figure 4: Neutrons star masses versus the radius for the two DD families: left
panel DD2 set, right panel DDMEX set. NS masses as a function of radius for
the DD2 (blue) and DDMEX (orange) families, allowing the three couplings that
define the isovector channel to vary, i.e.
$\Gamma_{\rho}(\rho_{0}),\,a_{\rho},\,y$. Both the median and the 90% CI
probability regions are plotted. The LIGO-Virgo collaboration and NICER
observations are also shown: The grey lines represent the constraints derived
from the binary components of GW170817, along with their corresponding 90% and
50% CI [3]. Also shown are the $1\sigma$ (68%) CI for the 2D posterior
distribution in the mass-radius domain for the millisecond pulsars PSR J0030 +
0451 (cyan and violet) [70, 69] and PSR J0740 + 6620 (purple and peach) [73,
5] from the NICER X-ray data.

###  III.3 The isovector channel: Bayesian inference

In the previous section we discussed the main properties of the DD2 and DDMEX
families, keeping fixed the isoscalar channel describing the symmetric nuclear
matter EOS and the symmetry energy at saturation: the variation of the
parameters $y$ and $a_{\rho}$ has allowed to vary the slope and the curvature
of the symmetry energy at saturation. As discussed, a direct consequence was
the effect these two parameters have on the onset of the nucleonic DU
processes, allowing a decrease in the onset density and consequently a
decrease in the mass of the NS with a central density coinciding with the DU
onset density.

In the following we present a Bayesian analysis to constrain the complete
isovector channel, i.e. the couplings $\Gamma_{\rho}(\rho_{0}),\,a_{\rho}$ and
$y$, considering as constraints the onset of nucleonic DU and the low density
chEFT neutron matter pressure. The corner plot shown in Fig. 5 summarises the
main results: The bidimensional distributions relate the isovector channel
parameters to the isovector nuclear matter properties $J_{\rm sym,0}\equiv
E_{\rm sym,0}$, $L_{0}$, $K_{\rm sym,0}$, and to the NS properties maximum
mass $M_{\rm max}$, radius for 1.4 $M_{\odot}$ NS $R_{1.4}$, onset density of
the DU processes $\rho_{\rm DU}$ and corresponding NS mass with this density
at the centre $M_{\rm DU}$. The two models previously discussed, DD2 and
DDMEX, were considered as a starting point, in particular keeping all the
parameters defining the isoscalar channel fixed. The results for the two
models are shown in different colours, pink-red for DD2 and light-blue for
DDMEX.

---
Figure 5: Corner plot obtained by varying the parameters responsible for the
isovector channel, $\Gamma_{\rho_{0}}$, $a_{\rho}$ and $y$, using the
parameterisation of DD2 (pink-red) and DDMEX (light-blue) for the isoscalar
channel. Isovector nuclear matter properties ($J_{\rm sym,0}\equiv E_{\rm
sym,0}$, $L_{\rm sym,0}\equiv L_{0}$, $K_{\rm sym,0}$) and NS properties (the
maximum mass $M_{\rm max}$, the radius for 1. 4 $M_{\odot}$ NS $R_{1.4}$, the
onset density of the DU processes $\rho_{\rm DU}$ and the corresponding NS
mass with this density at the center $M_{\rm DU}$) are shown. The different
line styles in the two-dimensional distribution (solid, dashed and dotted)
correspond to 68%, 95% and 99% CI. Figure 6: The symmetry energy, denoted as
$S(\rho)\equiv E_{\rm sym}(\rho_{B})$, is plotted against the number density
$\rho$. Both $S(\rho)$ and $\rho$ are measured in MeV and fm-3 respectively,
across various isovector variations of the DD2 and DDMEX models. The median
curve is depicted together with the 90% confidence interval.

The main conclusions that can be drawn are: 1) the parameter $y$ correlates
not only with the onset density of DU, and corresponding mass, but also with
the maximum star mass. The smaller $y$ the smaller the onset density of DU and
the larger the maximum star mass. Considering the constraints imposed it was
not possible to get values for the mass $M_{DU}$ below $\sim 2M_{\odot}$ for
the DDMEX family, but it is possible to go down to $\sim 1.6\,M_{\odot}$ for
the DD2 family at 99% CI. The parameter $y$ also affects the properties
$L_{0}$, $K_{\rm sym,0}$, in particular, a smaller value of $y$ pushes these
two quantities to larger values; 2) $a_{\rho}$ correlates strongly with the
curvature $K_{sym,0}$, giving larger values, that can even be positive, for
larger values of $a_{\rho}$. $K_{sym,0}$ is in average larger for the DDMEX
family, and at 68% CI it can take values of the order $\sim+30$ MeV. The
largest values of $M_{max}$ occur for the largest values of $a_{\rho}$; 3)
there is a clear correlation between $L_{0}$ and $R_{1.4}$ as shown in other
studies when the isoscalar channel is kept fixed [71, 64, 72, 74]. Notice,
however, that if the isoscalar channel is also allowed to vary this
correlation seems to disappear [75, 76]. For both models chEFT constraints
impose $40\lesssim L_{0}\lesssim 60$ MeV; 4) The coupling $\Gamma_{\rho_{0}}$,
is totally correlated with $E_{sym,0}$ by definition, and also correlated with
$K_{sym,0}$ the smaller values of this properties being associated with the
larger values of the coupling $\Gamma_{\rho_{0}}$; 5) the isovector channel
does not affect much the NS maximum mass, a maximum variation of
0.1$M_{\odot}$ is obtained. It is also observed that the DDMEX family predicts
masses $\sim 0.2\,M_{\odot}$ larger than the DD2 family; 6) finally let us
also point down the behavior of the families concerning the radius of the
1.4$M\odot$, $R_{1.4}$: in average DD2 predicts radii $\sim 400$m smaller than
DDMEX family; 7) notice that no clear correlation between $L_{0}$ and $K_{\rm
sym,0}$ is obtained, as discussed in the previous sections.

The comparison of the above results for the NS properties with the previous
discussion allows some interesting conclusions: the prediction of the MDU
obtained are compatible with the values calculated in Table 3, the lowest
$M_{\rm DU}$ obtained was 1.6 M⊙ for the DD2 family with a slope $\rm L_{0}$
greater than $60$ MeV with $68\%$CI and $\rm K_{sym,0}$ slightly lower than 0,
see DD2-1.6a, DD2-1.6b models in Tables 2 and 3. The models DD2-1.8a, DD2-1.8b
with a $M_{\rm DU}=1.8M_{\odot}$ have a slope below $60$ MeV. For DD2-2b,
DD2-2c with $M_{\rm DU}=2.0M_{\odot}$ with a slope around $55$MeV (for
$y>0.4$, $a_{\rho}$ around 1.0) and a negative value of $K_{\rm sym,0}\lesssim
0$ and the density $\rho_{\rm DU}$ is between $0.37$ fm-3 and $0.48$ fm-3.
Note, however, that DD2-2a has a slope above 60 MeV because $y$ takes a value
above 0.5 contrary to all the other DD2 parametrizations. In summary, the
smallest $M_{\rm DU}$ requires the largest values of $L_{0}$ and values of
$K_{\rm sym,0}$ close to zero. However, large values of $M_{\rm DU}$ may also
have large values of $L_{0}$, but in this case the $K_{\rm sym,0}$ parameter
is more negative: the larger $L_{0}$ the more negative $K_{\rm sym,0}$.

For reference we show the 90% CI distribution of the symmetry energy in Fig.
6. Also included, for comparison are the respective curves for DD2 and DDMEX.
As expected, the DU constrain imposes a much harder symmetry energy at high
densities.

In Fig. 4, the mass-radius 90% CI distribution resulting from the inference
calculation is plotted for both DD2 and DDMEX families. These distributions
have a large overlap for $M\lesssim 1.4\,M_{\odot}$, because this is the
region that is most affected by the chEFT constraints. Both models are
compatible with the present observations due to the large uncertainties still
associated with the measurements. Note that some of the parametrizations
proposed in Table 3 lie outside the 90% credible level distributions because
only extreme parametrizations could satisfy the conditions necessary to allow
DU processes inside the NS.

##  IV Conclusion

In order to study the properties of NSs, in this paper, we considered two
families of density-dependent relativistic mean-field models with stiff EOS,
namely the DD2 and DDMEX families. These models are able to generate NSs with
masses as large as 2.44 M⊙ (DD2) and 2.57 M⊙ (DDMEX). In particular, DDMEX was
constrained by the low mass compact object in the event GW190814, with a mass
of 2.50-2.67 M⊙ [44]. A common behaviour of many DD models is that they do not
predict the nucleonic DU processes inside nucleonic NSs ([26, 27]), due to the
exponential decrease of the $\rho$-meson coupling with density. To explain the
cooling curves of the thermal evolution of isolated non-magnetized and non-
rotating NSs or accreting NSs with these models, [29, 27] found it necessary
to include hyperons inside the star. This may have the consequence of making
the EOS too soft, and therefore not allowing for the existence of 2$M_{\odot}$
mass stars, if the hadronic EOS is soft.

To overcome this limitation, we adopt a different approach and introduce a
generalization of the $\rho$-meson coupling constrained by observations of NS
cooling as shown in [37], allowing for the nucleonic DU processes to occur
inside NS. The objective of the present work was to investigate the
consequence of this new parameter on the nuclear matter properties and NS
properties.

By adopting as a high-density constraint, the prediction of nucleonic DU
processes inside NS, we have generated a set of DD models that obey the
$\chi$EFT EOS neutron matter constraints and satisfy nuclear matter saturation
properties, namely saturation density, binding energy per particle,
incomprehensibility and symmetry energy, according to the presently accepted
values. For the isospin channel at low densities, we have considered an
interval compatible with $\chi$EFT calculations. In addition, the symmetry
energy of the models was also constrained in order to produce a maximum star
mass of at least 2 $M_{\odot}$. Notice, however, that the isospin channel does
not have a large influence on the maximum mass and we did not obtain deviation
larger than 0.1$M_{\odot}$ from the original maximum masses predicted by DD2
and DDMEX.

The models we obtained allowed us to predict an onset of DU cooling processes
inside stars with a mass $\gtrsim$1.6 $M_{\odot}$ ($\gtrsim$1.8 $M_{\odot}$)
for the DD2 (DDMEX) family with maximum masses and radii slightly larger than
the maximum M-R of the original DD models ($2.42M_{\odot}$ for DD2 and
$2.56M_{\odot}$ for DDMEX) and could still be compatible with the low-mass
object of the binary merger at the origin of GW190814, 2.50-2.67 $M_{\odot}$
[77]. We note that considering the isoscalar channels of the models DD2 and
DDMEX, the chEFT constraints excludes the appearance of the DU processes
inside NS of masses below 1.6 $M_{\odot}$. Within these models, if stars with
masses below 1.6$M_{\odot}$ show a fast cooling it will be necessary to
include hyperons, which could be envisaged in future work.

We have shown that to decrease the $M_{\rm DU}$, it is necessary to increase
$L_{0}$. A direct consequence on the NS properties of increasing $L_{0}$
keeping the properties of the isoscalar channel fixed is the increase of the
NS radius, the smaller the mass the larger the effect. At the maximum mass the
effect was small but there was a clear reduction of the central density
$\rho_{c}$.

It was shown that in order to allow for DU nucleonic processes inside NS
described by these two families of EOS, $K_{\rm sym,0}$ may take values close
to zero or positive.

In the literature it was often pointed out a correlation between $L_{0}$ and
$K_{\rm sym,0}$ [59, 60, 61]. However, for the modified DD EOS verifying the
high density constraints imposed by the DU process, the correlation
$L_{0}$-$K_{\rm sym,0}$ is washed out for both families of EOS. It is
interesting to understand which are the effects on the low density region,
below the saturation density, of considering the $L_{0}$-$K_{\rm sym,0}$
coefficients determined from the high density constraints. These affect the
densities that define the crust of the NS. More specifically, we are
interested in studying the impact of the high-density constrained DD models on
the critical and crust-core transitions properties.

## Acknowledgments

This research received partial funding from national sources through FCT
(Fundação para a Ciência e a Tecnologia, I.P, Portugal) for projects
UIDB/04564/2020 and UIDP/04564/2020, identified by DOI
10.54499/UIDB/04564/2020 and 10.54499/UIDP/04564/2020, respectively, as well
as for project 2022.06460.PTDC with DOI 10.54499/2022.06460.PTDC. The authors
acknowledge the Laboratory for Advanced Computing at the University of Coimbra
for providing HPC resources that have contributed to the research results
reported within this paper, URL: https://www.uc.pt/lca.

## References

  * Abbott  _et al._ [2017a] B. P. Abbott, R. Abbott, T. D. Abbott, and et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017a).
  * Abbott  _et al._ [2017b] B. P. Abbott, R. Abbott, T. D. Abbott, and et al., The Astrophysical Journal Letters 848, L12 (2017b).
  * Abbott  _et al._ [2018] B. P. Abbott, R. Abbott, T. D. Abbott, and  et al. (The LIGO Scientific Collaboration and the Virgo Collaboration), Phys. Rev. Lett. 121, 161101 (2018).
  * Abbott  _et al._ [2020a] B. P. Abbott, R. Abbott, T. D. Abbott, and et al., The Astrophysical Journal Letters 892, L3 (2020a).
  * Miller  _et al._ [2021] M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, W. C. G. Ho, J. M. Lattimer, M. Loewenstein, S. M. Morsink, P. S. Ray, M. T. Wolff, C. L. Baker, T. Cazeau, S. Manthripragada, C. B. Markwardt, T. Okajima, S. Pollard, I. Cognard, H. T. Cromartie, E. Fonseca, L. Guillemot, M. Kerr, A. Parthasarathy, T. T. Pennucci, S. Ransom, and I. Stairs, The Astrophysical Journal Letters 918, L28 (2021).
  * Annala  _et al._ [2022] E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, and A. Vuorinen, Phys. Rev. X 12, 011058 (2022).
  * Biswas [2022] B. Biswas, The Astrophysical Journal 926, 75 (2022).
  * Lattimer [2012] J. M. Lattimer, Annual Review of Nuclear and Particle Science 62, 485 (2012), https://doi.org/10.1146/annurev-nucl-102711-095018 .
  * Oertel  _et al._ [2017a] M. Oertel, M. Hempel, T. Klähn, and S. Typel, Rev. Mod. Phys. 89, 015007 (2017a).
  * Tolos and Fabbietti [2020] L. Tolos and L. Fabbietti, Progress in Particle and Nuclear Physics 112, 103770 (2020).
  * Hebeler [2021] K. Hebeler, Physics Reports 890, 1 (2021), three-nucleon forces: Implementation and applications to atomic nuclei and dense matter.
  * Raaijmakers  _et al._ [2021] G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, A. Schwenk, T. E. Riley, A. L. Watts, J. M. Lattimer, and W. C. G. Ho, The Astrophysical Journal Letters 918, L29 (2021).
  * Li  _et al._ [2021a] J. J. Li, A. Sedrakian, and M. Alford, Phys. Rev. D 104, L121302 (2021a).
  * Legred  _et al._ [2021] I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landry, Phys. Rev. D 104, 063003 (2021).
  * Pang  _et al._ [2021] P. T. H. Pang, I. Tews, M. W. Coughlin, M. Bulla, C. V. D. Broeck, and T. Dietrich, The Astrophysical Journal 922, 14 (2021).
  * Tang  _et al._ [2021] S.-P. Tang, J.-L. Jiang, M.-Z. Han, Y.-Z. Fan, and D.-M. Wei, Phys. Rev. D 104, 063032 (2021).
  * Hebeler  _et al._ [2013a] K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, The Astrophysical Journal 773, 11 (2013a).
  * Weber  _et al._ [2007] F. Weber, R. Negreiros, P. Rosenfield, and M. Stejner, Progress in Particle and Nuclear Physics 59 (2007), 10.1016/j.ppnp.2006.12.008.
  * Alford  _et al._ [2019] M. G. Alford, S. Han, and K. Schwenzer, Journal of Physics G: Nuclear and Particle Physics 46, 114001 (2019).
  * Blaschke [2020] A. A.-C. D. G. H. Blaschke, D.; Ayriyan, UniverseSpecial Issue Compact Stars in the QCD Phase Diagram and in the Multi-Messenger Era of Astronomy 6, 81 (2020).
  * Bao-An Li [2014] A. V. G. V. Bao-An Li, Ramos, The European Physical Journal A 50 (2014), https://doi.org/10.1140/epja/i2014-14009-x.
  * Cai and Li [2021] B.-J. Cai and B.-A. Li, Phys. Rev. C 103, 054611 (2021).
  * Taninah  _et al._ [2020] A. Taninah, S. Agbemava, A. Afanasjev, and P. Ring, Physics Letters B 800, 135065 (2020).
  * Wei  _et al._ [2020] B. Wei, Q. Zhao, Z.-H. Wang, J. Geng, B.-Y. Sun, Y.-F. Niu, and W.-H. Long, Chinese Physics C 44, 074107 (2020).
  * Brockmann and Toki [1992] R. Brockmann and H. Toki, Phys. Rev. Lett. 68, 3408 (1992).
  * Fortin  _et al._ [2016] M. Fortin, C. Providencia, A. R. Raduta, F. Gulminelli, J. L. Zdunik, P. Haensel, and M. Bejger, Phys. Rev. C 94, 035804 (2016), arXiv:1604.01944 [astro-ph.SR] .
  * Fortin  _et al._ [2021] M. Fortin, A. R. Raduta, S. Avancini, and C. Providência, Phys. Rev. D 103, 083004 (2021), arXiv:2102.07565 [nucl-th] .
  * Beznogov and Yakovlev [2015] M. V. Beznogov and D. G. Yakovlev, Mon. Not. Roy. Astron. Soc. 452, 540 (2015), arXiv:1507.04206 [astro-ph.SR] .
  * Providência  _et al._ [2018] C. Providência, M. Fortin, H. Pais, and A. Rabhi,  (2018), 10.3389/fspas.2019.00013, arXiv:1811.00786 [astro-ph.HE] .
  * Lim  _et al._ [2021] Y. Lim, A. Bhattacharya, J. W. Holt, and D. Pati, Phys. Rev. C 104, L032802 (2021).
  * Imam  _et al._ [2022] S. M. A. Imam, N. K. Patra, C. Mondal, T. Malik, and B. K. Agrawal, Phys. Rev. C 105, 015806 (2022).
  * Mondal and Gulminelli [2022] C. Mondal and F. Gulminelli, Phys. Rev. D 105, 083016 (2022).
  * Malik  _et al._ [2022a] T. Malik, M. Ferreira, B. K. Agrawal, and C. ca Providência, The Astrophysical Journal 930, 17 (2022a).
  * Hebeler  _et al._ [2013b] K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013b), arXiv:1303.4662 [astro-ph.SR] .
  * Essick  _et al._ [2021] R. Essick, I. Tews, P. Landry, and A. Schwenk, Phys. Rev. Lett. 127, 192701 (2021).
  * de Tovar  _et al._ [2021] P. B. de Tovar, M. Ferreira, and C. c. Providência, Phys. Rev. D 104, 123036 (2021).
  * Malik  _et al._ [2022b] T. Malik, B. K. Agrawal, and C. m. c. Providência, Phys. Rev. C 106, L042801 (2022b).
  * Scurto  _et al._ [2024] L. Scurto, H. Pais, and F. Gulminelli, Phys. Rev. D 109, 103015 (2024), arXiv:2402.15548 [nucl-th] .
  * Li  _et al._ [2021b] B.-A. Li, B.-J. Cai, W.-J. Xie, and N.-B. Zhang, Universe 7, 182 (2021b), arXiv:2105.04629 [nucl-th] .
  * Oertel  _et al._ [2017b] M. Oertel, M. Hempel, T. Klähn, and S. Typel, Rev. Mod. Phys. 89, 015007 (2017b), arXiv:1610.03361 [astro-ph.HE] .
  * Li and Han [2013] B.-A. Li and X. Han, Phys. Lett. B 727, 276 (2013), arXiv:1304.3368 [nucl-th] .
  * Lattimer and Lim [2013] J. M. Lattimer and Y. Lim, Astrophys. J. 771, 51 (2013), arXiv:1203.4286 [nucl-th] .
  * Typel  _et al._ [2010] S. Typel, G. Röpke, T. Klähn, D. Blaschke, and H. H. Wolter, Phys. Rev. C 81, 015803 (2010).
  * Huang  _et al._ [2020] K. Huang, J. Hu, Y. Zhang, and H. Shen, Astrophys. J. 904, 39 (2020), arXiv:2008.04491 [nucl-th] .
  * Typel and Terrero [2020] S. Typel and D. A. Terrero, The European Physical Journal A 56, 160 (2020).
  * Lalazissis  _et al._ [2005] G. A. Lalazissis, T. Nikisic, D. Vretenar, and P. Ring, Phys. Rev. C 71, 024312 (2005).
  * Antićc and Typel [2015] S. Antićc and S. Typel, Nuclear Physics A 938, 92 (2015).
  * Mishra  _et al._ [2015] R. N. Mishra, H. S. Sahoo, P. K. Panda, N. Barik, and T. Frederico, Phys. Rev. C 92, 045203 (2015).
  * Li  _et al._ [2007] T. Li, U. Garg, Y. Liu, R. Marks, B. Nayak, P. M. Rao, M. Fujiwara, H. Hashimoto, K. Kawase, K. Nakanishi, _et al._ , Physical review letters 99, 162503 (2007).
  * Chen  _et al._ [2005] L.-W. Chen, C. M. Ko, and B.-A. Li, Phys. Rev. Lett. 94, 032701 (2005), arXiv:nucl-th/0407032 .
  * Li  _et al._ [2008] B.-A. Li, L.-W. Chen, and C. M. Ko, Phys. Rept. 464, 113 (2008), arXiv:0804.3580 [nucl-th] .
  * Centelles  _et al._ [2009] M. Centelles, X. Roca-Maza, X. Vinas, and M. Warda, Phys. Rev. Lett. 102, 122502 (2009), arXiv:0806.2886 [nucl-th] .
  * Drischler  _et al._ [2016] C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev. C 93, 054314 (2016), arXiv:1510.06728 [nucl-th] .
  * Brendan T. Reed and Piekarewicz [2023] C. J. H. Brendan T. Reed, F. J. Fattoyev and J. Piekarewicz, (2023), arXiv:2305.19376v2 [Nuclear Theory (nucl-th)] .
  * Adhikari  _et al._ [2021] D. Adhikari _et al._ (PREX), Phys. Rev. Lett. 126, 172502 (2021), arXiv:2102.10767 [nucl-ex] .
  * Adhikari  _et al._ [2022] D. Adhikari _et al._ (CREX), Phys. Rev. Lett. 129, 042501 (2022), arXiv:2205.11593 [nucl-ex] .
  * Margueron  _et al._ [2018] J. Margueron, R. Hoffmann Casali, and F. Gulminelli, Phys. Rev. C 97, 025805 (2018), arXiv:1708.06894 [nucl-th] .
  * Mondal  _et al._ [2017] C. Mondal, B. K. Agrawal, J. N. De, S. K. Samaddar, M. Centelles, and X. Viñas, Phys. Rev. C 96, 021302 (2017), arXiv:1708.03846 [nucl-th] .
  * Providência, Constança  _et al._ [2014] Providência, Constança, Avancini, Sidney S., Cavagnoli, Rafael, Chiacchiera, Silvia, Ducoin, Camille, Grill, Fabrizio, Margueron, Jérôme, Menezes, Débora P., Rabhi, Aziz, and Vidaña, Isaac, Eur. Phys. J. A 50, 44 (2014).
  * Ducoin  _et al._ [2010] C. Ducoin, J. Margueron, and C. Providência, Europhysics Letters 91, 32001 (2010).
  * Ducoin  _et al._ [2011] C. Ducoin, J. Margueron, C. Providência, and I. Vidaña, Phys. Rev. C 83, 045810 (2011).
  * Lattimer  _et al._ [1991] J. M. Lattimer, C. J. Pethick, M. Prakash, and P. Haensel, Phys. Rev. Lett. 66, 2701 (1991).
  * Klähn  _et al._ [2006] T. Klähn, D. Blaschke, S. Typel, E. N. E. van Dalen, A. Faessler, C. Fuchs, T. Gaitanos, H. Grigorian, A. Ho, E. E. Kolomeitsev, M. C. Miller, G. Röpke, J. Trümper, D. N. Voskresensky, F. Weber, and H. H. Wolter, Phys. Rev. C 74, 035802 (2006).
  * Cavagnoli  _et al._ [2011] R. Cavagnoli, D. P. Menezes, and C. m. c. Providência, Phys. Rev. C 84, 065810 (2011).
  * Tolman [1939] R. C. Tolman, Phys. Rev. 55, 364 (1939).
  * Oppenheimer and Volkoff [1939] J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).
  * Baym  _et al._ [1971] G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
  * Malik  _et al._ [2024] T. Malik, H. Pais, and C. Providência, (2024), arXiv:2401.10842 [nucl-th] .
  * Miller  _et al._ [2019] M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, A. K. Harding, W. C. G. Ho, J. M. Lattimer, R. M. Ludlam, S. Mahmoodifar, S. M. Morsink, P. S. Ray, T. E. Strohmayer, K. S. Wood, T. Enoto, R. Foster, T. Okajima, G. Prigozhin, and Y. Soong, The Astrophysical Journal Letters 887, L24 (2019).
  * Riley  _et al._ [2019] T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M. Ludlam, S. Guillot, Z. Arzoumanian, C. L. Baker, A. V. Bilous, D. Chakrabarty, K. C. Gendreau, A. K. Harding, W. C. G. Ho, J. M. Lattimer, S. M. Morsink, and T. E. Strohmayer, The Astrophysical Journal Letters 887, L21 (2019).
  * Carriere  _et al._ [2003] J. Carriere, C. J. Horowitz, and J. Piekarewicz, Astrophys. J. 593, 463 (2003), arXiv:nucl-th/0211015 .
  * Providencia and Rabhi [2013] C. Providencia and A. Rabhi, Phys. Rev. C 87, 055801 (2013), arXiv:1212.5911 [nucl-th] .
  * Riley  _et al._ [2021] T. E. Riley, A. L. Watts, P. S. Ray, S. Bogdanov, S. Guillot, S. M. Morsink, A. V. Bilous, Z. Arzoumanian, D. Choudhury, J. S. Deneva, K. C. Gendreau, A. K. Harding, W. C. G. Ho, J. M. Lattimer, M. Loewenstein, R. M. Ludlam, C. B. Markwardt, T. Okajima, C. Prescod-Weinstein, R. A. Remillard, M. T. Wolff, E. Fonseca, H. T. Cromartie, M. Kerr, T. T. Pennucci, A. Parthasarathy, S. Ransom, I. Stairs, L. Guillemot, and I. Cognard, The Astrophysical Journal Letters 918, L27 (2021).
  * Salinas and Piekarewicz [2023] M. Salinas and J. Piekarewicz, Phys. Rev. C 107, 045802 (2023), arXiv:2301.09692 [nucl-th] .
  * Alam  _et al._ [2016] N. Alam, B. K. Agrawal, M. Fortin, H. Pais, C. Providência, A. R. Raduta, and A. Sulaksono, Phys. Rev. C 94, 052801 (2016), arXiv:1610.06344 [nucl-th] .
  * Char  _et al._ [2023] P. Char, C. Mondal, F. Gulminelli, and M. Oertel, Phys. Rev. D 108, 103045 (2023), arXiv:2307.12364 [nucl-th] .
  * Abbott  _et al._ [2020b] R. Abbott, T. D. Abbott, and et al., The Astrophysical Journal Letters 896, L44 (2020b).


11institutetext: Faculdade de Computação, Universidade Federal do Pará, Belém,
Brazil
11email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Centro de Informática, Universidade Federal
de Pernambuco, Recife, Brazil
22email<EMAIL_ADDRESS>33institutetext: Instituto Tecnológico Vale, Belém,
Brazil
33email<EMAIL_ADDRESS>
# Decoding machine learning benchmarks

Lucas F. F. Cardoso 11 0000-0003-3838-3214 Vitor C. A. Santos 11
0000-0002-7960-3079 Regiane S. Kawasaki Francês 11 0000-0003-3958-064X Ricardo
B. C. Prudêncio 22 0000-0001-7084-1233 Ronnie C. O. Alves 33
0000-0003-4139-0562

###### Abstract

Despite the availability of benchmark machine learning (ML) repositories
(e.g., UCI, OpenML), there is no standard evaluation strategy yet capable of
pointing out which is the best set of datasets to serve as gold standard to
test different ML algorithms. In recent studies, Item Response Theory (IRT)
has emerged as a new approach to elucidate what should be a good ML benchmark.
This work applied IRT to explore the well-known OpenML-CC18 benchmark to
identify how suitable it is on the evaluation of classifiers. Several
classifiers ranging from classical to ensembles ones were evaluated using IRT
models, which could simultaneously estimate dataset difficulty and
classifiers’ ability. The Glicko-2 rating system was applied on the top of IRT
to summarize the innate ability and aptitude of classifiers. It was observed
that not all datasets from OpenML-CC18 are really useful to evaluate
classifiers. Most datasets evaluated in this work (84%) contain easy instances
in general (e.g., around 10% of difficult instances only). Also, 80% of the
instances in half of this benchmark are very discriminating ones, which can be
of great use for pairwise algorithm comparison, but not useful to push
classifiers abilities. This paper presents this new evaluation methodology
based on IRT as well as the tool decodIRT, developed to guide IRT estimation
over ML benchmarks.

###### Keywords:

IRT Machine Learning Benchmarking OpenML Classification.

##  1 Introduction

Machine learning (ML) is a field of artificial intelligence which has been
rapidly growing in recent years, partly due to the diversity of applications
in different areas of knowledge. There are different types of machine learning
algorithms, from unsupervised to supervised ones [3]. In this work, we focus
on supervised algorithms, and more specifically on classification algorithms,
which are commonly adopted in many applications [1].

Empirically evaluating ML algorithms is crucial for assessing the advantages
and limitations of the available techniques. Algorithm evaluation is usually
performed by deploying datasets available in repositories [2]. Remarkably, the
OpenML repository has been developed as an online platform, allowing
researchers to share results, ML strategies and datasets used in benchmarking
experiments. Such platform improves reproducibility and minimizes double
effort to do the same experiment [2]. Another contribution is the OpenML
Curated Classification (OpenML-CC18), a benchmark for classification that has
72 standardized and cured datasets [7].

By relying on benchmark datasets, ML models can be trained and tested using a
chosen experimental methodology (e.g., cross-validation) and the performance
measures of interest (e.g., success rate). For increasing robustness,
different datasets can be adopted in the experiments and the performance
measures are averaged across the datasets. This approach however does not
allow a deeper analysis of classifier’s capacity, since the type of dataset
adopted in an experiment may bias the obtained results. For instance, the
performance of a ML algorithm can be overestimated if only easy datasets are
adopted in the experiments. Also, depending on the datasets, certain ML
algorithms may be favored, thus giving the false impression that one
classifier is the best in relation to the others [4].

Given the above context, it is important that the performance analysis of a
classifier takes into account the complexity of the datasets as well. Previous
work [5, 4, 6] has addressed this issue by adopting concepts from Item
Response Theory (IRT) to provide a more robust approach to devise benchmarks.
IRT is commonly used in psychometric tests to assess the performance of
individuals on a set of items (e.g., questions) with different levels of
difficulty. IRT has been extended to ML evaluation by treating classifiers as
individuals and test instances as items. Thus ML algorithms can be ranked by
IRT according to their ability to correctly respond the most difficult
instances.

The main contributions of this work are: I) an empirical evaluation of IRT to
evaluate ML benchmarks, adopting the OpenML-CC18 benchmark as a case study;
II) a global assessment of ML algorithms in benchmarks comparison, improving
the identification of the strengths and weaknesses among them on the basis of
IRT estimators; and III) the decodIRT tool, developed to guide the IRT
estimation over ML benchmarks automatically along with a robust rating system
to establish proper boundaries between ML algorithms. Regarding the third
contribution, rating systems are widely used in other tasks (e.g., chess
playing) to indicate the strength of an individual in a competition [8]. Among
the existing systems, the Glicko-2 [9] system was used in order to create a
ranking that is able to summarize the algorithms’ results generated by the
IRT. Thus it pointed out the most suitable ML algorithms for the
classification tasks proposed by the OpenML-CC18 benchmark.

The remaining of this paper is organized as follows: Section 2 contextualizes
the main subjects explored in this work, specifically w.r.t. OpenML, Item
Response Theory and the Glicko-2 system. Section 3 presents the methodology
used, explains how decodIRT tool works and how the Glicko-2 system can be
applied to summarize IRT results. Section 4 presents an overall discussion.
Section 5 concludes the article and future developments.

##  2 Background

###  2.1 OpenML

OpenML is a repository in which machine learning researchers can share data
sets, descriptions of experiments and obtained results in as much detail as
possible. Hence, a better organization and use of such information is allowed,
thus creating a collaborative environment for sharing experiments in a global
scale [2]. In OpenML anyone can download a dataset, execute a machine learning
method of their choice and share the results obtained with other users,
generating discussions and new information about the data sets and the applied
algorithms.

In addition, the platform also provides several sets of reference datasets.
OpenML-CC18 111Link to access OpenML-CC18: https://www.openml.org/s/99 is one
of those reference sets that includes 72 existing datasets in OpenML from
mid-2018 and that meet several requirements in order to compile a complete
reference set. OpenML-CC18 also includes datasets frequently used in
benchmarks published in recent years [7]. Several associated metadata is
available such as: number of instances, number of features, number of classes
and proportion between the minority and majority class of the group.

OpenML actually has a lot to contribute to research in the field of machine
learning. The current work aims to apply IRT to evaluate this gold standard to
provide greater strength for OpenML-CC18 and ML community. The analysis of the
datasets within the IRT perspective allows adding to the reference set new
relevant metadata, such as the difficulty of the datasets and discriminatory
potential of the data. Such information can be very useful for choosing a
benchmarking set and will be better discussed in the following sections.

###  2.2 Item Response Theory

Andrade, Tavares and Valle (2000) [10] describe IRT as a set of mathematical
models that aim to represent the probability that an individual will correctly
answer an item according to the item’s parameters and the respondent’s
ability, so that the greater the ability, the greater the chance of a correct
response.

The item can be classified according to its associated response. It can be
dichotomous, when the responses are just right or wrong, or non-dichotomous,
if there are more than two possible answers. Logistical models for dichotomous
items are the most used in literature and practice. There are basically three
types of models, which differ by the number of item parameters that are
estimated. In the 3-parameter logistic model, known as 3PL, the probability of
an individual $j$ correctly respond an item $i$ given his ability is defined
by the following equation:

$P(U_{ij}=1|\theta_{j})=c_{i}+(1-c_{i})\frac{1}{1+e^{-a_{i}(\theta_{j}-b_{i})}}$
(1)

Where:

  * $\bullet$

$U_{ij}$ is the dichotomous response that can take the values 1 or 0, being 1
when the individual j hits the item i and 0 when he misses;

  * $\bullet$

$\theta_{j}$ is the ability of the individual j;

  * $\bullet$

$b_{i}$ is the item’s difficulty parameter and indicates the location of the
logistic curve;

  * $\bullet$

$a_{i}$ is the item’s discrimination parameter, i.e., how much the item i
differentiates between good and bad respondents. This parameter indicates the
slope of the logistic curve. The higher its value, the more discriminating the
item is;

  * $\bullet$

$c_{i}$ is the guessing parameter, representing the probability of a casual
hit. It is the probability that a respondent with low ability hits the item.

The 2PL can be derived by simplifying the above model and dropping the
guessing parameter, i.e., $c_{i}=0$. Finally, the 1PL parameter can be defined
dropping the discrimination parameter, i.e., assuming $a_{i}=1$.

Although it is possible to obtain negative values of discrimination, they are
not expected by IRT. The negative values are not expected because it means
that the probability of a hit is higher for individuals with lower ability
values, rather than higher than what is normally expected [10].

The logistic curves are estimated from responses collected for a group of
items and respondents. According to Martínez-Plumed et al (2016) [4], there
are three possible situations for estimation. In the first situation, the
parameters of the items are known, but the ability of the respondents is not
known. The second possible situation is just the opposite: the ability of each
respondent is defined, but the parameters of the items are not known. In the
third case, which is also the most common, both the parameters of the items
and the abilities of the respondents are unknown. This work is in the third
case and for this situation, the following two-step interactive method
proposed by Birnbaum is used [12] for estimation from the group of responses:

  * •

At first, the parameters of each item are calculated with only the answers of
each individual. The initial values of respondent’s ability can be the number
of correct answers obtained. In the case of classifiers, this work used the
obtained accuracy as the initial ability;

  * •

Once the parameters of the items are obtained, the ability of the individuals
can be estimated. For both item’s parameters and respondent’s ability, simple
estimation techniques can be used, such as the Maximum Likelihood Estimate
[4].

IRT is usually applied for educational purposes, in which respondents are
students and items are questions. Recently, IRT has been extended to AI, and
more specifically to ML, in which respondents are assumed to be AI techniques
and items are AI tasks [5, 4, 6]. To analyze datasets and learning algorithms
through IRT. In this work, it was used the instances of a dataset as the items
and the classifiers were assumed as the respondents. Among the existing IRT
models, the 3PL logistic model was used because it is the most complete. The
IRT parameters were used to assess the datasets directly, informing the
percentage of instances with high difficulty, with great discriminative power
and with high chance of casual hit. This allows to have an insight into the
complexity of the datasets evaluated and how different classifiers behave in
the face of the challenge of classifying different datasets.

###  2.3 Glicko-2 System

The Glicko-2 system is an extension of Glicko, originally developed by Mark E.
Glickman [9] to measure the “strength” of chess players. Today, the system is
used worldwide by several organizations such as the Australian Chess
Federation [11].

In the Glicko-2 system, each player has three variables used to measure his
statistical strength: the R rating, the RD rating deviation and the $\sigma$
volatility. The R rating is the numerical value itself that measures the
strength of a player. As it cannot be said that this value accurately
represents the strength of an individual, the RD is used to indicate the
reliability of an R rating. The Glicko-2 system estimates the players’
strength in a 95% confidence interval, as follows: $[R-2RD,R+2RD]$. Thus,
there is a 95% chance of the real player’s strength is within the calculated
range. The lower a player’s RD, the more accurate his rating value is. In
addition, the $\sigma$ volatility makes it possible to measure the degree of
expected rating fluctuation, i.e., the higher the volatility value, the
greater the chances of a player’s rating fluctuating within his RD range.
Therefore, the lower the value of $\sigma$, the more reliable and expected a
player’s rating is. For example, among players with low volatility, it is
possible to state more precisely who is the strongest based on their ratings
[11]

To estimate the rating values of the players, the Glicko-2 system uses the so-
called classification period. This period is a sequence of games played by the
individual in question. At the end of this sequence, the Glicko-2 system can
update the player’s parameters. To do this, the Glicko-2 system uses the
opponents’ rating and RD values, along with the score of the game result
(e.g., 0 points for defeat and 1 points for victory). If there is no previous
player data, the Glicko system uses default initial values for new players,
1500 for rating, 350 for RD and 0.06 for $\sigma$ [9].

##  3 Materials and methods

In order to build the IRT models and analyze ML benchmarks, it was developed a
tool called decodIRT 222Link to the source code:
https://github.com/LucasFerraroCardoso/IRT_OpenML. The tool was written mostly
in Python together with R language features. The decodIRT has the main
objective of assisting the analysis of existing datasets on the OpenML
platform as well as the proficiency of different classifiers. For this, it
relies on the probability of correctness derived by the IRT model as well as
the estimated items’ parameters and respondents’ ability.

As it can be seen in Figure 1, the decodIRT tool consists of three scripts
designed to be used in sequence. The first script is in charge of downloading
the OpenML datasets and generating several ML models and putting them to
classify the datasets. It will generate a response matrix, which contains the
result of the classification of each test instance. The response matrix will
be the input for the second script, which in turn is in charge of calculating
the item parameters. The last script will use the data generated by the
previous scripts to rank the datasets using the item parameters and use the
IRT to analyze the ML models.

Figure 1: Flowchart of the decodIRT execution.

###  3.1 decodIRT_OtML

The first script will download the existing OpenML datasets and execute
several classifiers implemented in the Scikit-learn [13] package. These
classifiers will generate the response data that will be used to estimate the
item parameters. After downloading the dataset, the stratified split for
training and testing is performed, to maintain the proportions of the classes.
By default, the 70/30 split is used. However, for very large datsets the split
is adapted to leave a maximum of 500 instances for testing. Martínez-Plumed et
al. (2019) [6] explain that for a very large number of items, packages that
estimate IRT values may be stuck at the local minimum or may not converge.
Therefore, for very large datasets it is recommended to use less than a
thousand instances to generate the parameters.

In order to generate diverse response values, at first, 120 models of Neural
Networks (MLP) are executed, varying only the depth of the networks in which
the depth starts at 1 and gradually increases to 120. Each of the MLPs is
trained using the 10-fold cross-validation.

After that, classifications are performed with a second set of classifiers
that will be evaluated by the tool, which are: Naive Bayes Gaussian standard,
Naive Bayes Bernoulli standard, KNN of 2 neighbors, KNN of 3 neighbors, KNN of
5 neighbors, KNN of 8 neighbors , Standard Decision Trees, Random Forests (RF)
with 3 trees, Random Forests with 5 trees, Standard Random Forests, Standard
SVM and Standard MLP. The models defined with "default" means that the default
parameter values that Scikit-learn provide has been used for each ML
algorithm.

Based on the work of Martínez-Plumed et al. (2016) [4], It was inserted 7
artificial classifiers, three random (classify randomly), a majority
classifier (classifies all instances with the predominant class), a minority
classifier (classifies all instances with the non-predominant class), a very
bad one (it misses all classifications) and a great one (it gets all
classifications right). Artificial classifiers serve to set proper performance
indicators regarding ideal classification boundaries.

###  3.2 decodIRT_MLtIRT

In the second script, the IRT item parameters are calculated based on the
classifier responses that were generated in the previous step. For this, it
was necessary to communicate with the R language through the use of Rpy2
library. The R language was used due to the Ltm package, which implements a
framework containing several mechanisms for the calculation and analysis of
IRT for dichotomous data, as well as the generation of item parameters.

As mentioned previously, a maximum of 500 instances were used to generate the
item parameters. This limit was chosen since for larger values it was not
possible to generate the parameters for some datasets.

###  3.3 decodIRT_analysis

The last step is in charge of doing the analysis on the data generated by the
previous steps. By definition, the code will generate in all executions a
ranking of the datasets on the percentage of test instances, whose parameter
value exceeds the defined limit. For example, 50% of instances have a
difficulty value above 1.

If the user does not define any limit for the parameters, standard limit
values based on Adedoyin, Mokobi et al. (2013) [14] are used. In Adedoyin’s
work it is said that to consider an item (instance) as being difficult its
value of Difficulty should normally be higher than 1. Items with high
discriminative capacity should normally have Discrimination values greater
than 0.75, and items with high guessing values are normally greater than 0.2.
The analysis of the percentages of each item parameter is one of the interests
of this work.

Following the Birnbaum method, to calculate the probability of success for
each classifier, the proficiency of each classifier is estimated first. To do
so, Python’s Catsim package is used, which implements both functions to
calculate the probability of success, as well as functions to estimate
proficiency using the responses of each classifier and the item parameters.

In order to analyze more generally how the different classifiers perform on
the datasets. The concept of True-Score [15] was also implemented. The True-
Score is the sum of all the probabilities of a student’s (classifier) correct
answer to the questions of a test. The purpose of using True-Score is to
calculate a final grade, just as we get when we take a test. Thus, it is
expected to identify the "aptitude" that the different classifiers have before
being presented a dataset. A conceptual view of the innate ability,
engineering of the model, through its default hyperparameters settings, it
will be put to the test by the Glicko-2 system.

###  3.4 Ranking of classifiers by the Glicko-2 system

Glicko-2 [9] has been used to estimate the strength of classifiers and to
generate a global performance assessment by simulating a competition among all
ML methods. The simulation explores a round-robin tournament, thus all
classifiers face each other in the championship and at the end there is a
ranking that will reflect how each ML algorithm performed over the
competition.

Each dataset is seen as the championship phases and is used as a
classification period, so the classifiers face each other in each dataset. The
True-Score values obtained by each classifier are used to decide the winner of
the match and this is done as follows: if the score is higher than that of
your opponent, it is counted as victory, if it is lower it is counted as
defeat and if it is the same it is given as a tie. For that, the scoring
system used in chess competitions is used, with 1 point for victory, 0 for
defeat and 0.5 for a tie. Thus, at the end of each dataset, all classifiers
will have their rating, RD and volatility values. These new values are used as
the initial values for the next dataset. Once all datasets are finalized, the
rating values generated are used to create the ranking of classifiers.

##  4 Results and discussion

One of the main objectives of this work is to use the IRT item parameters to
give more reliability to the use of the OpenML-CC18 benchmark. The item
parameters explored are those of difficulty and discrimination. Although there
are 72 datasets available in OpenML-CC18, only 60 were evaluated by decodIRT.
This is due to two reasons: (1) the size of datasets, being 11 too large and
require a long time for the execution of all ML models; (2) It was not
possible to generate the item parameters for the “Pc4” dataset. The latter is
still an open issue. Thus, the OpenML-CC18 benchmark evaluated in this work
corresponds to 83.34 % of the total reference set. These datasets are then
used to performance analysis 333All classification results can be obtained at
https://github.com/LucasFerraroCardoso/IRT_OpenML/tree/master/benchmarking.

###  4.1 Decoding OpenML-CC18 Benchmark

Figure 2: It shows the percentages of difficult and very discriminative
instances arranged in ascending order. There is a certain percentage of
discrimination and a percentage of difficulty that are in the same position on
the X axis do not necessarily correspond to the same dataset. “tic-tac-toe”,
“credit-approval” and “optdigits” are respectively the datsets with the most
difficult instances. While “banknote-authentication”, “analcatdata_authorship”
and “texture” are the most discriminative.

It was possible to observe that there is a inverse relationship between
difficulty and discrimination parameters. The rankings highlights that the
most difficult datasets are also the least discriminative and vice versa
(Figure 2). Thus, the difficulty posed by these datasets, despite make them
more challenging, they are less suitable to differentiate classifiers. On the
other hand, the most discriminatory datasets are not suitable for testing
classifiers abilities.

Among the analyzed datasets, 49 have less than 27% of instances considered
difficult, while only 7 have more than 50% difficult. This means that of the
60 datasets observed, 81.67% of the total have more than 70% of their
instances considered easy and only 11.67% are difficult. Therefore, the
OpenML-CC18 benchmark, for model comparison purposes, should be adopted with
caution. Figure 2 illustrates graphically the classification challenge
presented in this benchmark. It is clear that the situation is reversed when
compared to the difficulty (blue line) values. In which only 25% of datasets
have less than 50% of their instances considered to be less discriminatory and
31 of the 60 datasets have at least 80% of very discriminative instances. From
this observation one can infer that although the OpenML-CC18 is not considered
as difficult as it was expected, it does have appropriate datasets to
differentiate good and bad classifiers. In addition, knowing which datasets
are more difficult, allows the user to choose his benchmark set more
specifically, in case he wants to test the classification power of his
algorithm, disregarding the need to test with all datasets. In summary,
focusing on datasets that present instances with a higher positive
discrimination values.

###  4.2 Classifiers performance on OpenML-CC18

Taking into account only the values of True-Score 444All data generated can be
accessed at
https://github.com/LucasFerraroCardoso/IRT_OpenML/tree/master/BRACIS of the
classifiers for each dataset, it possible to see a pattern in the overall
performance of the classifiers. As expected, the optimal classifier has the
highest score and the other artificial classifiers, the worst results. In
addition, there are cases where the real classifiers match the score of the
optimum and even surpass it. However, in some specific cases the score is
reversed and the best classifiers have the worst scores (see Figure 3).

The score of the optimal classifier be surpassed and the last case mentioned
above can happen due to the occurrence of many instances having negative
values of discrimination. Therefore, it can be inferred that datasets having
many negative instances may not be good for creating benchmarks. A future work
would be to carefully analyze how the characteristics of a dataset can cause
such situations, in addition to observing how these characteristics can
positively or negatively influence the performance of the classifiers. OpenML
already has an extensive set of metadata which can be used for this purpose.

Figure 3: The True-Score values obtained for the “jm1” dataset.

Although the True-Score values generated already allow observing the
performance of the classifiers, there is still a lot of data to be evaluated
individually, towards a natural assessment of the innate ability of these
learning algorithms. Thus, the Glicko-2 rating system was used to summarize
this information. Table 1 shows the final rating ranking that was obtained.

Table 1: Classifier rating ranking Classifier | Rating | RD | Volatility
---|---|---|---
optimal | 1732.56 | 33.25 | 0.0603
MLPClassifier | 1718.65 | 31.20 | 0.0617
RandomForestClassifier | 1626.60 | 30.33 | 0.0606
RandomForestClassifier(5_estimators) | 1606.69 | 30.16 | 0.0621
RandomForestClassifier(3_estimators) | 1575.26 | 30.41 | 0.0646
DecisionTreeClassifier() | 1571.46 | 31.16 | 0.0674
SVM | 1569.48 | 32.76 | 0.0772
KNeighborsClassifier(3) | 1554.15 | 30.74 | 0.0646
GaussianNB | 1530.86 | 31.25 | 0.0683
KNeighborsClassifier(2) | 1528.41 | 30.40 | 0.0638
KNeighborsClassifier(5) | 1526.10 | 30.27 | 0.0630
BernoulliNB | 1494.87 | 32.64 | 0.0770
KNeighborsClassifier(8) | 1457.78 | 30.25 | 0.0638
minority | 1423.01 | 30.66 | 0.0631
rand2 | 1374.78 | 30.27 | 0.0605
rand3 | 1337.27 | 30.95 | 0.0600
rand1 | 1326.38 | 31.42 | 0.0610
majority | 1301.08 | 31.74 | 0.0666
pessimal | 1270.46 | 31.74 | 0.0603

The positions of the artificial classifiers in the ranking were as expected.
The excellent and very bad classifiers assumed the first and last positions,
respectively, while the other artificial classifiers performed less than all
real classifiers. Among the real classifiers, MLP was the winner, being
slightly better than RF classifiers.

When looking at the position of the real classifiers, it is surprising how the
MLP got very close to the rating of the optimal classifier, while the other
classifiers were more than 100 points below the first place. MLP’s “strength”
can also be confirmed by looking at the volatility values of the classifiers,
in which all have low values. Such values also help to give more confidence
with respect to the position of each classifier in the ranking and allow us to
infer that the rating values in fact represent the “strength” of the
classifiers with precision. Considering the rating fluctuation of each
classifier within its respective RD range, the final positions may change. The
MLP, for example, if it had its “strength” at the lowest value by your RD
range, would have a new rating of 1656.25. This would allow classifiers who
are 3rd and 4th to overtake if their rating values fluctuate upwards. However,
from the 5th position, no classifier can reach even if its rating fluctuates
to the maximum of its range.

It can be observed that there are groups of classifiers that have equivalent
“strength”, in which it is not possible to define, with high confidence, which
one of the them present highest aptitude in OpenML-CC18 challenge. This
“strength” is translated into the aptitude of these learning algorithms. Since
tests were performed with several datasets presenting distinct IRT estimators,
and classifiers have same configuration model settings, it is assumed that the
all classification results indeed reflects the learning algorithm aptitude, by
design. It is important to note that optimization can have an effect of fine-
tuning the decision boundaries in the most difficult datasets, and
consequently a better performance. On the other hand, it would not allow the
assessment of innate, engineering skills of these classifiers. A future work
could stress how far one can go in ML optimization, extrapolating the
boundaries between ability and aptitude.

In order to bring more credibility to the generated rating values, the
Friedman [16] test was carried out. The objective is to identify whether the
rating values, in fact, allow to differentiate the “innate ability” of the
classifiers. The Friedman test was performed using the rating values of the
real classifiers. The execution of the Friedman Test resulted in a p-value of
approximately $9.36\text{\times}{10}^{-80}$. Subsequently, the Nemenyi [17]
test was applied to compare and identify which of the distributions differ
from each other. Figure 4 presents a Heatmap of the Nemenyi test.

Figure 4: Heatmap generated from the Nemenyi Test, using only the rating
distributions of the real classifiers.

When analyzing the Heatmap, one can observe that the three classifiers having
the highest ratings have also high p values, so they do not differ from each
other. And even though they are the more performant, they all have high values
for at least one different classifier. The other classifiers also have
performance similarities, even though they are from different families in some
cases. Therefore, it is not evident that there is a clear separation of
classifiers into different groups, since some classifiers that are at the top
of the ranking do not differ from all classifiers that are in lower positions.
This leads us to believe that, although the Friedman Test points out that
there are different ML groups, these differences are not statistically
significant enough to identify the highest skilled one.

##  5 Final considerations

This work explores IRT estimation for the evaluation ML benchmarks. ML
benchmarks are commonly used to explore how far ML algorithms can go while
handling datasets. The OpenML-CC18 is a gold standard. Although the ML
community make broad use of it, such utilization might be taken with caution.
From the 60 datasets evaluated in this work, only 12% of their instances are
considered difficult; 80% of the instances in half of this benchmark are very
discriminating ones which can be of great source for comparisons analysis, but
not useful to push classifiers abilities. The benchmark assessment methodology
is provided and it can be reproduced using the decodIRT tool. Even though
classifier abilities are highlighted by IRT, there were also an issue
regarding the innate ability, aptitude, whether one can set the borders
between the ML algorithm (by design) and training (optimization). The IRT
results were explored by rating systems such as the ones used to evaluate the
strength of chess players to establish the ML winner, and consequently
providing a glimpse towards a aptitude score of ML algorithms.

## References

  * [1] Domingos, Pedro (2012). A few useful things to know about machine learning. Communications of the ACM 55.10, pp 78-87.
  * [2] Vanschoren, J., Van Rijn, J. N., Bischl, B., Torgo, L. (2014). OpenML: networked science in machine learning. ACM SIGKDD Explorations Newsletter, 15(2), 49-60.
  * [3] Monard, M. C., Baranauskas, J. A. (2003). Conceitos sobre aprendizado de máquina. Sistemas inteligentes-Fundamentos e aplicações, 1(1), 32.
  * [4] Martínez-Plumed, F., Prudêncio, R. B., Martínez-Usó, A., Hernández-Orallo, J. (2016). Making sense of item response theory in machine learning. In Proceedings of the Twenty-second European Conference on Artificial Intelligence (pp. 1140-1148). IOS Press.
  * [5] Prudêncio, R. B., Hernández-Orallo, J., Martınez-Usó, A. (2015). Analysis of instance hardness in machine learning using item response theory. In Second International Workshop on Learning over Multiple Contexts in ECML 2015. Porto, Portugal.
  * [6] Martínez-Plumed, F., Prudêncio, R. B., Martínez-Usó, A., Hernández-Orallo, J. (2019). Item response theory in AI: Analysing machine learning classifiers at the instance level. Artificial Intelligence, 271, 18-42.
  * [7] Bischl, B., Casalicchio, G., Feurer, M., Hutter, F., Lang, M., Mantovani, R. G., … Vanschoren, J. (2017). OpenML benchmarking suites and the OpenML100. arXiv preprint arXiv:1708.03731.
  * [8] Samothrakis, S., Perez, D., Lucas, S. M., Rohlfshagen, P. (2014). Predicting dominance rankings for score-based games. IEEE Transactions on Computational Intelligence and AI in Games, 8(1), 1-12.
  * [9] Glickman, M. E. (2012). Example of the Glicko-2 system. Boston University, 1-6.
  * [10] de Andrade, D. F., Tavares, H. R., da Cunha Valle, R. (2000). Teoria da Resposta ao Item: conceitos e aplicações. ABE, Sao Paulo.
  * [11] Veček, N., Mernik, M., Črepinšek, M. (2014). A chess rating system for evolutionary algorithms: A new method for the comparison and ranking of evolutionary algorithms. Information Sciences, 277, 656-679.
  * [12] Birnbaum, A. Lord (1968). "Some latent trait models and their use in inferring an examinee’s ability." Statistical theories of mental test scores.
  * [13] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., … & Vanderplas, J. (2011). Scikit-learn: Machine learning in Python. the Journal of machine Learning research, 12, 2825-2830.
  * [14] Adedoyin, O. O., Mokobi, T. (2013). Using IRT psychometric analysis in examining the quality of junior certificate mathematics multiple choice examination test items. International Journal of Asian Social Science, 3(4), 992-1011.
  * [15] Lord, F. M., Wingersky, M. S. (1984). Comparison of IRT true-score and equipercentile observed-score" equatings". Applied Psychological Measurement, 8(4), 453-461.
  * [16] Pereira, D. G., Afonso, A., Medeiros, F. M. (2015). Overview of Friedman’s test and post-hoc analysis. Communications in Statistics-Simulation and Computation, 44(10), 2636-2653.
  * [17] Nemenyi, P. (1962). Distribution-free multiple comparisons. In Biometrics (Vol. 18, No. 2, p. 263). 1441 I ST, NW, SUITE 700, WASHINGTON, DC 20005-2210: International Biometric Soc.