# Comment on “Distinction of Electron Dispersion in Time-Resolved
Photoemission Spectroscopy”

E. E. Krasovskii  Universidad del Pais Vasco/Euskal Herriko Unibertsitatea,
20080 Donostia/San Sebastián, Basque Country, Spain  Donostia International
Physics Center (DIPC), 20018 Donostia/San Sebastián, Basque Country, Spain
IKERBASQUE, Basque Foundation for Science, 48013 Bilbao, Basque Country, Spain
R.O. Kuzian  Donostia International Physics Center (DIPC), 20018 Donostia/San
Sebastián, Basque Country, Spain  Institute for Problems of Materials Science
NASU, Krzhizhanovskogo 3, 03180 Kiev, Ukraine

###### Abstract

In a recent paper [Phys. Rev. Lett. 125, 043201 (2020)] (Ref. 1) Liao et al.
propose a theory of the interferometric photoemission delay based on the
concepts of the photoelectron phase and photoelectron effective mass. The
present comment discusses the applicability and limitations of the proposed
approach based on an ab initio analysis supported by vast literature. Two
central assumptions of the paper are questioned, namely that the photoelectron
can be characterized by a phase (have a well-defined phase velocity), and that
it can always be ascribed an effective mass. Theories based on these concepts
are concluded to be inapplicable to real solids, which is illustrated by the
example of the system addressed in Ref. 1. That the basic assumptions of the
theory are never fulfilled in nature discredits the underlying idea of the
“time-domain interferometric solid-state energy-momentum-dispersion imaging
method” suggested in Ref. 1. Apart from providing a necessary caution to
experimentalists, the present comment also gives an insight into the
photoelectron wave function and points out problems and pitfalls inherent in
modeling real crystals.

##  .1 Critical notes on the use of the concepts of the effective mass and
propagating phase in photoemission from solids

In Ref. [1], it is proposed to relate the interferometric photoemission delay
to the photoelectron (PE) transport time, based on the PE phase $\xi$ and
effective mass $m_{\text{eff}}$. Generally, the concept of effective mass is
applicable neither to the PE transport in real solids nor to the PE
interaction with light. While formally the effective mass tensor
$m_{\text{eff}}$ can be defined for any bulk Bloch state [2], for the PE it is
undefined in the following cases: (i) When the PE energy lies in a bulk gap,
so the outgoing wave ${\Psi_{\text{\sc pe}}}$ contains no propagating waves.
(ii) When several propagating Bloch states contribute comparably to
${\Psi_{\text{\sc pe}}}$. (iii) When the PE mean free path (MFP) is small
enough so the weight of the evanescent waves (complex Bloch vector
${k_{\perp}}$) in the surface region relevant for photoemision is comparable
to that of the propagating wave(s) (real ${k_{\perp}}$). The ubiquity of the
${\mathbf{k}}_{\parallel}$-projected gaps and the multi-Bloch-wave character
of ${\Psi_{\text{\sc pe}}}$ is well documented for simple [3, 4, 5, 6], noble
and $d$-metals [7, 8], layered dichalcogenides [9, 10], etc. (see illustration
in Fig. 1).

These fundamental and essential features are not allowed for in the basic Eq.
(1) of Ref. [1], which depends on the ratio “between the effective energy-
dependent PE masses”. In the cases (i)–(iii) this equation becomes
inapplicable because either $m_{\text{eff}}$ is undefined or PE is
characterized by more than one $m_{\text{eff}}$. This is characteristic of all
real solids, as illustrated for Mg in Fig. 1(a) by the complex band structure
(CBS) [11] responsible for the PE transport. The role of each Bloch wave in
photoemission is characterized by its contribution to the time-reversed LEED
state $\Phi_{\text{\sc leed}}^{*}$ (low energy electron diffraction) [3]. We
consider the same material (see Ref. 12 for a detailed study of photoemission
from Mg up to photon energies of 320 eV) and parameters as in Ref. [1]:
$\hbar\omega_{\text{\sc ir}}=1.6$ eV, $2q=50$, and $2p$ initial states of 50
eV binding energy. The $2q+1$ final state (31.6 eV) consists of two
propagating waves, while the $2q-1$ state (28.4 eV) comprises only one.
However, for 28.4 eV, down to a depth of at least three atomic layers (15
a.u.) the genuinely evanescent waves 111In the presence of absorbing
potential, all ${k_{\perp}}$ have a nonzero imaginary part; the waves that
have complex wave vectors also with a real potential, $V_{\rm i}=0$, are
referred to as genuinely evanescent waves, see Refs. [4] and [11]. generated
by the surface [14] strongly contribute to ${\Psi_{\text{\sc pe}}}$, as is
evident from the density profile of $\Phi_{\text{\sc leed}}^{*}$ in Fig. 1(c).
The depth of 15 a.u. is larger than the MFP, so the evanescent waves are
important both in the PE transport and in optical transitions forming the
sidebands. Thus, even though there is only one real ${k_{\perp}}$, an
effective mass cannot be ascribed to such wave. Furthermore, the strong
admixture of the evanescent waves makes the notion of a “time delay
$\delta\tau_{a_{s}}$ accumulated by PEs traveling over one lattice constant”
inapplicable and Eqs. (7) and (8) of Ref. [1] meaningless.

Figure 1: (a) Black lines: real band structure of Mg along $\Gamma A\Gamma$.
Red: conducting CBS given by the dispersion ${\rm Re}\,{k_{\perp}}(E)$ of
$\Phi_{\text{\sc leed}}^{*}$ constituents. Line thickness is proportional to
the flux carried by the respective Bloch wave [9]. LEED is calculated for
absorbing potential $V_{\rm i}=0.2$ eV [7]. Arrows show the $2q\pm 1$ final
states. (b), (c) Density profile $\rho(z)$ of $\Phi_{\text{\sc leed}}^{*}$
(with $V_{\rm i}=0$) for $2q+1$ (b) and $2q-1$ (c): full density (black) and
single-plane-wave ${\mathbf{G}}_{\parallel}=0$ contribution (red). Gray area
shows the contribution from the surface harmonics
${\mathbf{G}}_{\parallel}\neq 0$.

Let us now imagine a situation when neither of the obstacles (i)–(iii) occurs,
and the outgoing wave in the crystal is a single Bloch wave:

${\Psi_{\text{\sc
pe}}}({\mathbf{r}})=\sum\limits_{{\mathbf{G}}_{\parallel}}\sum\limits_{g}\psi_{{\mathbf{G}}_{\parallel}g}\exp\left[i({k_{\perp}}+g)z+i{\mathbf{G}}_{\parallel}{\mathbf{r}}_{\parallel}\right].$
(1)

Here $g$ are surface-normal and ${\mathbf{G}}_{\parallel}$ are surface-
parallel reciprocal lattice vectors. The theory of Ref. [1] assumes that a
propagation phase $\xi=k\Delta z$ can be introduced for the PE wave [Eq. (2)
of Ref. [1]]. This assumption is justified only if one of the plane waves
$\\{g,{\mathbf{G}}_{\parallel}\\}$ dominates the expansion (1). In reality,
this is almost never the case, see, e.g., the density profiles in Figs. 1(b)
and 1(c): The weight of ${\mathbf{G}}_{\parallel}\neq 0$ waves is of the same
order as (and often larger than) the weight of the
${\mathbf{G}}_{\parallel}=0$ component. Thus, because each of the
${\mathbf{G}}_{\parallel}$-waves acquires its own phase $\xi$ the propagation
phase is ill-defined for a Bloch wave.

The conclusion of Ref. [1] that “interferometric spectroscopy addresses the
material-dependent change of PE phase velocities inside the solids”
contradicts the well-documented multi-plane-wave nature of the Bloch states,
and the notion of effective mass conflicts both with the multi-Bloch-wave
structure of the PE and with the surface sensitivity of photoemission. The
concepts of PE phase and effective mass, which are plausible for one-
dimensional models, turn out to be inapplicable to real solids, which makes
the method to extract photoemission delays based on Eq. (1) of Ref. [1]
unreliable.

##  .2 Addendum: Critique of the Authors’ Response

A shorter version of this note has been published as a Comment [15]. For the
benefit of the reader, we quote the main points of the Authors’ Response [16]
to the Comment.

  * ∙

Ref. 16: “… case (i) is irrelevant since the absence of “propagating waves” in
direction normal to the surface prohibits electron emission.”

In fact, electron emission in the absence of propagating waves is not
prohibited. This is known since early days of photoelectron spectroscopy [17]
and is referred to as “band gap emission”. It is well-documented
experimentally [18] and is discussed in classical textbooks [19, 20]: “it is a
matter of common experience that EDCs in the UPS regime are also observed when
the final state is in the gap.”—Hüfner 2003 [20], p. 383.

  * ∙

Ref. 16: “With regard to case (ii), we point out that outgoing final waves are
specified by their momentum (vector), not energy. Should several propagating
Block states with the same energy contribute, they must have different
momenta, both inside [as shown in Fig. 1(a) of the Comment at 31.6 eV] and
outside the solid, and thus do not contribute to the same outgoing wave. Such
superpositions can be disentangled in time- and angle-resolved photoemission
spectra.”

This statement is incorrect. No experiment can distinguish between
photoelectrons excited into different Bloch waves at the same energy and same
${\mathbf{k}}_{\parallel}$. This is especially obvious in the one-step theory
of photoemission, in which the final state is represented by the time-reversed
LEED state [17, 21, 22, 3], so the photoemission matrix element is a coherent
sum of the contributions from different Bloch waves, see a typical example in
Fig. 2 below. This gives rise to an interference between transitions to
different Bloch waves, and the respective interference terms have been
thoroughly studied in Ref. [9]. The presence of the interference terms makes
it principally impossible to “disentangle” the individual contributions in
angle-resolved photoemission spectra.

Figure 2:  Transitions from the $k^{+}$ and $k^{-}$ constituents of the
initial Te 5$p_{z}$ band of TiTe2 to Bloch-wave constituents of the unoccupied
band structure. Full curves display real band structure and dashed curves real
lines of complex band structure (evanescent waves). Thick lines show
conducting constituents of the final states. The line thickness is
proportional to the current carried (absorbed) by this partial wave. The
transitions indicated by the red straight lines give rise to interference
terms in the matrix element, see a detailed discussion in Ref. [9].

Apart from the general criticism above, we have to point at a serious error in
the response by Liao et al.: They claim that different propagating Bloch
states must have different momenta also outside the solid. Obviously, the wave
that is defined both inside and outside the solid is not a Bloch wave, and the
momentum outside the solid is uniquely determined by the kinetic energy and
the emission direction.

  * ∙

Ref. 16: “In case (iii), the damped amplitude of the evanescent wave does not
directly influence the phase accumulation during propagation.”

This is not an answer to the case (iii) remark. Case (iii) is illustrated in
Fig. 1(c), and it concerns the issue of effective mass rather than phase
accumulation. Case (iii) is a generalization of case (ii), the difference
being that in case (ii) there are several propagating waves, and in case (iii)
there may be one propagating and many evanescent waves. For photoemission,
propagating and evanescent waves are equally important, so there is no
physical difference between the two cases. However, while case (ii) can be
excluded at some energies, case (iii) should be taken into account at every
energy, see the detailed discussions in Refs. [14, 23, 24].

  * ∙

Ref. 16: “… photoelectron transport in the solid target while, to our
knowledge, photoelectron transport times have not been clearly related to the
directly observable phase (difference) information. Equation (1) in Ref. [1]
provides this missing link based on the decisively different effect of
photoelectron energy-momentum dispersion on the active electron’s phase and
group velocity. ”

It is exactly the validity of the concept of wave phase in application to
photoelectron transport that we question in our Comment, in particular, the
relevance of Eq. (1) in Ref. [1]. The authors made no attempt to respond to
our criticism that the propagation phase is ill-defined for a Bloch wave.

  * ∙

Ref. 16: “Clearly, the degree to which predictions from one-dimensional models
and Eq. (1) in Ref. [1] are applicable to real (strong-field dressed three-
dimensional solids) remains to be seen.”

It is not clear how the theory of Ref. [1] can make predictions for specific
materials, but this is not our point. Whatever remains to be seen, ignoring
the limitations of the model and using the concepts of effective mass and wave
phase in cases when they are undefined are essential flaws of the theory.

  * ∙

Ref. 16: “… one-dimensional numerical models have been applied with some
success to model time-resolved photoemission from surfaces in the past”

One-dimensional models may be useful and instructive, but they cannot be used
beyond their limitations.

  * ∙

Ref. 16: “… we (i) disagree with the conclusion of Krasovskii and Kuzian (KK)
that the concept of photoelectron phase and effective mass … are inapplicable
to real solids and (ii) find the statement by KK that for the “cases (i)—(iii)
Eq. (1) becomes inapplicable” overstated and inappropriate.”

It is important to reveal the limitations and delineate the applicability of
the model. Here, we have questioned the main claim of Ref. [1], namely that
“attosecond time-resolved interferometric spectroscopy addresses the material-
dependent change of PE phase velocities inside the solids”, as well as the
main message for experimentalists that the suggested method can be used “to
extract accurately photoemission-time delays for imaging electronic dispersion
in real solids”.

##  .3 Acknowledgments:

This work was supported by the Spanish Ministry of Science and Innovation
(Project No. PID2019-105488GB-I00) and by the National Academy of Sciences of
Ukraine (Project No. III-4-19).

## References

  * Liao  _et al._ [2020] Q. Liao, W. Cao, Q. Zhang, K. Liu, F. Wang, P. Lu, and U. Thumm, Distinction of electron dispersion in time-resolved photoemission spectroscopy, Phys. Rev. Lett. 125, 043201 (2020).
  * Slater [1949] J. C. Slater, Electrons in perturbed periodic lattices, Phys. Rev. 76, 1592 (1949).
  * Kuzian and Krasovskii [2020] R. O. Kuzian and E. E. Krasovskii, One-step theory of photoelectron escape time: Attosecond spectroscopy of mg(0001), Phys. Rev. B 102, 115116 (2020).
  * Krasovskii  _et al._ [2008] E. E. Krasovskii, W. Schattke, P. Jiříček, M. Vondráček, O. V. Krasovska, V. N. Antonov, A. P. Shpak, and I. Bartoš, Photoemission from Al(100) and (111): Experiment and ab initio theory, Phys. Rev. B 78, 165406 (2008).
  * Krasovskii  _et al._ [2010] E. E. Krasovskii, V. M. Silkin, V. U. Nazarov, P. M. Echenique, and E. V. Chulkov, Dielectric screening and band-structure effects in low-energy photoemission, Phys. Rev. B 82, 125102 (2010).
  * Strocov [2018] V. Strocov, Photoemission response of 2D electron states, Journal of Electron Spectroscopy and Related Phenomena 229, 100 (2018).
  * Krasovskii and Schattke [1999] E. E. Krasovskii and W. Schattke, Calculation of the wave functions for semi-infinite crystals with linear methods of band theory, Phys. Rev. B 59, R15609 (1999).
  * Lobo-Checa  _et al._ [2011] J. Lobo-Checa, J. E. Ortega, A. Mascaraque, E. G. Michel, and E. E. Krasovskii, Effect of photoelectron mean free path on the photoemission cross-section of Cu(111) and Ag(111) Shockley states, Phys. Rev. B 84, 245419 (2011).
  * Krasovskii  _et al._ [2007a] E. E. Krasovskii, V. N. Strocov, N. Barrett, H. Berger, W. Schattke, and R. Claessen, Band mapping in the one-step photoemission theory: Multi-Bloch-wave structure of final states and interference effects, Phys. Rev. B 75, 045432 (2007a).
  * Krasovskii  _et al._ [2007b] E. E. Krasovskii, K. Rossnagel, A. Fedorov, W. Schattke, and L. Kipp, Determination of the hole lifetime from photoemission: Ti 3d states in TiTe2, Phys. Rev. Lett. 98, 217604 (2007b).
  * Heine [1963] V. Heine, On the general theory of surface states and scattering of electrons in solids, Proceedings of the Physical Society 81, 300 (1963).
  * Krasovskii [2020] E. E. Krasovskii, Character of the outgoing wave in soft x-ray photoemission, Phys. Rev. B 102, 245139 (2020).
  * Note [1] In the presence of absorbing potential, all ${k_{\perp}}$ have a nonzero imaginary part; the waves that have complex wave vectors also with a real potential, $V_{\rm i}=0$, are referred to as genuinely evanescent waves, see Refs. [4] and [11].
  * Appelbaum and Hamann [1972] J. A. Appelbaum and D. R. Hamann, Self-consistent electronic structure of solid surfaces, Phys. Rev. B 6, 2166 (1972).
  * Krasovskii and Kuzian [2021] E. E. Krasovskii and R. O. Kuzian, Comment on “distinction of electron dispersion in time-resolved photoemission spectroscopy”, Phys. Rev. Lett. 126, 109303 (2021).
  * Liao  _et al._ [2021] Q. Liao, W. Cao, Q. Zhang, K. Liu, F. Wang, P. Lu, and U. Thumm, Liao et al. reply:, Phys. Rev. Lett. 126, 109304 (2021).
  * Feibelman and Eastman [1974] P. J. Feibelman and D. E. Eastman, Photoemission spectroscopy–correspondence between quantum theory and experimental phenomenology, Phys. Rev. B 10, 4932 (1974).
  * Courths  _et al._ [1979] R. Courths, S. Hüfner, and H. Schulz, High-resolution normal emission UV-photoelectron spectra and band structure of Cu, Z. Physik B Condensed Matter 35, 107 (1979).
  * Feuerbacher and Fitton [1977] B. Feuerbacher and B. Fitton, Photoemission spectroscopy, in  _Electron Spectroscopy for Surface Analysis_, edited by H. Ibach (Springer Berlin Heidelberg, Berlin, Heidelberg, 1977) pp. 151–203.
  * Hüfner [2003] S. Hüfner,  _Photoelectron Spectroscopy Principles and Applications_ , Advanced Texts in Physics (Springer-Verlag Berlin Heidelberg, 2003).
  * Freericks  _et al._ [2009] J. K. Freericks, H. R. Krishnamurthy, and T. Pruschke, Theoretical description of time-resolved photoemission spectroscopy: Application to pump-probe experiments, Phys. Rev. Lett. 102, 136401 (2009).
  * Braun  _et al._ [2015] J. Braun, R. Rausch, M. Potthoff, J. Minár, and H. Ebert, One-step theory of pump-probe photoemission, Phys. Rev. B 91, 035119 (2015).
  * Bross [1982] H. Bross, The scattering of conduction electrons in Cu by a stacking fault, Journal of Physics F: Metal Physics 12, 2883 (1982).
  * Wachutka [1986] G. Wachutka, New layer method for the investigation of the electronic properties of two-dimensional periodic spatial structures: First applications to copper and aluminum, Phys. Rev. B 34, 8512 (1986).


# Diquarks and the production of charmed baryons

Hyeongock Yun<EMAIL_ADDRESS>Department of Physics and Institute of
Physics and Applied Physics, Yonsei University, Seoul 03722, Korea  Sungsik
Noh<EMAIL_ADDRESS>Department of Physics and Institute of Physics
and Applied Physics, Yonsei University, Seoul 03722, Korea  Sanghoon Lim
<EMAIL_ADDRESS>Department of Physics, Pusan National University, Pusan,
Korea  Taesoo Song<EMAIL_ADDRESS>GSI Helmholtzzentrum für
Schwerionenforschung GmbH, Planckstrasse 1, 64291 Darmstadt, Germany  Juhee
Hong<EMAIL_ADDRESS>Department of Physics and Institute of Physics and
Applied Physics, Yonsei University, Seoul 03722, Korea  Aaron Park
<EMAIL_ADDRESS>Department of Physics and Institute of Physics and
Applied Physics, Yonsei University, Seoul 03722, Korea  Su Houng Lee
<EMAIL_ADDRESS>Department of Physics and Institute of Physics and
Applied Physics, Yonsei University, Seoul 03722, Korea  Benjamin Dönigus
<EMAIL_ADDRESS>Institut für Kernphysik, Johann Wolfgang Goethe-
Universität Frankfurt, Max-von-Laue-Str. 1, 60438 Frankfurt, Germany

###### Abstract

Utilizing a quark model characterized by parameters that effectively replicate
the masses of ground state hadrons, we illustrate that $(us)$ or $(ds)$
diquarks exhibit greater compactness in comparison to $(ud)$ diquarks.
Concretely, the binding energy of the $(us)$ diquark - defined as the
diquark’s mass minus the combined masses of its individual quarks - is found
to be stronger than that of the $(ud)$ diquark. This heightened attraction
present in $(us)$ diquarks could lead to enhanced production of $\Xi_{c}/D$
particles in high-energy pp or ultrarelativistic heavy-ion collisions.

##  I Introduction

While the main objective of the heavy-ion collision experiments at CERN is to
study the properties of quark-gluon plasma [1], recent results including those
from pp and pPb collisions have provided new opportunities to study QCD and
hadron properties in general [2]. Of particular interest are the recent
measurements of heavy baryon to heavy meson ratios in pp, pPb and PbPb
collisions because these probe the quark-quark interaction in the
hadronization processes [3, 4, 5].

The baryon-to-meson ratio enhancement for both the light and heavy quark
sectors observed in heavy-ion collision can be explained by medium effects in
the hadronization process well encoded in the coalescence model [6, 7, 8, 9].
What is interesting to see is that $\Lambda_{c}/D$ enhancement over that
expected from the fragmentation process exists even in pp collision, which can
only be explained when one assumes an additional production mechanism for the
baryon[10]. In the thermal model additional resonances not listed by the
Particle Data Group [11] are needed [12, 2]. Even more striking, the models
that explain the $\Lambda_{c}/D$ enhancement underestimate the recent
measurement of the $\Xi_{c}/D$ ratio[13].

In this work, we will show that a strong $(qs)$ diquark correlation, where
$q=u,d$, provides a new effect that provides an additional production
mechanism to enhance the $\Xi_{c}/D$ ratio.

The two light quarks in the ground state baryon with one heavy quark are in
the color anti-triplet channel with either isospin zero (spin zero) or isospin
one (spin one), which respectively form the lowest triplet or sextet charmed
baryon representations in flavor SU(3), respectively. Hence, if strong diquark
correlations exist, a strong overlap of the diquarks and a heavy quark could
lead to additional enhancement in the heavy baryon-to-meson ratio. On the
other hand, data seems to suggest that there should be an additional
enhancement in the $\Xi_{c}/D$ ratio compared to $\Lambda_{c}/D$ ratio in pp
collision at both $\sqrt{s}=5.02$ TeV and 13 TeV [4, 5]. The diquark in the
$\Xi_{c}$ and $\Lambda_{c}$ are $(qs)$ and $(ud)$ diquarks, respectively. As
we will see, the quark model that well reproduces the ground state hadron
masses shows that there is a stronger attraction in the $(qs)$ diquark
compared to that in the $(ud)$ diquark. This will lead to a small but non-
negligible enhancement in $\Xi_{c}$ productions. Enhanced production in the
presence of attraction can be shown to be true in general using S-matrix
theory for hadron production [14]. Here, we will employ a coalescence model to
estimate the extent of the additional enhancement anticipated in the
$\Xi_{c}/D$ ratio as a result of the strong $(qs)$-diquark correlation, which
is represented by a corresponding diquark distribution inherent in the
hadronization process

The paper is organized as follows. In Sec. II, we discuss the strength and the
binding energies of the $(us)$ and $(ud)$ diquarks. We then construct a
phenomenological coalescence model to estimate the additional $\Xi_{c}/D$
ratio production in pp, pPb, and PbPb collisions. Finally, we give the
summary.

##  II Diquarks in quark model

It has been long noticed that diquarks play important roles in hadron
structures and reactions [15]. To study the flavor dependence on their
binding, let us study the diquark configurations in a quark model.
Specifically, we evaluate the masses of diquarks by using a non-relativistic
quark model described by the following Hamiltonian[16].

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum^{n}_{i=1}\left(m_{i}+\frac{{\mathbf{p}}^{2}_{i}}{2m_{i}}\right)-\frac{3}{4}\sum^{n}_{i<j}\frac{\lambda^{c}_{i}}{2}\,\,\frac{\lambda^{c}_{j}}{2}\left(V^{C}_{ij}+V^{CS}_{ij}\right),\qquad$
(1)

where $n=2(3)$ for diquarks (baryons). The internal quark potentials
$V^{C}_{ij}$ and $V^{CS}_{ij}$ in Eq. (1) are as follows.

$\displaystyle V^{C}_{ij}$ $\displaystyle=$
$\displaystyle-\frac{\kappa}{r_{ij}}+\frac{r_{ij}}{a^{2}_{0}}-D,$ (2)
$\displaystyle V^{CS}_{ij}$ $\displaystyle=$
$\displaystyle\frac{\hbar^{2}c^{2}\kappa^{\prime}}{m_{i}m_{j}c^{4}}\frac{e^{-\left(r_{ij}\right)^{2}/\left(r_{0ij}\right)^{2}}}{(r_{0ij})r_{ij}}\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j}\,,$
(3)

where

$\displaystyle r_{0ij}$ $\displaystyle=$ $\displaystyle
1/\left(\alpha+\beta\frac{m_{i}m_{j}}{m_{i}+m_{j}}\right)\,,$
$\displaystyle\kappa^{\prime}$ $\displaystyle=$
$\displaystyle\kappa_{0}\left(1+\gamma\frac{m_{i}m_{j}}{m_{i}+m_{j}}\right)\,.$

This model can be used to fit the ground state hadron masses including light,
charm, and bottom quarks[16, 17]. The model has also been extensively used to
study possible compact exotic configurations[16, 17, 18]. Here, we fit the
model parameters in the Hamiltonian to the ground state masses listed in Table
LABEL:baryons. The list is limited to baryon states to better fit the quark-
quark interactions with different flavors. The selected baryons in Table
LABEL:baryons contain the diquark structure of interest in this work. The
fitted model parameters are as follows.

$\displaystyle\kappa=130.0\,\textrm{MeV fm},\quad
a_{0}=0.0318119\,\textrm{(MeV${}^{-1}$fm)${}^{1/2}$},$ $\displaystyle
D=975\,\textrm{MeV},\quad m_{u}=m_{d}=335\,\textrm{MeV},$ $\displaystyle
m_{s}=652\,\textrm{MeV},\quad m_{c}=1940\,\textrm{MeV},$
$\displaystyle\alpha=1.2499\,\textrm{fm${}^{-1}$},\,\,\beta=0.0008314\,\textrm{(MeV
fm)${}^{-1}$},$
$\displaystyle\gamma=0.00168\,\textrm{MeV${}^{-1}$},\,\,\kappa_{0}=185.144\,\textrm{MeV}.$
(4)

Table 1: Masses of baryons obtained (Column 3) from the model calculation in this work with the fitting parameter set given in Eq. (4). Column 4 shows the variational parameters $a_{1}$ and $a_{2}$. Particle | Experimental | Mass |  Variational
---|---|---|---
Value (MeV) | (MeV) |  Parameters (${\rm fm}^{-2}$)
$\Lambda$ | 1115.7 | 1116.3 |  $a_{1}$ = 3.0, $a_{2}$ = 2.9
$\Lambda_{c}$ | 2286.5 | 2272.2 |  $a_{1}$ = 3.1, $a_{2}$ = 3.9
$\Sigma_{c}$ | 2452.9 | 2446.2 |  $a_{1}$ = 2.2, $a_{2}$ = 4.0
$\Sigma_{c}^{*}$ | 2517.5 | 2531.4 |  $a_{1}$ = 2.0, $a_{2}$ = 3.5
$\Sigma$ | 1192.6 | 1202.2 |  $a_{1}$ = 2.2, $a_{2}$ = 3.3
$\Sigma^{*}$ | 1383.7 | 1401.1 |  $a_{1}$ = 1.9, $a_{2}$ = 2.5
$\Xi$ | 1314.9 | 1331.8 |  $a_{1}$ = 3.6, $a_{2}$ = 3.1
$\Xi^{*}$ | 1531.8 | 1544.1 |  $a_{1}$ = 3.1, $a_{2}$ = 2.3
$\Xi_{c}$ | 2467.8 | 2474.2 |  $a_{1}$ = 3.5, $a_{2}$ = 5.0
$\Xi_{c}^{*}$ | 2645.9 | 2654.9 |  $a_{1}$ = 2.6, $a_{2}$ = 4.6
$\Xi^{\prime}_{c}$ | 2579.2 | 2570.2 |  $a_{1}$ = 2.8, $a_{2}$ = 5.2
$\Omega_{c}$ | 2695.2 | 2684.7 |  $a_{1}$ = 3.9, $a_{2}$ = 6.1
$\Omega_{c}^{*}$ | 2765.9 | 2768.8 |  $a_{1}$ = 3.6, $a_{2}$ = 5.3
$p$ | 938.27 | 951.42 |  $a_{1}$ = 2.5, $a_{2}$ = 2.5
$\Delta$ | 1232 | 1246.9 |  $a_{1}$ = 1.8, $a_{2}$ = 1.8

The two variational parameters $a_{1}$ and $a_{2}$ appearing in Table
LABEL:baryons are the scaling factors of the Gaussian wave function for the
relative distance between the two light quarks within the diquark and the
relative distance between the center of the diquark and the heavy quark,
respectively. The standard deviation of the masses obtained in Table
LABEL:baryons is
$\sigma=(\frac{1}{N-1}\sum_{i=1}^{N}\left(M^{Thr}_{i}-M^{Exp}_{i}\right)^{2})^{1/2}=11.9~{}{\rm
MeV},$ where $M^{Thr}_{i}$ indicates the mass obtained by the model
calculation and $M^{Exp}_{i}$ the experimentally measured mass.

Table 2: Mass, binding energy, and size of $ud$, $us$, and $uc$ diquarks in a color $\mathbf{\bar{3}}$ and spin 0 state. The binding energy is defined as $M_{Diquark}-m_{u}-m_{q}$ where $m_{q}$ can be $m_{d}$, $m_{s}$, or $m_{c}$. Diquark |  | Mass(MeV) |  | Binding energy(MeV) |  | Size(fm)
---|---|---|---|---|---|---
$ud$ |  | 680.3 |  | 10.22 |  | 0.761
$us$ |  | 970.7 |  | -16.30 |  | 0.714
$uc$ |  | 2244.8 |  | -30.23 |  | 0.700
Table 3: Individual contribution to the mass of $ud$, $us$, and $uc$ diquarks in Table 2. In the following table, $m_{q}$ can be $m_{d}$, $m_{s}$, or $m_{c}$. All units are MeV. Contribution |  | $ud$ |  | $us$ |  | $uc$
---|---|---|---|---|---|---
$m_{u}+m_{q}-\frac{1}{2}D$ |  | 182.5 |  | 499.5 |  | 1787.5
$\sum^{2}_{i=1}\frac{{\mathbf{p}}^{2}_{i}}{2m_{i}}$ |  | 383.6 |  | 329.9 |  | 265.8
$V^{C}$ |  | 267.1 |  | 236.6 |  | 227.5
$V^{CS}$ |  | -152.9 |  | -95.34 |  | -36.0
Total |  | 680.3 |  | 970.7 |  | 2244.8

We now calculate the masses of the $(ud)$ and $(us)$ diquarks $M_{Diquark}$
using these parameters. The results are summarized in Table 2. In the table,
one can see that the binding is stronger for the heavier diquarks. The binding
energy is defined as $M_{Diquark}-m_{u}-m_{i}$ where $m_{i}$ can be $m_{d}$,
$m_{s}$, or $m_{c}$. Also, one can see that the size of the diquark becomes
smaller when the component quark becomes heavier. To analyze the origin of the
stronger binding for the heavier diquarks, we show each part of the
Hamiltonian contributing to the masses of diquarks in Table 3. We first note
that the attraction coming from the color-spin interaction $V^{CS}$ becomes
smaller as the quark masses increase. This is due to the inverse quark mass
dependence in $V^{CS}$. On the other hand, both the Coulomb and confining
potential $V^{C}$ decrease as the quark masses increase. This is attributed to
the reduced size of heavier diquarks. Heavier quarks experience a more
pronounced Coulomb attraction and a weaker linearly rising potential, a
phenomenon commonly observed when dealing with heavier quarks, as discussed
previously in Ref. [19, 20]. Altogether, one finds that the diquarks become
more bound as the quarks involved become heavier. As a result, in Table 2,
$(us)$ and $(uc)$ diquarks are bound while $(ud)$ diquark is not.

###  II.1 Diquark mass at finite temperature

In reality, one has to introduce a model scenario to implement the strong
$(us)$ correlation at the hadronization point in the presence of quarks. One
potential scenario involves assuming the presence of local quark matter at the
point of hadronization, investigating the persistence of diquark correlations,
and subsequently assessing the impact of diquarks on the process of
hadronization.

For that purpose, we will first analyze the diquarks when the Coulomb and the
confining part of the potential are modified as given in Ref.[21] at the
chemical freeze-out temperature. Then the color potential in Eq. (2) will be
modified as follows.

$\displaystyle V^{C}_{ij}$ $\displaystyle=$
$\displaystyle-\kappa\bigg{[}m_{D}+\frac{e^{-m_{D}r_{ij}}}{r_{ij}}\bigg{]}$
(5)
$\displaystyle+\frac{1}{a^{2}_{0}}\bigg{[}\frac{2}{m_{D}}-\frac{e^{-m_{D}r_{ij}}(2+m_{D}r_{ij})}{m_{D}}\bigg{]}-D,$

which reduces to Eq. (2) when the Debye screening mass $m_{D}\rightarrow 0$.
Here we will consider $m_{D}$ given in Ref. [21] at several temperatures above
the critical point.

Table 4: Differences in binding energies $B_{(us)}-B_{(ud)}$ where $B_{(us)}$ and $B_{(ud)}$ indicate the binding energies of $(us)$ and $(ud)$ diquarks, respectively. The mass and size of $(us)$ diquark are also presented. The table presents the results obtained from five different schemes discussed in the text. Type | $B_{(us)}-B_{(ud)}$ |  $m_{us}$(MeV) | Size(fm)
---|---|---|---
Scheme 1 | -26.52 | 970.7 | 0.714
Scheme 2 | -11.14 | 710.5 | 0.753
Scheme 3 | -9.67 | 673.2 | 0.800
Scheme 4 | -7.81 | 613.8 | 0.889
Scheme 5 | -6.93 | 580.2 | 0.954

For the light-quark system, we will further consider the effects of chiral
symmetry restoration and thermal masses. To this end, in evaluating the masses
of diquarks within our quark model approach, we will study cases where we use
the thermal masses for the light quarks, typically taken to be 300 MeV for the
$u,d$ quarks, and 400 MeV for the strange quark, which is about 100 MeV
heavier than the light quark mass as is the case for the bare quark masses.

We now solve for the diquark masses and their bindings for the following five
cases. 1) Scheme 1(Table 2): use the potential and quark masses as fitted to
the baryon mass spectrum. That is, use zero temperature potential in Eq. (2)
with $m_{s}=652$ MeV and $m_{u}=335$ MeV. 2) Scheme 2: use the zero
temperature potential in Eq. (2) but with $m_{s}=400$ MeV and $m_{u}=300$ MeV.
3) Scheme 3: use the thermal potential in Eq. (5) with $m_{s}=400$ MeV,
$m_{u}=300$ MeV, and $m_{D}=T_{c}=156$ MeV. 4) Scheme 4: use the thermal
potential in Eq. (5) with $m_{s}=400$ MeV, $m_{u}=300$ MeV, and $m_{D}\sim
1.5\times T_{H}$ MeV[21], where we take $T_{H}=181$ MeV to be the
hadronization temperature in pp collision from the flow analysis given later.

We also introduce Scheme 5, aimed at describing the phase where chiral
symmetry is restored. From a phenomenological point of view, the value of the
potential at $r\rightarrow\infty$ is related to the creation of the quark-
antiquark pair, whose mass is related to chiral symmetry breaking[22]. When
$m_{D}$ is small, the potential is very large due to the linearly rising
potential. That is the point where the usual quark model is applied to
calculate the ground state hadron masses; the potential at large values does
not affect the properties of the ground state hadrons as they are smaller than
1 fm in size. On the other hand, $m_{D}$ increases with temperature and when
$m_{D}\sim 0.325$ GeV, $V(r\rightarrow\infty)=0$. Lattice gauge theory fit to
$m_{D}$ indeed indicates that 0.325 GeV is reached near the critical
temperature within the lattice error[21]. Therefore, we consider Scheme 5,
where we take $m_{D}=325$ MeV and take the quark masses to be their thermal
masses $m_{s}=400$ MeV and $m_{u}=300$ MeV.

Table 4 shows the result for the differences in the binding energy between the
$(us)$ and $(ud)$ diquarks, the mass of the $(us)$ diquark, and the size of
the $(us)$ diquark for the five schemes. As can be seen in the differences in
the bindings, for all cases, the $(us)$ diquark has a stronger attraction than
the $(ud)$ diquark. We further find that within the uncertainties given in the
present analysis, the $(us)$ diquark mass could be between 580 to 970 MeV.
Subsequently, in the following section, we will assess the influence of the
$(us)$ diquark using the coalescence model and incorporating the estimated
range of $(us)$ diquark masses.

###  II.2 Diquarks and hadron overlap

The strong correlation within the $(us)$ diquark leads to increased production
when a light quark $q$ and an $s$-quark are in close proximity; this
correlation is stronger. These effects should be incorporated into the
modeling of hadron production. As previously discussed, the light diquark
within the baryon triplets $\Lambda_{c}$ and $\Xi_{c}$ is primarily the so-
called ’good diquark,’ characterized by an anti-triplet color configuration,
zero spin, and anti-triplet flavor. An explicit quark model calculation
reveals that the SU(3) flavor symmetry breaking in $\Xi_{c}$ results in a
mixing of less than 0.01% from the spin-1 diquark component. This implies that
we can indeed assume the $(us)$ diquark within the $\Xi_{c}$ to be a good
diquark, denoted as $(us)_{S=0}$ with the subscript indicating a spin of 0.
Assuming multiple quarks are present in the hadronization process, the strong
correlation between a light quark $q$ and a strange quark $s$ characterized by
$(qs)_{S=0}$ will provide an additional production mechanism for $\Xi_{c}$.

It should be noted that the strong correlation of $(qs)_{S=0}$ has a smaller
effect on the production of baryon octets with strange quarks as the overlaps
are smaller. This is because the spin-1/2 baryon octet is composed of a mixed
flavor symmetry so that any two quarks within a baryon will be in either a
color anti-triplet and flavor symmetric configuration or a color anti-triplet
and flavor antisymmetric configuration with equal probability. In other words,
the quark model predicts that the probability of a diquark being in the ”good”
diquark state inside any baryon octet state is 1/2, while inside a heavy
baryon, it is 1. Therefore, focusing on the $(qs)$ diquark, we find the
following constraint.

$\displaystyle\sum_{q=u,d}|\langle(sq)_{{S=0}}|B_{8}\rangle|^{2}$
$\displaystyle\leq$ $\displaystyle\frac{1}{2},$ (6)
$\displaystyle\sum_{q=u,d}|\langle(sq)_{{S=0}}|\Xi_{c}\rangle|^{2}$
$\displaystyle=$ $\displaystyle 1,$ (7)

where $|B_{8}\rangle$ ($|\Xi_{c}\rangle$ ) means any diquark in the flavor
octet baryon ($\Xi_{c}$) state. Table 5 shows how all the probabilities of
diquarks inside the baryon octet states are distributed. For example,
$|\langle(us)_{{S=0}}|\Sigma^{+}\rangle|^{2}=\frac{1}{2}$ and
$|\langle(us)_{{S=0}}|\Sigma^{0}\rangle|^{2}=|\langle(ds)_{{S=0}}|\Sigma^{0}\rangle|^{2}=\frac{1}{4}$,
while $|\langle(us)_{{S=0}}|\Lambda\rangle|^{2}=\frac{1}{6}$. Also,
$|\langle(qs)_{{S=0}}|\Xi\rangle|^{2}=\frac{1}{2}$. Consequently, the presence
of a $(us)_{S=0}$ will have a lesser impact on the production of hyperons
compared to the anticipated effect on the production of $\Xi_{c}$. Thus, we
will solely focus on the supplementary production of $\Xi_{c}$ using the
coalescence model.

Table 5: Probabilities of each diquark component for the baryon octet states are depicted. Each diquark configuration is represented in the second row by the respective color-spin factor, considering the quark masses of the two quarks $i$ and $j$ in the diquark as $-\frac{1}{m_{i}m_{j}}\lambda_{i}^{c}\lambda_{j}^{c}\sigma_{i}\cdot\sigma_{j}$. The color-spin factors in the second column are determined by multiplying the probability of a diquark with the corresponding color-spin factors containing mass terms multiplied by 3. These values are summed across all diquarks listed in the row for the corresponding baryon |  | Good diquarks | Bad diquarks
---|---|---|---
Baryon | Color-spin factor | $-\frac{8}{m_{q}^{2}}$ | $-\frac{8}{m_{q}m_{s}}$ | $-\frac{8}{m_{s}^{2}}$ | $\frac{8}{3m_{q}^{2}}$ | $\frac{8}{3m_{q}m_{s}}$ | $\frac{8}{3m_{s}^{2}}$
$p,n$ | $-\frac{8}{m_{q}^{2}}$ | $\frac{1}{2}$ | 0 | 0 | $\frac{1}{2}$ | 0 | 0
$\Lambda$ | $-\frac{8}{m_{q}^{2}}$ | $\frac{1}{3}$ | $\frac{1}{6}$ | 0 | 0 | $\frac{1}{2}$ | 0
$\Sigma$ | $\frac{8}{3m_{q}^{2}}-\frac{32}{3m_{q}m_{s}}$ | 0 | $\frac{1}{2}$ | 0 | $\frac{1}{3}$ | $\frac{1}{6}$ | 0
$\Xi$ | $\frac{8}{3m_{s}^{2}}-\frac{32}{3m_{q}m_{s}}$ | 0 | $\frac{1}{2}$ | 0 | 0 | $\frac{1}{6}$ | $\frac{1}{3}$

###  II.3 Diquarks with flow

To incorporate the strong correlation in the $(qs)$ diquark within a
coalescence model, we will assume the presence of $(qs)$ diquarks in the
quark-gluon plasma[23] and estimate their total number using the thermal
model. Whether to explicitly include all or part of the diquark configurations
is a subtle question. In coalescence models of hadron production from the
quark-gluon plasma, the 3-body coalescence formula for $\Lambda_{c}$
production, as presented in Refs. [7, 8], involves phase space integrals of
two relative coordinates. Consequently, it can be reformulated as a quark-
quark to diquark 2-body coalescence, along with an additional diquark-heavy
quark 2-body coalescence. Hence, with suitable normalization factors, one can
effectively replace 3-body coalescence with the 2-body coalescence with
diquarks if there is no binding energy. This seems to be consistent with
coalescence approaches, where the bulk part of the heavy baryon production
comes from the recombination of the heavy quark and a surrounding diquark[9,
10]. What we want to emphasize is that there is an additional binding in the
(qs) diquark, which has to be taken into account through the existence of
bound diquarks, when using the coalescence model[23], to explain the missing
strength in the $\Xi_{c}/D^{0}$ production.

The expected diquark number at central rapidity can be estimated using the
statistical model at the hadronization temperature $T_{H}$ [24].

$\displaystyle N^{stat}_{[us]}$
$\displaystyle=V_{H}\frac{g_{[us]}}{2\pi^{2}}\int^{\infty}_{0}\frac{p^{2}dp}{\gamma^{-1}_{h}e^{E_{h}/T_{H}}\pm
1},$ (8)

where $V_{H}$ is the volume at the chemical freeze-out point. Later, we will
use the limiting volume $V_{H}=20$ fm3 extracted for a small system in pp
collision at 7 TeV [25]. We further augment the distribution with the flow.
The collective radial expansion of the fireball, created in heavy-ion
collisions can be understood well within a hydrodynamic picture. A more
phenomenological way to capture this isotropic expansion is known as blast-
wave model [26]. The model assumes a spectrum of purely thermal sources which
are boosted in transverse direction. The velocity distribution in $0\leq r\leq
R$ is assumed to be

$\displaystyle\beta_{r}=\left(\frac{r}{R}\right)^{n}\beta_{s},$ (9)

where $\beta_{s}$ is the surface velocity, a free parameter of the fit. In
many applications, a linear profile is assumed and $n$ is fixed equal to
unity. The quality of the fit can be improved if $n$ is considered as an
additional free parameter. The resulting values for the kinetic freeze-out
temperature $T_{kin}$ and $\beta_{s}$ are generally anti-correlated. The so-
obtained spectral shape is a superposition of the contributions due to the
individual thermal sources and is given by

$\frac{1}{m_{\rm T}}\frac{{\rm d}N}{{\rm d}m_{\rm T}}\propto m_{\rm
T}\int_{0}^{R}I_{0}\left(\frac{p_{\rm
T}\sinh\rho}{T_{kin}}\right)K_{1}\left(\frac{m_{\rm
T}\cosh\rho}{T_{kin}}\right)r\,{\rm d}r\,,$ (10)

where $I_{0}(x)$ and $K_{1}(x)$ are Bessel functions, $m_{\rm T}=\sqrt{p_{\rm
T}^{2}+m^{2}}$ and $\rho=\tanh^{-1}\beta_{r}$.

Figure 1:  Normalised blast-wave distributions for three different possible
diquark masses (680 MeV, 826 MeV and 970 MeV) and for three collisions systems
(pp minimum bias, p–Pb non-single diffractive and central Pb–Pb).

Blast-wave fits to existing data have been performed by the ALICE
Collaboration at the LHC for different collision energies and collision
systems, spanning a large range of charged particle multiplicities per
pseudorapidity $\langle{\mathrm{d}}N_{\mathrm{c}h}/{\mathrm{d}}\eta\rangle$
per event. A summary of this is discussed in Ref. [27]. The typical approach
is to fit the transverse-momentum spectra of charged $\pi$, K, and p
altogether and by that extract a common freeze-out temperature $T_{kin}$ and
mean velocity $\langle\beta\rangle$. In addition, a shape factor $n$ is used.
This set of common blast-wave parameters can be used to predict the shape of
spectra of unmeasured particles (spectra). It has been used by the ALICE
Collaboration to encounter possible biases in case only parts of the whole
$p_{\mathrm{T}}$ spectrum was accessible [28] or something similar [29]. We
use it here to predict the transverse momentum shape for calculations where
flow is not included. The normalized blast-wave function can be multiplied
onto the corresponding number for the constituent with the same mass to obtain
the transverse momentum distribution. The parameters for the blast-wave
functions for four multiplicity values, corresponding to minimum bias pp
collision at $\sqrt{s}=13$ TeV, minimum bias pp collision at $\sqrt{s}=5$ TeV,
non-single diffractive p–Pb at $\sqrt{s_{\mathrm{NN}}}=5$ TeV and central
Pb–Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5$ TeV, are given in Table 6.

Collision system | $\langle{\mathrm{d}}N_{\mathrm{c}h}/{\mathrm{d}}\eta\rangle$ |  $T_{kin}$ (GeV) | $\langle\beta\rangle$ | n
---|---|---|---|---
pp 13 TeV MB | 6.88 | 0.184 | 0.270 | 3.878
pp 5.02 TeV MB | 4.30 | 0.181 | 0.198 | 6.248
p–Pb 5.02 TeV NSD | 17.81 | 0.177 | 0.423 | 1.846
Pb–Pb 5.02 TeV 0-10% | 1781 | 0.113 | 0.659 | 0.650
Table 6: Values used for the blast-wave predictions shown in Fig. 1

We can now predict the transverse momentum shape of the diquarks with
different masses in the blast-wave model with parameters given in Table 6 for
different multiplicities. The results are shown in Fig. 1 for three different
diquark masses.

##  III $\Xi_{c}/D^{0}$ ratios

To estimate the additional enhancement of $\Xi_{c}/D^{0}$ ratio coming from
$(us)$ diquark, we use the 2-dimensional coalescence model.

$\displaystyle\frac{\mathrm{d}^{2}N_{\Xi_{c}^{+}}}{\mathrm{d}^{2}P_{\mathrm{T}}}=\frac{g_{\Xi_{c}^{+}}}{g_{[us]}g_{c}}$
$\displaystyle\int
d^{2}x_{1}d^{2}x_{2}d^{2}p_{1\mathrm{T}}d^{2}p_{2\mathrm{T}}\
\frac{\mathrm{d}^{2}N_{[us]}}{A_{L}\mathrm{d}^{2}p_{[us]\mathrm{T}}}\frac{\mathrm{d}^{2}N_{c}}{A_{L}\mathrm{d}^{2}p_{c\mathrm{T}}}\times
W_{\Xi_{c}^{+}}(\vec{r},\vec{k})\delta^{(2)}(\vec{P}_{\mathrm{T}}-\vec{p}_{[us]\mathrm{T}}-\vec{p}_{c\mathrm{T}})$
$\displaystyle=\frac{g_{\Xi_{c}^{+}}}{g_{[us]}g_{c}}$
$\displaystyle\left(2\sqrt{\pi}\sigma\right)^{2}\frac{1}{A}\int
d^{2}p_{1\mathrm{T}}d^{2}p_{2\mathrm{T}}\
\frac{\mathrm{d}^{2}N_{[us]}}{\mathrm{d}^{2}p_{[us]\mathrm{T}}}\frac{\mathrm{d}^{2}N_{c}}{\mathrm{d}^{2}p_{c\mathrm{T}}}\times\exp{\left[-\sigma^{2}p^{\prime
2}\right]}\delta^{(2)}(\vec{P}_{\mathrm{T}}-\vec{p}_{[us]\mathrm{T}}-\vec{p}_{c\mathrm{T}}),$
(11)

where the $g$ values and $A$ are statistical factors and the coalescence area
at freeze-out point. In this model, we used the following two-dimensional
Gaussian-type Wigner function.

$\displaystyle W_{\Xi_{c}^{+}}=4\exp{\left(-\frac{r^{\prime
2}}{\sigma^{2}}-\sigma^{2}k^{\prime 2}\right)}.$ (12)

The parameter $\sigma$ is related to the $\Xi_{c}^{+}$ radius by
$\sigma=\sqrt{8/3}r_{\Xi_{c}^{+}}$ and we used $r_{\Xi_{c}^{+}}=0.222$ fm
obtained by quark model calculation. Furthermore, the primed momenta in
Eq.(11) and Eq.(12) are taken in the center-of-mass frame of the $\Xi_{c}$, as
discussed in Ref. [8].

To determine $A$, we use the method in Ref. [30]. The yield of d and 3He in
PbPb collisions are determined at the chemical freeze-out point [1].
Therefore, we can obtain the coalescence area $A$ at the chemical freeze-out
point by fitting the production data at pp collisions using the same
coalescence formula given in Eq. (11) for the deuteron and the corresponding
three-body formula for 3He. It should be noted that in the case of pp
collisions, where the kinetic freeze-out temperature exceeds the critical
temperature, we consider the chemical freeze-out temperature to be the same as
the kinetic freeze-out temperature.

Figure 2: Transverse momentum distribution of d (left panel) and 3He (right
panel) measured by the ALICE Collaboration in pp collisions at $\sqrt{s}=7$
TeV [31]. The blue curves are descriptions from the two-dimensional
coalescence model described in the text.

Fig. 2 shows coalescence model results for the d and 3He in pp collisions at
$\sqrt{s}=7$ TeV. In the model calculation, we also accounted for the
correction arising from the relatively smaller system size in pp collisions
compared to that of d and 3He[32]. The correction factor can be derived from
experimental data on the ratios of d/p and 3He/p between different collision
systems [31, 33]. The correction factor for deuteron to proton in pp
collisions at $\sqrt{s}=7$ TeV
($\left<\mathrm{d}N_{ch}/\mathrm{d}\eta\right>=6.01$) is 0.4065[31], while
that for 3He/p at similar multiplicity is 0.13 [34]. These correction factors
are multiplied after the model calculation using the following method.

We used measured proton transverse momentum distribution in pp collisions at
$\sqrt{s}=7$ TeV [35] multiplied by $r_{f}$ to remove the feed-down
contribution in the formation of the deuteron and 3He. For evaluating the
feed-down fraction, we employed the statistical hadronization model (as
expressed in Eq. 8) at a chemical freeze-out temperature of $181$ MeV. This
temperature aligns with the freeze-out temperature deduced from the Blast-wave
fit for pp collisions at 7 TeV and is anticipated to be analogous at 5.02 TeV.
Then, one finds that $r_{f}=$26.6$\%$ of protons participate in the
coalescence. The outcome achieved at $T=181$ MeV lies within the range of
uncertainties stemming from our range of diquark masses

Furthermore, we employ $A_{\rm{pp}}^{7\ \rm{TeV}}=2.52$ fm2. This value is
derived by scaling the coalescence area in PbPb collisions at
$\sqrt{s_{\rm{NN}}}=2.76$ TeV ($A_{\rm{PbPb}}^{2.76\ \rm{TeV}}$=608 fm2 [30])
with the charged particle multiplicity ratio between pp collisions at
$\sqrt{s}=7$ TeV and PbPb collisions at $\sqrt{s_{\rm{NN}}}=2.76$ TeV. The
charged particle multiplicity is
$\left<\mathrm{d}N_{ch}/\mathrm{d}\eta\right>=6.01$ in pp collision at
$\sqrt{s}=7$ TeV [36] and
$\left<\mathrm{d}N_{ch}/\mathrm{d}\eta\right>=1447.5\pm 39$ in PbPb collisions
at $\sqrt{s_{\rm{NN}}}=2.76$ TeV 0-10$\%$ event [37].

The consistency of our approach can be justified by the model description
shown in Fig. 2, where we could have determined the two independent parameters
$r_{f}$ and $A$ by fitting the experimental data for the deuteron and 3He,
which turns out to yield almost the same values. Hence, we will use $A=2.52$
fm2 in Eq. (11).

It is important to observe that both $V_{H}$ and $A$ scale in relation to the
multiplicity, which is empirically proportional to $s^{0.103}$[2]. However,
these quantities are used in Eq. (11) as the ratio $V_{H}/A$. Consequently,
even though we established this ratio at 7 TeV, its value is anticipated to be
comparably consistent across pp collisions at both 5.02 and 13 TeV.

###  III.1 Charm quark transverse momentum distribution

To obtain the transverse momentum distribution of charm quark in pp collisions
at $\sqrt{s}=5.02$ TeV, we used the transverse momentum distribution of
$D^{0}$ meson [38] measured at the same collision energy. First, we assume
that the shape of the charm quark distribution is the same as that of the
$D^{0}$ meson. We then normalized this distribution using the ratio of the
charm quark cross section $d\sigma^{c\bar{c}}/dy_{|y|<0.5}=1165$ $\mu$b [39]
to fitted $D^{0}$ meson cross-section $d\sigma^{D^{0}}/dy_{|y|<0.5}=447$
$\mu$b [38] . Figure (3) shows the measured $D^{0}$ meson spectrum, the
corresponding fit (left) and the charm quark distribution (right).

Figure 3: Fitted $p_{T}$ distribution of $D^{0}$ and charm quark.

###  III.2 Diquark distribution

For the $(us)$ diquark $p_{T}$ distribution we adopt the blast-wave
distribution. The shape of the diquark distribution was obtained from the
blast-wave model given in Eq. (10) with the normalization constant $\alpha$ to
fit the statistical model prediction for the total number after integrating
over the transverse momentum. Finally, the $p_{T}$ distribution of
$\Xi_{c}^{0}$ composed of the $c$ quark and $(us)$ diquark becomes

$\displaystyle\frac{\mathrm{d}^{2}N_{\Xi_{c}^{+}}}{\mathrm{d}^{2}P_{\mathrm{T}}\mathrm{d}y}=\frac{g_{\Xi_{c}^{+}}}{g_{[us]}g_{c}}\left(2\sqrt{\pi}\sigma\right)^{2}\frac{\alpha}{2\pi
A}\ \int_{0}^{R}r\,{\rm d}r\int d^{2}p_{1\mathrm{T}}\ m_{1\mathrm{T}}\times
I_{0}\left(\frac{p_{1\mathrm{T}}\sinh\rho}{T_{kin}}\right)K_{1}\left(\frac{m_{1\mathrm{T}}\cosh\rho}{T_{kin}}\right)\int
d^{2}p_{2\mathrm{T}}\frac{\mathrm{d}^{2}N_{c}}{\mathrm{d}^{2}p_{c\mathrm{T}}\mathrm{d}y}$
$\displaystyle\times\exp{\left[-\sigma^{2}p^{\prime
2}\right]}\delta^{(2)}(\vec{P}_{\mathrm{T}}-\vec{p}_{\mathrm{1T}}-\vec{p}_{c\mathrm{T}}).$
(13)

Here, we use $A$=2.52 fm2 and $V_{H}=20$ fm3 determined at $\sqrt{s}=7$ TeV.
However, since both quantities scale as multiplicity, we use the same ratio
$V_{H}/A$ to calculate the $p_{\rm{T}}$ distribution of $\Xi_{c}$ in pp
collisions at $\sqrt{s}=5.02$ TeV. Furthermore, it is important to mention
that the smaller size of $\Xi_{c}$ allows us to disregard the correction
factor arising from the small coalescence size effect as it amounts to less
than 2$\%$ correction.

###  III.3 Result

####  III.3.1 feed-down

Since we are assuming $(qs)$ diquarks with spin 0 in the quark-gluon plasma,
the coalescence with charm quark with non-zero relative orbital angular
momentum $L$ will produce excited states with higher total spin. The
coalescence of higher orbitals is sequentially suppressed, as shown in Ref.
[40]. Therefore, we have considered the contribution only from states
generated with $L=1$, which will produce $\Xi_{c}$ excited states with
negative parity and total spin 1/2 and 3/2, both of which are three-star
states. These states will contribute to the production of the $\Xi_{c}$ states
through the feed-down process via the $\Xi_{c}^{\prime}$ intermediate state,
which eventually decay to $\Xi_{c}$ electromagnetically.

The yield ratios for the production of these excited states to the ground
state $\Xi_{c}$ can be calculated by either the statistical hadronization
model or the coalescence model, which results in similar ratios[40]. Here, we
use the ratios given by the statistical model at several temperatures used in
the result section as given in Table 7.

$T$ (MeV) | $N_{\Xi_{c}(2790)}/N_{\Xi_{c}}$ | $N_{\Xi_{c}(2815)}/N_{\Xi_{c}}$ | Total
---|---|---|---
165 | 0.1707 | 0.2971 | 0.4678
181 | 0.2024 | 0.3569 | 0.5593
184 | 0.2078 | 0.3673 | 0.5751
Table 7: Feed-down fraction for selected chemical freeze-out temperatures.

These ratios are added to the total production of the $\Xi_{c}$ state. Since
coalescence from higher orbitals as well as higher nodal modes will be further
suppressed, we can neglect contributions from any other possible excited
states at this stage.

####  III.3.2 Total contribution

In Figure 4, we present the transverse momentum distribution of the
$\Xi_{c}^{0}$/$D^{0}$ ratio. The $p_{\mathrm{T}}$ distribution of $\Xi_{c}$
comprising a $(us)$ diquark and a charm quark is computed using Eq. (13) for
various $(us)$ diquark masses permitted by different schemes. Experimental
data is employed for the $D^{0}$ meson, while the upper and lower bounds are
derived from the associated experimental uncertainties in its distribution.

To quantify the systematic uncertainty, we consider several scenarios for the
chemical and kinetic freeze-out temperatures, denoted as $T_{c}\sim T_{H}$ and
$T_{K}$, respectively. Additionally, we explore cases with different diquark
masses. The temperature for kinetic decoupling is generally different from
that for the phase transition where parton coalescence takes place because
hadronic scattering delays the kinetic decoupling after hadronization. But if
the size of nuclear matter, or in other words, the particle multiplicity is
small as in pp collisions, the effects of hadronic scattering will be less
important. This is supported by the blast wave model, where the temperature
for the kinetic freeze-out $T_{K}$ is not far from $T_{c}$. Nevertheless, to
probe the uncertainties associated with the effective kinetic freeze-out
temperature, we compare the two possible kinetic freeze-out temperatures
$T_{K}=184$ and 181 MeV as given in Table 6 for pp collision depending on
collision energy. As can be seen in the upper panel of Figure 4, the fit with
a slightly larger temperature has a higher flow velocity, consequently
shifting the distribution slightly towards higher $p_{T}$. In the lower
figure, we consider three possible scenarios for the pp collision: i) We first
take $m_{[qs]}=580$ MeV and $T_{c}=T_{K}=181$ MeV as given by the blast-wave
model. ii) $T_{C}=165$ MeV and $T_{K}=181$ MeV as given in Ref. [41] and the
blast-wave model, respectively. iii) To probe the uncertainty in the diquark
mass, we take $m_{[qs]}=680$ MeV and $T_{c}=T_{K}=181$ MeV. The lower Figure 4
shows the transverse momentum distribution for the $\Xi_{c}/D^{0}$ ratio
calculated for the possible scenarios. As can be seen in the plot, the only
noticeable difference comes from the values of $T_{c}$ and $m_{[qs]}$, which
lead to variations in the overall magnitude of the production. One can observe
the characteristic peak structure remains around $p_{T}\sim 2$ GeV.

Lastly, in Figure 5, we depict the outcome for a $(us)$ diquark mass of 580
MeV and $T_{K}=181$ MeV with two different $T_{c}=181$ and 165 MeV after
incorporating the production ratios calculated in Ref. [41]. It should be
noted that in Refs. [41] and [8], the overall normalization in the coalescence
formula for charmed hadrons is fixed by requiring that all zero transverse
momentum charm quarks contribute to the production of the charmed hadrons.
Therefore, if there is additional production of $\Xi_{c}$ originating from
$(qs)$ diquarks, the normalization must be readjusted. However, charm
coalescence is primarily driven by the production of $D$, $D^{*}$, $D_{s}$,
$D_{s}^{*}$, and $\Lambda_{c}$ as discussed in Ref. [40]. Consequently, we did
not modify the normalization factor for the $\Xi_{c}/D^{0}$ ratio from Ref.
[41] when incorporating it into Figure 5. Furthermore, while these
calculations are derived for $\sqrt{s}=5.02$ TeV, any corrections anticipated
in higher-energy pp collisions are likely to be negligible for the ratios.

As observed in the illustration, assuming a diquark mass of 580 MeV, the
introduced contribution seems to provide the additional strength needed to
reproduce the experimental observation.

Figure 4: $\Xi_{c}/D^{0}$ ratio due to $(qs)$ diquark with different masses .
The upper panel shows the effects of different $T_{K}$, while the lower panel
shows different combinations of $T_{K}$, $T_{c}$ and diquark mass. Although
the ratio was obtained for 5.02 TeV, the value will only be slightly different
for 13 or 7 TeV. Figure 5: $\Xi_{c}/D^{0}$ ratio after adding the enhancement
due to the $(qs)$ diquark with mass 580 MeV to the coalescence and
fragmentation calculation in Ref. [41].

##  IV Discussion and summary

We have demonstrated that within a quark model employing parameters that
effectively reproduce ground state hadron masses, a $(us)$ diquark in the
color anti-triplet spin-0 channel exhibits stronger binding compared with a
$(ud)$ diquark in the same color and spin state.

This heightened attraction persists even when the potentials and constituent
quark masses are adjusted to approximate values near the critical temperature.
The presence of such robust $(us)$ or $(ds)$ diquarks is anticipated to
amplify the production of charmed baryons where the dominant diquark
configurations align with these ”good” diquarks.

Furthermore, using the coalescence model, we have illustrated that such
diquark correlations might potentially amplify $\Xi_{c}$ production, thereby
offering a plausible explanation for the observed anomalous enhancement in
$\Xi_{c}/D^{0}$ production ratios within high-energy pp collisions. This
effect is postulated to arise from the existence of a $(us)$ diquark with an
approximate mass of 580 MeV at the point of hadronization.

Analogous enhancements are predicted in pPb or PbPb collisions. Comprehensive
estimations will be detailed in an upcoming publication.

## Acknowledgements

This work was supported by Samsung Science and Technology Foundation under
Project Number SSTF-BA1901-04, and by the Korea National Research Foundation
under grant number No. 2023R1A2C3003023. B.D. acknowledges the support from
Bundesministerium für Bildung und Forschung through ErUM-FSP T01
(Förderkennzeichen 05P21RFCA1). The work of A.P. was supported by the Korea
National Research Foundation under the grant number 2021R1I1A1A01043019. S.
Lim acknowledges support from the National Research Foundation of Korea grants
funded by the Korean government under project number NRF-2008-00458.

##  Appendix A Parametrization of the transverse momentum distribution

An exponential function and a power-law type function were used for low
$p_{\mathrm{T}}$ and high $p_{\mathrm{T}}$, respectively. Here, $p_{1}=1$ GeV
and $|y|<0.5$.

$\displaystyle\frac{\mathrm{d}N}{\mathrm{d}^{2}p_{\mathrm{T}}\mathrm{d}y}$
$\displaystyle=ae^{-b(p_{\mathrm{T}}/p_{1})^{c}}\qquad(p_{\mathrm{T}}<p_{c}),$
$\displaystyle=\frac{d}{\left[1+(p_{\mathrm{T}}/p_{0})^{2}\right]^{e}}\qquad(p_{\mathrm{T}}>p_{c})$
(14)

Particle |  $\sqrt{s}$ (GeV) | a (GeV)-2 | b | c | d (GeV)-2 | e |  $p_{0}$ (GeV) |  $p_{c}$ (GeV)
---|---|---|---|---|---|---|---|---
pl | 7 | $6.733\times 10^{-2}$ | 1.904 | 1.456 | $3.477\times 10^{-2}$ | 3.421 | 1.436 | 1.2
pu | 7 | $7.738\times 10^{-2}$ | 1.829 | 1.404 | $4.546\times 10^{-2}$ | 3.462 | 1.438 | 1.2
$D^{0}_{l}$ | 5.02 | $1.004\times 10^{-3}$ | 0.8601 | 1.156 | 4.851$\times 10^{-4}$ | 3.067 | 2.864 | 2.0
$D^{0}_{u}$ | 5.02 | $1.360\times 10^{-3}$ | 0.7729 | 1.172 | 8.940 $\times 10^{-4}$ | 2.967 | 2.644 | 2.0
Table 8: The fitting parameters of p and $D^{0}$ in pp collisions. The
subscripts in the first column correspond to the fits for the upper bound
($u$) and lower bound ($l$).

## References

  * Andronic  _et al._ [2018] A. Andronic, P. Braun-Munzinger, K. Redlich, and J. Stachel, Nature 561, 321 (2018), arXiv:1710.09425 [nucl-th] .
  * ALI [2022] (2022), arXiv:2211.04384 [nucl-ex] .
  * Acharya  _et al._ [2023] S. Acharya _et al._ (ALICE), Phys. Lett. B 839, 137796 (2023), arXiv:2112.08156 [nucl-ex] .
  * Acharya  _et al._ [2021a] S. Acharya _et al._ (ALICE), JHEP 10, 159 (2021a), arXiv:2105.05616 [nucl-ex] .
  * Acharya  _et al._ [2021b] S. Acharya _et al._ (ALICE), Phys. Rev. Lett. 127, 272001 (2021b), arXiv:2105.05187 [nucl-ex] .
  * Greco  _et al._ [2003] V. Greco, C. M. Ko, and P. Levai, Phys. Rev. Lett. 90, 202302 (2003), arXiv:nucl-th/0301093 .
  * Plumari  _et al._ [2018] S. Plumari, V. Minissale, S. K. Das, G. Coci, and V. Greco, Eur. Phys. J. C 78, 348 (2018), arXiv:1712.00730 [hep-ph] .
  * Cho  _et al._ [2020] S. Cho, K.-J. Sun, C. M. Ko, S. H. Lee, and Y. Oh, Phys. Rev. C 101, 024909 (2020), arXiv:1905.09774 [nucl-th] .
  * Beraudo  _et al._ [2022] A. Beraudo, A. De Pace, M. Monteno, M. Nardi, and F. Prino, Eur. Phys. J. C 82, 607 (2022), arXiv:2202.08732 [hep-ph] .
  * Beraudo  _et al._ [2024] A. Beraudo, A. De Pace, D. Pablos, F. Prino, M. Monteno, and M. Nardi, Phys. Rev. D 109, L011501 (2024), arXiv:2306.02152 [hep-ph] .
  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Andronic  _et al._ [2021] A. Andronic, P. Braun-Munzinger, M. K. Köhler, A. Mazeliauskas, K. Redlich, J. Stachel, and V. Vislavicius, JHEP 07, 035 (2021), arXiv:2104.12754 [hep-ph] .
  * Minissale  _et al._ [2023] V. Minissale, S. Plumari, and V. Greco, EPJ Web Conf. 276, 02013 (2023).
  * Huovinen  _et al._ [2017] P. Huovinen, P. M. Lo, M. Marczenko, K. Morita, K. Redlich, and C. Sasaki, Phys. Lett. B 769, 509 (2017), arXiv:1608.06817 [hep-ph] .
  * Anselmino  _et al._ [1993] M. Anselmino, E. Predazzi, S. Ekelin, S. Fredriksson, and D. B. Lichtenberg, Rev. Mod. Phys. 65, 1199 (1993).
  * Park  _et al._ [2019] W. Park, S. Noh, and S. H. Lee, Nucl. Phys. A 983, 1 (2019), arXiv:1809.05257 [nucl-th] .
  * Noh  _et al._ [2021] S. Noh, W. Park, and S. H. Lee, Phys. Rev. D 103, 114009 (2021), arXiv:2102.09614 [hep-ph] .
  * Park  _et al._ [2017] W. Park, A. Park, S. Cho, and S. H. Lee, Phys. Rev. D 95, 054027 (2017), arXiv:1702.00381 [hep-ph] .
  * Karliner and Rosner [2014] M. Karliner and J. L. Rosner, Phys. Rev. D 90, 094007 (2014), arXiv:1408.5877 [hep-ph] .
  * Karliner and Rosner [2017] M. Karliner and J. L. Rosner, Phys. Rev. Lett. 119, 202001 (2017), arXiv:1707.07666 [hep-ph] .
  * Lafferty and Rothkopf [2020] D. Lafferty and A. Rothkopf, Phys. Rev. D 101, 056010 (2020), arXiv:1906.00035 [hep-ph] .
  * Gubler  _et al._ [2020] P. Gubler, T. Song, and S. H. Lee, Phys. Rev. D 101, 114029 (2020), arXiv:2003.09073 [hep-ph] .
  * Lee  _et al._ [2008] S. H. Lee, K. Ohnishi, S. Yasui, I.-K. Yoo, and C.-M. Ko, Phys. Rev. Lett. 100, 222301 (2008), arXiv:0709.3637 [nucl-th] .
  * Vovchenko  _et al._ [2019] V. Vovchenko, B. Dönigus, and H. Stoecker, Phys. Rev. C 100, 054906 (2019), arXiv:1906.03145 [hep-ph] .
  * Mazeliauskas and Vislavicius [2020] A. Mazeliauskas and V. Vislavicius, Phys. Rev. C 101, 014910 (2020), arXiv:1907.11059 [hep-ph] .
  * Schnedermann  _et al._ [1993] E. Schnedermann, J. Sollfrank, and U. W. Heinz, Phys. Rev. C 48, 2462 (1993), arXiv:nucl-th/9307020 .
  * Acharya  _et al._ [2020] S. Acharya _et al._ (ALICE), Eur. Phys. J. C 80, 693 (2020), arXiv:2003.02394 [nucl-ex] .
  * Acharya  _et al._ [2018a] S. Acharya _et al._ (ALICE), Nucl. Phys. A 971, 1 (2018a), arXiv:1710.07531 [nucl-ex] .
  * Adam  _et al._ [2016a] J. Adam _et al._ (ALICE), Phys. Lett. B 752, 267 (2016a), arXiv:1506.07499 [nucl-ex] .
  * Yun  _et al._ [2023] H. Yun, D. Park, S. Noh, A. Park, W. Park, S. Cho, J. Hong, Y. Kim, S. Lim, and S. H. Lee, Phys. Rev. C 107, 014906 (2023), arXiv:2208.06960 [hep-ph] .
  * Acharya  _et al._ [2018b] S. Acharya _et al._ (ALICE), Phys. Rev. C 97, 024615 (2018b), arXiv:1709.08522 [nucl-ex] .
  * Sun  _et al._ [2019] K.-J. Sun, C. M. Ko, and B. Dönigus, Phys. Lett. B 792, 132 (2019), arXiv:1812.05175 [nucl-th] .
  * Adam  _et al._ [2016b] J. Adam _et al._ (ALICE), Phys. Rev. C 93, 024917 (2016b), arXiv:1506.08951 [nucl-ex] .
  * Acharya  _et al._ [2022a] S. Acharya _et al._ (ALICE), JHEP 01, 106 (2022a), arXiv:2109.13026 [nucl-ex] .
  * Adam  _et al._ [2015] J. Adam _et al._ (ALICE), Eur. Phys. J. C 75, 226 (2015), arXiv:1504.00024 [nucl-ex] .
  * Aamodt  _et al._ [2010] K. Aamodt _et al._ (ALICE), Eur. Phys. J. C 68, 345 (2010), arXiv:1004.3514 [hep-ex] .
  * Abelev  _et al._ [2013] B. Abelev _et al._ (ALICE), Phys. Rev. C 88, 044910 (2013), arXiv:1303.0737 [hep-ex] .
  * Acharya  _et al._ [2019] S. Acharya _et al._ (ALICE), Eur. Phys. J. C 79, 388 (2019), arXiv:1901.07979 [nucl-ex] .
  * Acharya  _et al._ [2022b] S. Acharya _et al._ (ALICE), Phys. Rev. D 105, L011103 (2022b), arXiv:2105.06335 [nucl-ex] .
  * Cho  _et al._ [2017] S. Cho _et al._ (ExHIC), Prog. Part. Nucl. Phys. 95, 279 (2017), arXiv:1702.00486 [nucl-th] .
  * Minissale  _et al._ [2021] V. Minissale, S. Plumari, and V. Greco, Phys. Lett. B 821, 136622 (2021), arXiv:2012.12001 [hep-ph] .