# Two-dimensional easy-plane SU$(3)$ magnet with the transverse field:
Anisotropy-driven multicriticality

Yoshihiro Nishiyama  Department of Physics, Faculty of Science, Okayama
University, Okayama 700-8530, Japan

###### Abstract

The two-dimensional easy-plane SU$(3)$ magnet subjected to the transverse
field was investigated with the exact-diagonalization method. So far, as to
the $XY$ model (namely, the easy-plane SU$(2)$ magnet), the transverse-field-
driven order-disorder phase boundary has been investigated with the exact-
diagonalization method, and it was claimed that the end-point singularity
(multicriticality) at the $XX$-symmetric point does not accord with
large-$N$-theory’s prediction. Aiming to reconcile the discrepancy, we extend
the internal symmetry to the easy-plane SU$(3)$ with the anisotropy parameter
$\eta$, which interpolates the isotropic ($\eta=0$) and fully anisotropic
($\eta=1$) cases smoothly. As a preliminary survey, setting $\eta=1$, we
analyze the order-disorder phase transition through resorting to the fidelity
susceptibility $\chi_{F}$, which exhibits a pronounced signature for the
criticality. Thereby, with $\eta$ scaled carefully, the $\chi_{F}$ data are
cast into the crossover-scaling formula so as to determine the crossover
exponent $\phi$, which seems to reflect the extension of the internal symmetry
group to SU(3).

##  1 Introduction

The one-dimensional $XY$ model subjected to the transverse field $H$ and
anisotropy $\eta$ with the Hamiltonian ${\cal
H}_{XY}=-\sum_{i}[(1+\eta)\sigma^{x}_{i}\sigma^{x}_{i+1}+(1-\eta)\sigma^{y}_{i}\sigma^{y}_{i+1}+H\sigma^{z}_{i}]$
(${\mathbf{\sigma}}_{i}$: Pauli matrices at site $i$) is attracting much
attention [1, 2, 3, 4] in the context of the quantum information theory [5,
6]. A key ingredient is that the model covers both $XX$\- ($\eta=0$) and
Ising-symmetric ($\eta=1$) cases, and there appear rich characters as to the
transverse-field-driven order-disorder phase transition. Recently, the
multicriticality at $\eta=0$, i.e., the end-point singularity of the phase
boundary toward the $XX$-symmetric point, was explored in depth from the
quantum-information-theoretical viewpoint [7].

Meanwhile, its extention to the two-dimensional counterpart has been made [8,
9]. According to the large-$N$ theory [9], namely, for sufficiently large
internal symmetry group, the phase boundary should exhibit a reentrant (non-
monotonic) behavior. The reentrant behavior leads to a counterintuitive
picture such that the disorder phase is induced by the lower internal symmetry
group. On the contrary, as to the $XY$ (namely, easy-plane SU(2)) model, the
exact-diagonalization study [8, 9, 10] claimed that the phase boundary rises
up linearly (monotonically) in proximity to the multicritical point; the
similar phase diagram was reported in Ref. [11], where a topological index is
specified to each phase surrounding the multicritical point. Rather
technically, it has to be stressed that the exact-diagonalization method
allows direct access to the ground state (namely, infinite imaginary-time
system size), and hence, we do not have to care about the anisotropy between
the real-space and imaginary-time directions rendered by the dynamical
critical exponent $z=2(\neq 1)$ [12] at the multicritical point.

The aim of this paper is to reconcile the discrepancy between the large-$N$\-
[9] and easy-plane-SU$(2)$-based [8, 9, 10] results; namely, so far, only the
extremum cases have been considered. For that purpose, we extend the internal
symmetry group to the easy-plane SU(3) [13], and explore how the extention of
the internal symmetry group affects the multicriticality. Additionally, as a
probe to detect the phase transition [14, 15, 16, 17, 18, 19, 20, 21], we
resort to the fidelity [22, 23, 24, 25]

$F(H,H+\Delta H)=|\langle H|H+\Delta H\rangle|,$ (1)

where the vectors, $|H\rangle$ and $|H+\Delta H\rangle$, denote the ground
states with the proximate interaction parameters, $H$ and $H+\Delta H$,
respectively. The fidelity (1) is readily accessible via the exact-
diagonalization method, which yields the ground-state vector $|H\rangle$
explicitly. According to the elaborated exact-diagonalization study of the
two-dimensional $XXZ$ and Ising models [17], the fidelity-mediated analysis
admits a reliable estimate for the criticality, although the available system
size $N\leq 20$ is rather restricted.

To be specific, we present the Hamiltonian for the two-dimensional easy-plane
SU(3) magnet subjected to the transverse field

${\cal H}=-J\sum_{\langle
ij\rangle}[(1+\eta)(S^{x}_{i}S^{x}_{j}+S^{y}_{i}S^{y}_{j})+(1-\eta)(Q^{zx}_{i}Q^{zx}_{j}+Q^{yz}_{i}Q^{yz}_{j})]+H\sum_{i=1}^{N}Q^{z^{2}}_{i},$
(2)

with the $S=1$-spin operator ${\mathbf{S}}_{i}$ placed at each square-lattice
point $i=1,2,\dots,N$. Likewise, the quadrupolar moments at site $i$,
$Q^{x^{2}-y^{2}}_{i}=(S^{x}_{i})^{2}-(S^{y}_{i})^{2}$,
$Q^{z^{2}}_{i}=\sqrt{3}(S^{z}_{i})^{2}-\frac{2}{\sqrt{3}}$,
$Q^{xy}_{i}=S^{x}_{i}S^{y}_{i}+S^{y}_{i}S^{x}_{i}$,
$Q^{yz}_{i}=S^{y}_{i}S^{z}_{i}+S^{z}_{i}S^{y}_{i}$, and
$Q^{zx}_{i}=S^{z}_{i}S^{x}_{i}+S^{x}_{i}S^{z}_{i}$, are incorporated. These
eight operators constitute the SU$(3)$ algebra [26] just like the Gell-Mann
matrices. The summation $\sum_{\langle ij\rangle}$ runs over all possible
nearest-neighbor pairs $\langle ij\rangle$, and the coupling constant $J$ sets
the unit of energy, $J=1$, throughout this study. The parameters, $\eta$ and
$H$, denote the anisotropy of the internal symmetry and the transverse field,
respectively. The anisotropy $\eta\neq 0$ gives rise to the asymmetry between
the $S^{x,y}_{i}$ and $Q^{yz,zx}_{i}$ sectors. Irrespective of the anisotropy
$\eta$, the Hamiltonian commutes with the $z$-axis-rotation generator
$\sum_{i=1}^{N}S^{z}_{i}$, and hence, the Hamiltonian (2) describes the easy-
plane sector of the SU(3) magnet.

In Fig. 1, we present a schematic phase diagram of the easy-plane SU(3) magnet
(2) for the anisotropy $\eta$ and the transverse field $H$; the overall
character would resemble that of the $XY$ model [8]. The phase boundary
$H_{c}(\eta)$ separates the order ($H<H_{c}(\eta)$) and disorder
($H>H_{c}(\eta)$) phases. In the order phase, there appears the spontaneous
magnetization of the in-plane moment $(S^{x}_{i},S^{y}_{i})$
[$(Q^{zx}_{i},Q^{yz}_{i})$] for the $\eta>(<)0$ regime. Therefore, the phase
boundary $H=H_{c}(\eta)$ should belong to the three-dimensional (3D) $XY$
universality class except at the multicritical point $\eta=0$. At $\eta=0$,
the Hamiltonian (2) commutes with the transverse moment
$\sum_{i=1}^{N}Q^{z^{2}}_{i}$. Hence, as the transverse field $H$ increases,
the successive level crossings take place as to the ground-state energy up to
$H<H_{c}(0)(=8/\sqrt{3})$ [7], and above the threshold $H>H_{c}(0)$, the
transverse moment $\sum_{i=1}^{N}Q^{z^{2}}_{i}$ saturates eventually; see
Appendix. This transition mechanism is the same as that of the magnetization
plateau (saturation of magnetization), and the power-law singularities have
been investigated in depth [12]. Around the threshold $H=H_{c}(0)$, the
efficiency of the quantum Monte Carlo sampling suffers from the slowing-down
problem [27], and the aforementioned exact-diagonalization study [8]
circumvented this difficulty.

Then, there arises a problem how the phase boundary $H_{c}(\eta)$ terminates
at the multicritical point $\eta=0$; see Fig. 2. The exact-diagonalization
analysis of the $XY$ model indicates that the phase boundary $H_{c}(\eta)$
rises up linearly (monotonically) [8, 9] around the multicritical point
$\eta=0$. Actually, the power-law singularity of the phase boundary [28, 29]

$H_{c}(\eta)-H_{c}(0)\sim|\eta|^{1/\phi},$ (3)

is characterized by the crossover exponent $\phi\approx 1$ as for the $XY$
magnet [8, 9, 10]. On the one hand, the large-$N$ theory [9] suggests that the
phase boundary $H_{c}(\eta)$ shows a reentrant behavior. That is, the phase
boundary $H_{c}(\eta)$ exhibits a non-monotonic dependence on the anisotropy
$\eta$. So far, only the limiting cases such as the $XY$ [8, 9, 10] and
spherical [9] models have been considered, and no information has been
provided as to the multicriticality in between. As the internal symmetry group
is enlarged gradually from SU(2), the phase boundary may become curved
convexly, accompanying with the suppressed crossover exponent $\phi<1$. In
this paper, considering the SU(3) version of the easy-plane magnet, we
investigate the crossover exponent $\phi$ quantitatively by the agency of the
fidelity $F(H,H+\Delta H)$ (1).

The rest of this paper is organized as follows. In Sec. 2, we present the
numerical results. Prior to the detailed analysis of the multicriticality, we
investigate the case of the fully anisotropic limit $\eta=1$, aiming to
demonstrate the performance of our simulation scheme. In the last section, we
address the summary and discussions.

##  2 Numerical results

In this section, we present the numerical results for the easy-plane SU$(3)$
magnet subjected to the transverse field (2). We employed the exact-
diagonalization method for the finite-size cluster with $N\leq 5\times 5$
spins. The linear dimension of the cluster is given by $L=\sqrt{N}$, which
sets the fundamental length scale in the subsequent finite-size-scaling
analyses. As a probe to detect the phase transition, we utilized the fidelity
succeptibility [14, 15, 16, 17, 18, 19, 20, 21]

$\chi_{F}(H)=-\frac{1}{L^{2}}\partial^{2}_{\Delta H}F(H,H+\Delta H)|_{\Delta
H=0},$ (4)

with the fidelity $F(H,H+\Delta H)$ (1). As was demonstrated in Ref. [17] for
the two-dimensional $XXZ$ and Ising models under the transverse field with
$N\leq 20$ spins, the fidelity susceptibility $\chi_{F}$ (4) admits a reliable
estimate of the criticality, even though the available system size $N\leq 20$
is rather restricted.

In fairness, it has to be mentioned that the similar scheme was applied to the
one-dimensional $XY$ magnet under the transverse field [7]. In this pioneering
study [7]. the authors took a direct route toward the multicritical point,
$(\eta,H)\to(0,H_{c}(0))$. The $\chi_{F}$ data exhibit the intermittent peaks,
reflecting the successive level crossings along the ordinate axis $\eta=0$, as
shown in Fig. 1. In this paper, to avoid such a finite-size artifact, we took
a different route to the multicritical point, keeping the anisotropy $\eta$ to
a finite value. That is, based on the the crossover-scaling theory [28, 29],
the anisotropy $\eta$ is scaled properly, as the system size $L$ changes. As a
byproduct, we are able to estimate the crossover exponent $\phi$
quantitatively, which characterizes the power-law singularity of the phase
boundary; see Fig. 2. Before commencing detailed crossover-scaling analyses of
$\chi_{F}$, we devote ourselves to the fully anisotropic case $\eta=1$ so as
to examine the performance of our simulation scheme.

###  2.1  Phase transition point $H_{c}(1)$ at $\eta=1$: Fidelity-
susceptibility analysis

In this section, as a preliminary survey, setting the anisotropy parameter to
the fully anisotropic limit $\eta=1$, we investigate the order-disorder phase
transition of the easy-plane SU(3) magnet under the transverse field (2) via
the fidelity susceptibility $\chi_{F}$ (4). At this point $\eta=1$, the model
(2) reduces to the spin-$S=1$ $XY$ model with the single-ion anisotropy $D$,
for which a variety of preceding results are available [30, 31, 32, 33, 34,
35].

In Fig. 3, we present the fidelity susceptibility $\chi_{F}$ (4) for various
values of the transverse field $H$, and the system sizes, ($+$) $L=3$,
($\times$) $4$, and ($*$) $5$, with the fixed $\eta=1$. We see that the
fidelity susceptibility exhibits a pronounced signature for the order-disorder
phase transition around $H\approx 6$.

In Fig. 4, we present the approximate critical point $H_{c}^{*}(L)$ for
$1/L^{1/\nu}$ with $\eta=1$ fixed. Here, the approximate critical point
$H_{c}^{*}(L)$ denotes the location of the fidelity-susceptivity peak

$\partial_{H}\chi_{F}(L)|_{H=H_{c}^{*}(L)}=0,$ (5)

for each system size $L$. The power of the abscissa scale $1/\nu$ comes from
the scaling dimension of the parameter $H$ [20], and the correlation-length
critical exponent $\nu$ is set to the value of the 3D-$XY$ universality class,
$\nu=0.6717$ [36, 37]; the validity of this proposition is examined in the
next section. From Fig. 4, the least-squares fit to these data yields an
estimate

$H_{c}(\eta=1)=6.53(5),$ (6)

in the thermodynamic limit $L\to\infty$.

This is a good position to address an overview of the related studies. In
Table 1, we recollect a number of preceding results for the critical point
$H_{c}(1)$ at the fully anisotropic case $\eta=1$. As mentioned above, at
$\eta=1$, our model (2) reduces to the spin-$S=1$ $XY$ model with the single-
ion anisotropy $D$, and the $D$-based results are converted via the relation
$H=2D/\sqrt{3}$ so as to match our notation. As shown in table 1, the critical
point $H_{c}(1)$ was estimated as $H_{c}=6.300$ [30], $6.628$ [31, 32], and
$6.524(23)$ [33] by means of the bosonic mean-field approximation (BMFA),
self-consistent harmonic approximation (SCHA), and quantum Monte Carlo (QMC)
methods, respectively. (The BMFA estimate is read off from Fig. 1 of Ref.
[30].) As a probe to detect the phase transition, various quantifiers, such as
the spontaneous magnetization, energy gap, and correlation length, were
utilized in the respective studies. As indicated, we also incorporated the
Heisenberg-model [34, 35] result [33], because “the anisotropy does not have
large effects” [30] as to the critical point $H_{c}$. Our exact-
diagonalization (ED) result $H_{c}=6.53(5)$ [Eq. (6)] obtained via the
fidelity susceptibility $\chi_{F}$ appears to be accordant with these
preceding studies. Particularly, our result $H_{c}=6.53(5)$ [Eq. (6)] agrees
with the large-scale-QMC-simulation result $6.524(23)$ [33], validating the
$\chi_{F}$-mediated simulation scheme. Encouraged by this finding, we further
explore the critical behavior of $\chi_{F}$ in the next section.

###  2.2  Criticality of the fidelity susceptibility at $\eta=1$

In this section, we investigate the criticality of the fidelity susceptibility
$\chi_{F}$ (4) at the fully anisotropic limit, $\eta=1$. To begin with, we set
up the finite-size-scaling formula for the fidelity susceptibility [20]

$\chi_{F}=L^{\alpha_{F}/\nu}f\left((H-H_{c})L^{1/\nu}\right),$ (7)

with a scaling function $f$, and the fidelity-susceptibility critical exponent
$\alpha_{F}$; namely, the index $\alpha_{F}$ describes the singularity
$\chi_{F}\sim|H-H_{c}|^{-\alpha_{F}}$ at the critical point $H_{c}$. According
to Ref. [20], the critical exponent $\alpha_{F}$ satisfies the relation

$\alpha_{F}/\nu=\alpha/\nu+1,$ (8)

with the specific-heat critical index $\alpha$; namely, the index $\alpha$
describes the singularity of the specific heat as $C\sim|H-H_{c}|^{-\alpha}$.
We postulate that the the criticality belongs to the 3D-$XY$ universality
class [36, 37] for the anisotropic regime $\eta\neq 0$. Putting the existing
values [36], $\alpha=0.0151$ and $\nu=0.6717$, for the 3D-$XY$ universality
class into the scaling relation (8), we arrive at

$\alpha_{F}/\nu=0.9775.$ (9)

Notably, this index $\alpha_{F}/\nu=0.9775$ is larger than that of the
specific heat, $\alpha/\nu=-0.0225$ [36]; actually, the latter takes a
negative value. The scaling parameters, $\alpha_{F}/\nu$ and $\nu$, appearing
in the expression (7) are all fixed, and we are able to carry out the scaling
analysis of $\chi_{F}$ unambiguously.

In Fig. 5 we present the scaling plot,
$(H-H_{c})L^{1/\nu}$-$\chi_{F}L^{-\alpha_{F}/\nu}$, for various $H$ and system
sizes, ($+$) $L=3$, ($\times$) $4$, and ($*$) $5$, with the fixed $\eta=1$.
Here, the scaling parameters are set to $H_{c}=6.53$ [Eq. (6)], $\nu=0.6717$
[36], and $\alpha_{F}/\nu=0.9775$ [Eq. (9)]. The scaled data appear to
collapse into the scaling curve satisfactorily, confirming that the simulation
data already enter into the scaling regime. We stress that no ad hoc parameter
adjustment is undertaken in the present scaling analysis.

We address a number of remarks. First, the scaling plot, Fig. 5, indicates
that the phase transition belongs to the 3D-$XY$ universality class. So far,
as for the Heisenberg antiferromagnet with the single-ion anisotropy $D$, it
has been claimed that the $D$-driven singularity belongs to the 3D-$XY$
universality class [33, 35]. Our simulation result shows that the quantum $XY$
ferromagnet is also under the reign of the 3D-$XY$ universality class. Second,
from the scaling plot, Fig. 5, we see that corrections to finite-size scaling
are rather suppressed. Actually, as first noted by Ref. [17], the fidelity
susceptibility detects the underlying singularity clearly out of the
subdominant contributions. Last, a key ingredient is that $\chi_{F}$’s
singularity $\alpha_{F}/\nu=0.9775$ [Eq. (9)] is substantially larger than
that of the specific heat, $\alpha/\nu=-0.0225$ [36]; actually, the latter
takes a negative value as for the 3D-$XY$ universality. In this sense, the
former is appropriate as a quantifier for the criticality.

###  2.3  Crossover-scaling plot of the fidelity susceptibility around
$\eta=0$: Analysis of the crossover exponent $\phi$

We then turn to the analysis of the multicriticality at $\eta=0$. For that
purpose, we extend the above scaling formalism (7) to the crossover-scaling
formula [28, 29]

$\chi_{F}=L^{\dot{\alpha}_{F}/\dot{\nu}}g\left((H-H_{c}(\eta))L^{1/\dot{\nu}},\eta
L^{\phi/\dot{\nu}}\right),$ (10)

with yet another controllable parameter $\eta$, the accompanying crossover
exponent $\phi$, and a scaling function $g$. Here, the indices,
$\dot{\alpha}_{F}$ and $\dot{\nu}$, denote the fidelity-susceptibility and
correlation-length critical exponents, respectively, right at the
multicritical point $\eta=0$. The meaning of the second argument of the
crossover-scaling formula (10), $\eta L^{\phi/\dot{\nu}}$, is as follows. The
correlation-length critical exponent $\dot{\nu}$ describes the singularity of
the length scale $L$ as $L\sim|H-H_{c}(0)|^{-\dot{\nu}}$. This relation leads
to the physically convincing expression $\eta
L^{\phi/\dot{\nu}}\sim\eta/|H-H_{c}(0)|^{\phi}$. Now, it is apparent that the
crossover exponent $\phi$ describes the mutual relationship between $H$ and
the new entity $\eta$.

Before carrying out the crossover-scaling analyses, we fix the values of the
critical indices appearing in the formula (10). As mentioned in Introduction,
the phase transition along the ordinate axis $\eta=0$ is essentially the same
as that of the magnetization plateau, and the power-law singularities have
been investigated in considerable detail [12] as follows. The corelation-
length critical exponent was determined as $\dot{\nu}=1/2$ [12, 38, 39]. This
index immediately yields $\dot{\alpha}_{F}/\dot{\nu}=3$ via the scaling
relation $\dot{\alpha}_{F}/\dot{\nu}=\dot{\alpha}/\dot{\nu}+z$ [20]. Here, the
dynamical critical exponent takes the value $z=2$ [12], and the specific-heat
critical exponent $\dot{\alpha}=1/2$ comes from the idea that the field-
induced magnetization, $\sim\sqrt{H-H_{c}}$ [12], is regarded as the internal
energy. (Note that the first derivative of the ground-state energy corresponds
to the “internal energy” in the classical statistical mechanical context.) The
above consideration now completes the prerequisite for the crossover-scaling
analysis. The remaining index $\phi$ is to be adjusted so as to achieve an
alignment of the crossover-scaled $\chi_{F}$ data.

In Fig. 6, we present the crossover-scaling plot,
$(H-H_{c}(\eta))L^{1/\dot{\nu}}$-$\chi_{F}L^{-\dot{\alpha}_{F}/\dot{\nu}}$,
for various system sizes, ($+$) $L=3$, ($\times$) $4$, and ($*$) $5$, with
$\dot{\nu}=1/2$ and $\dot{\alpha}_{F}/\dot{\nu}=3$ determined above. Here, we
fixed the second argument of the scaling formula (10) to a constant value
$\eta L^{\phi/\dot{\nu}}=5$ with an optimal crossover exponent $\phi=0.8$, and
the parameter $H_{c}(\eta)$ was determined via the same scheme as that of Sec.
2.1. From Fig. 6, we see that the crossover-scaled data fall into the scaling
curve satisfactorily. Particularly, the $L=4$ ($\times$) and $5$ ($*$) data
are about to overlap each other. Such a feature suggests that the choice
$\phi=0.8$ is a feasible one.

Likewise, setting the crossover exponent to $\phi=0.7$, in Fig. 7, we present
the crossover-scaling plot, $(H-H_{c}(\eta))L^{2}$-$\chi_{F}L^{-3}$, for
various system sizes $L=3,4,5$; the symbols are the same as those of Fig. 6.
Here, the second argument of the crossover-scaling formula (10) is fixed to
$\eta L^{2\phi}=3.62$ with a proposition $\phi=0.7$. For such a small value of
$\phi=0.7$, the crossover-scaled data get scattered; particularly, the left-
side slope becomes dispersed, as compared to that of Fig. 6. Similarly, under
the setting $\phi=0.9$, in Fig. 8, we present the crossover-scaling plot,
$(H-H_{c}(\eta))L^{2}$-$\chi_{F}L^{-3}$, for various system sizes $L=3,4,5$;
the symbols are the same as those of Fig. 6. Here, the second argument of the
crossover-scaling formula is set to a constant value $\eta L^{2\phi}=6.9$ with
$\phi=0.9$. For such a large value of $\phi=0.9$, the crossover-scaled data
become scattered; actually, the hill-top data split up. Considering that the
above cases, Figs. 7 and 8, yield the lower and upper bounds, respectively,
for the estimate of $\phi$, we conclude that the crossover exponent lies
within

$\phi=0.8(1).$ (11)

Our result excludes the possibility that the phase boundary $H_{c}(\eta)$
rises up linearly $\phi\approx 1$, as observed for the $XY$ (namely, easy-
plane SU(2)) magnet [8, 9, 10]. Rather, as for the extended internal symmetry
group, the phase boundary is curved convexly, accompanying with a slightly
suppressed crossover exponent (11).

A number of remarks are in order. First, the underlying physics behind the
crossover-scaling plot, Fig. 6, differs from that of the fixed-$\eta$ scaling,
Fig. 5. Actually, the former scaling dimension $\dot{\alpha}_{F}/\dot{\nu}=3$
is much larger than that of the latter, $\alpha_{F}/\nu=0.9775$ [Eq. (9)].
Therefore, the index $\phi$ has to be adjusted carefully, and the collapse of
the crossover-scaled data points, Fig. 6, is by no means coincidental. In this
sense, the present crossover-scaling analysis captures the characteristics of
the multicritical behavior. Second, because the fidelity-susceptibility
approach does not rely on any presumptions as to the order parameter involved,
it detects both singularities, $\dot{\alpha}/\dot{\nu}=3$ and
$\alpha/\nu=0.9775$ [Eq. (9)], in a unified manner. Moreover, as mentioned in
Sec. 2.2, the fidelity susceptibility is less affected by corrections to
scaling [17], and it picks up the underlying singularity out of the
subdominant contributions. Last, as explained in Introduction, so far, only
the extremum cases, namely, the large-$N$ [9] and easy-plane-SU(2) [8, 9, 10]
symmetry groups, have been considered, and the multicritical behavior in
between remains unclear. Our result $\phi=0.8(1)(<1)$ [Eq. (11)] indicates
that the phase boundary $H_{c}(\eta)$ is curved convexly, as the internal
symmetry group is extended to the easy-plane SU(3). Because the reentrant
behavior occurs only in $1<d\leq 2.065\dots$ dimensions even for the large-$N$
case [9], it is reasonable that the slightly enlarged SU(3) group does not
lead to the reentrant behavior immediately in $d=2$ dimensions.

##  3 Summary and discussions

The two-dimensional easy-plane SU(3) magnet under the transverse field (2) was
investigated with the exact-diagonalization method. Because the method allows
direct access to the ground state, we do not have to care about the anisotropy
between the real-space and imaginary-time directions rendered by the dynamical
critical exponent $z=2(\neq 1)$ [12] at the multicritical point $\eta=0$. Our
main concern is to reconcile the discrepancy between the large-$N$ theory [9]
and the exact-diagonalization analysis of the $XY$ magnet [8, 9, 10] as to the
multicriticality (Fig. 2). For that purpose, we consider the SU(3) version of
the easy-plane magnet (2), and performed the exact-diagonalization simulation
by the agency of the fidelity susceptibility $\chi_{F}$ (4). The fidelity
susceptibility has an advantage in that it detects the phase transition
sensitively, as compared to that of the specific heat [17, 20]. As a
preliminary survey, setting the anisotropy parameter to $\eta=1$, we estimate
the critical point via $\chi_{F}$ as $H_{c}(1)=6.53(5)$ [Eq. (6)]. This result
is comparable to the preceding results, $6.300$ [30], $6.628$ [31, 32], and
$6.524(23)$ [33], estimated with the BMFA, SCHA, and QMC methods,
respectively. Particularly, our result $H_{c}(1)=6.53(5)$ agrees with the
large-scale-QMC-simulation result $6.524(23)$ [33], validating the
$\chi_{F}$-mediated simulation scheme. Thereby, we cast the $\chi_{F}$ data
into the crossover-scaling formula (10) with the anisotropy $\eta$ scaled
properly. Adjusting the crossover exponent $\phi$ carefully, we attain an
alignment of the crossover-scaled data points for $\phi=0.8(1)$ [Eq. (11)].
This result $\phi<1$ indicates that the phase boundary $H_{c}(\eta)$ gets
curved convexly around the multicritical point $\eta=0$, as the internal
symmetry group is extended to SU(3). Actually, only in $1<d\leq 2.065\dots$
dimensions, the reentrant behavior is realized even for the large-$N$ case
[9]. Hence, it is reasonable that the slightly enlarged internal symmetry does
not lead to the reentrant behavior immediately.

We conjecture that the phase boundary may exhibit a quadratic curvature,
$\phi=0.5$, eventually for an extremely large internal symmetry group, and
above this threshold, the reentrant behavior should set in. It is thus
tempting to consider the fractional-dimensional ($1<d<2$) system realized
effectively by the power-law-decaying interactions [40]. In such a fractional-
dimensional system, the reentrant behavior may come out even for the $XY$
model. This problem is left for the future study.

This work was supported by a Grant-in-Aid for Scientific Research (C) from
Japan Society for the Promotion of Science (Grant No. 20K03767).

Table 1:  A number of preceding results for the critical point $H_{c}(1)$ at $\eta=1$ are recollected. At $\eta=1$, the model (2) reduces to the spin-$S=1$ $XY$ model with the single-ion anisotropy $D$. These $D$-based results are converted so as to match our notation via the relation $H=2D/\sqrt{3}$. So far, a variety of technieus, such as the bosonic mean-field approximation (BMFA) [30], self-consistent harmonic approximation (SCHA) [31, 32], and quantum Monte Carlo (QMC) [33] methods, have been employed, and these results are comparable to the present exact diagonalization (ED) result. (The BMFA estimate is read off from Fig. 1 of Ref. [30].) In the respective studies, as a probe to detect the phase transition, various quantifiers, such as the spontaneous magnetization, energy gap, correlation length and fidelity susceptibility, were utilized. Because “the anisotropy does not have large effects” [30], the large-scale QMC result [33] for the Heisenberg model [34, 35] is shown as well.  Method | Quantifier | Model | $H_{c}(\eta=1)$
---|---|---|---
BMFA [30] | spontaneous magnetization | $XY$ | $6.300$
SCHA [31, 32] | energy gap | $XY$ | $6.628$
QMC [33] | correlation length | Heisenberg | $6.524(23)$
ED (this work) | fidelity susceptibility | $XY$ | $6.53(5)$
Figure 1:  A schematic phase diagram of the easy-plane SU$(3)$ magnet (2) for
the anisotropy $\eta$ and transverse field $H$ is presented. The phase
boundary $H_{c}(\eta)$ terminates at the multicritical point
$H=H_{c}(\eta=0)=8/\sqrt{3}$ at $\eta=0$, and this end-point singularity is
our concern. For small $H<H_{c}(\eta)$, the in-plane order,
$(S^{x}_{i},S^{y}_{i})$ [$(Q^{yz}_{i},Q^{zx}_{i})$] develops in the
$\eta>(<)0$ side, whereas for large $H>H_{c}(\eta)$, the disorder phase
extends. At the isotropic point $\eta=0$ in between, the transverse moment
$\sum_{i=1}^{N}Q^{z^{2}}_{i}$ commutes with the Hamiltonian, and the ground-
state level crossing occurs successively [7] up to $H<H_{c}(0)$. Above this
threshold $H>H_{c}(0)$, the magnetization plateau [12] (saturation of the
moment $\sum_{i=1}^{N}Q^{z^{2}}_{i}$) sets in.  Figure 2:  The
multicriticality (end-point singularity) of the phase boundary $H_{c}(\eta)$
at $\eta=0$ is characterized by the crossover exponent $\phi$ [28, 29] such as
$H_{c}(\eta)-H_{c}(0)\sim|\eta|^{1/\phi}$ (3). As for the $XY$ model (namely,
the easy-plane SU$(2)$ magnet), the exact-diagonalization simulation suggests
that the phase boundary rises up linearly with $\phi\approx 1$ [8, 9, 10]. On
the one hand, the large-$N$ analysis admits the reentrant (non-monotonic)
behavior [9]. As indicated, the reentrant behavior leads to a counterintuitive
picture such that the disorder phase is induced by the lower internal symmetry
group $\eta\neq 0$ around the multicritical point. It is anticipated that for
the SU$(3)$ case, the phase boundary is curved convexly with a slightly
suppressed crossover exponent $\phi<1$.  Figure 3:  The fidelity
susceptibility $\chi_{F}$ (4) is plotted for various values of the transverse
field $H$ and the system sizes, ($+$) $L=3$, ($\times$) $4$, and ($*$) $5$,
with the fixed anisotropy parameter $\eta=1$. The fidelity susceptibility
indicates a notable peak around the critical point $H\approx 6$.  Figure 4:
The approximate critical point $H_{c}^{*}(L)$ (5) is plotted for $1/L^{1/\nu}$
with the fixed anisotropy parameter $\eta=1$. Here, the correlation-length
critical exponent is set to the 3D-$XY$-universality value $\nu=0.6717$ [36,
37]. The least-squares fit to these data yields an estimate $H_{c}=6.53(5)$ in
the thermodynamic limit $L\to\infty$.  Figure 5:  The scaling plot,
$(H-H_{c})L^{1/\nu}$-$\chi_{F}L^{-\alpha_{F}/\nu}$, is presented for various
system sizes, ($+$) $L=3$, ($\times$) $4$, and ($*$) $5$, with the fixed
$\eta=1$; see the scaling formula (7). Here, the scaling parameters are set to
$H_{c}=6.53$ [Eq. (6)], $\nu=0.6717$ [36, 37] and $\alpha_{F}/\nu=0.9775$ [Eq.
(9)]. The scaled data collapse into the scaling curve satisfactorily,
suggesting that the criticality belongs to the 3D-$XY$ universality class.
Figure 6:  The crossover-scaling plot,
$(H-H_{c}(\eta))L^{1/\dot{\nu}}$-$\chi_{F}L^{-\dot{\alpha_{F}}/\dot{\nu}}$, is
presented for various system sizes, ($+$) $L=3$, ($\times$) $4$, and ($*$)
$5$, with $\dot{\nu}=1/2$ [12] and $\dot{\alpha}_{F}/\dot{\nu}=3$; see text
for details. Here, the second argument of the crossover-scaling formula (10)
is fixed to $\eta L^{\phi/\dot{\nu}}=5$ with the crossover exponent
$\phi=0.8$. The crossover-scaled data fall into the scaling curve
satisfactorily; particularly, the $L=4$ ($\times$) and $5$ ($*$) data are
about to overlap each other under the setting, $\phi=0.8$.  Figure 7:  The
crossover-scaling plot,
$(H-H_{c}(\eta))L^{1/\dot{\nu}}$-$\chi_{F}L^{-\dot{\alpha_{F}}/\dot{\nu}}$, is
presented for various system sizes, ($+$) $L=3$, ($\times$) $4$, and ($*$)
$5$, with $\dot{\nu}=1/2$ and $\dot{\alpha}_{F}/\dot{\nu}=3$. Here, the second
argument of the crossover-scaling formula (10) is fixed to $\eta
L^{\phi/\dot{\nu}}=3.62$ with the crossover exponent $\phi=0.7$. The left-side
slope gets scattered under such a small value of $\phi=0.7$.  Figure 8:  The
crossover-scaling plot,
$(H-H_{c}(\eta))L^{1/\dot{\nu}}$-$\chi_{F}L^{-\dot{\alpha_{F}}/\dot{\nu}}$, is
presented for various system sizes, ($+$) $L=3$, ($\times$) $4$, and ($*$)
$5$, with $\dot{\nu}=1/2$ and $\dot{\alpha}_{F}/\dot{\nu}=3$. Here, the second
argument of the crossover-scaling formula (10) is fixed to $\eta
L^{\phi/\dot{\nu}}=6.9$ with the crossover exponent $\phi=0.9$. The hill-top
data split up for such a large value of $\phi=0.9$.

##  Appendix A  Transition point $H_{c}=8/\sqrt{3}$ at $\eta=0$

At the isotropic point $\eta=0$, the location of the multicritical point
$H_{c}(0)=8/\sqrt{3}$ is obtained analytically as follows. For sufficiently
large $H$, the ground state is given by the direct product
$\otimes_{i}^{N}|0\rangle_{i}$ of the local base $|0\rangle_{i}$, which
satisfies $S^{z}_{i}|m\rangle_{i}=m|m\rangle_{i}$ ($m=-1,0,1$) at each site
$i$. Due to the strict selection rule at $\eta=0$, magnons’ pair creation is
prohibited. Thus, the single magnon at site $j$,
$|j\rangle=|1\rangle_{j}\otimes(\otimes_{i\neq j}^{N}|0\rangle_{i})$,
propagates coherently through the transfer amplitude $-2J$ over the nearest
neighbors, obeying the dispersion relation $-4J(\cos k_{x}+\cos
k_{y})+\sqrt{3}H$ (${\mathbf{k}}$: wave number) above the ground state.
Therefore, at $H_{c}(\eta=0)=8/\sqrt{3}$, the band gap closes in a way
reminiscent of the metal-insulator transition. This transition mechanism is
precisely the same as that of the magnetization plateau [12] (saturation of
the magnetization). A notable point is that the dynamical critical exponent
takes $z=2$, because the quadratic band bottom touches the ground-state energy
level. In other words, the symmetry between the real-space and imaginary-time
directions is violated by $z\neq 1$. It is a benefit of the exact-
diagonalization method that the method allows direct access to the ground
state, for which the imaginary-time system size is infinite.

## References

## References

  * [1] J. Maziero, H. C. Guzman, L. C. Céleri, M. S. Sarandy, and R. M. Serra, Phys. Rev. A 82 (2010) 012106.
  * [2] Z.-Y. Sun, Y.-Y. Wu, J. Xu, H.-L. Huang, B.-F. Zhan, B. Wang, and C.-B. Duanpra, Phys. Rev. A 89 (2014) 022101.
  * [3] G. Karpat, B. Çakmak, and F. F. Fanchini, Phys. Rev. B 90 (2014) 104431.
  * [4] Q. Luo, J. Zhao, and X. Wang, Phys. Rev. E 98 (2018) 022106.
  * [5] A. Steane, Rep. Prog. Phys. 61 (1998) 117.
  * [6] C.H. Bennett and D.P. DiVincenzo, Nature 404 (2000) 247.
  * [7] V. Mukherjee, A. Polkovnikov, and A. Dutta, Phys. Rev. B 83 (2011) 075118.
  * [8] M. Henkel, J. Phys. A: Mathematical and Theoretical 17 (1984) L795.
  * [9] S. Wald and M. Henkel, J. Stat. Mech.: Theory and Experiment (2015) P07006.
  * [10] Y. Nishiyama, Eur. Phys. J. B 92 (2019) 167.
  * [11] S. Jalal, R. Khare, and S. Lal, arXiv:1610.09845.
  * [12] V. Zapf, M. Jaime, and C. D. Batista, Rev. Mod. Phys. 86 (2014) 563.
  * [13] J. D’Emidio and R. K. Kaul, Phys. Rev. B 93 (2016) 054406.
  * [14] H. T. Quan, Z. Song, X. F. Liu, P. Zanardi, and C. P. Sun, Phys. Rev. Lett. 96 (2006) 140604.
  * [15] P. Zanardi and N. Paunković, Phys. Rev. E 74 (2006) 031123.
  * [16] H.-Q. Zhou, and J. P. Barjaktarevic̃, J. Phys. A: Math. Theor. 41 (2008) 412001.
  * [17] W.-C. Yu, H.-M. Kwok, J. Cao and S.-J. Gu, Phys. Rev. E 80 (2009) 021108.
  * [18] W.-L. You and Y.-L. Dong, Phys. Rev. B 84 (2011) 174426.
  * [19] D. Rossini and E. Vicari, Phys. Rev. E 98 (2018) 062137.
  * [20] A. F. Albuquerque, F. Alet, C. Sire, and S. Capponi, Phys. Rev. B 81 (2010) 064418.
  * [21] D. Schwandt, F. Alet, and S. Capponi, Phys. Rev. Lett. 103 (2009) 170501.
  * [22] A. Uhlmann, Rep. Math. Phys. 9 (1976) 273.
  * [23] R. Jozsa, J. Mod. Opt. 41 (1994) 2315.
  * [24] A. Peres, Phys. Rev. A 30 (1984) 1610.
  * [25] T. Gorin, T. Prosen, T. H. Seligman, and M. Žnidarič, Phys. Rep. 435 (2006) 33.
  * [26] A. Läuchli, F. Mila, and K. Penc, Phys. Rev. Lett. 97 (2006) 087205.
  * [27] V. Z. Kashurnikov, N. V. Prokof’ev, B. V. Scistunov, and M. Troyer, Phys. Rev. B 59 (1999) 1162.
  * [28] E.K. Riedel and F. Wegner, Z. Phys. 225 (1969) 195.
  * [29] P. Pfeuty, D. Jasnow, and M. E. Fisher, Phys. Rev. B 10 (1974) 2088.
  * [30] H.-T. Wang and Y. Wang, Phys. Rev. B 71 (2005) 104429.
  * [31] A.S.T. Pires, J. Mag. Mag. Mat. 323 (2011) 1977.
  * [32] A.R. Moura, J. Mag. Mag. Mat. 369 (2014) 62.
  * [33] T. Roscilde and S. Haas, Phys. Rev. Lett. 99 (2007) 047205.
  * [34] C. J. Hamer, O. Rojas, and J. Oitmaa, Phys. Rev. B 81 (2010) 214424.
  * [35] Z. Zhang, K. Wierschem, I. Yap, Y. Kato, C. D. Batista, and P. Sengupta, Phys. Rev. B 87 (2013) 174405.
  * [36] M. Campostrini, M. Hasenbusch, A. Pelissetto, and E. Vicari, Phys. Rev. B 74, 144506 (2006).
  * [37] E. Burovski, J. Machta, N. Prokof’ev, and B. Svistunov, Phys. Rev. B 74, 132502 (2006).
  * [38] M. Adamski, J. Jȩdrzejewski, and T. Krokhmalskii, arXiv:1502.05268.
  * [39] C. Hoeger, G. V. Gehlen, and V. Rittenberg, J. Phys. A: Math. Gen. 18 (1985) 1813.
  * [40] N. Defenu, A. Trombettoni, and S. Ruffo, Phys. Rev. B 96 (2017) 104432.


# Charged-particle multiplicity dependence of charm-baryon-to-meson ratio in
high-energy proton-proton collisions

Yu Chen  Min He  Department of Applied Physics, Nanjing University of Science
and Technology, Nanjing 210094, China

###### Abstract

We propose that the charged-particle multiplicity dependence of the charm-
baryon-to-meson ratio observed in high-energy $pp$ collisions can be explained
by canonical treatment of quantum charges in the statistical hadronization
model (SHM). Taking the full particle listings of PDG complemented by
additional charm-baryon states from relativistic quark model predictions, we
evaluate the canonical partition function and the charm-hadron chemical
factors that measure the canonical suppression arising from the requirement of
strict conservation of quantum charges. We demonstrate that, while charm
number conservation induces common suppression on the production of both
charm-baryons and -mesons, baryon (strangeness) number conservation causes
further suppression on charm-baryons (charm-strange mesons) relative to
nonstrange charm-mesons, thereby resulting in a decreasing $\Lambda_{c}/D^{0}$
($D_{s}/D^{0}$) ratio toward smaller multiplicity events. The charm-hadron
thermal densities thus computed are then used as pertinent weights to perform
charm-quark fragmentation simulations yielding $p_{T}$-dependent
$\Lambda_{c}/D^{0}$ and $D_{s}/D^{0}$ ratios at varying multiplicities in fair
agreement with ALICE measurements.

###### keywords:

Heavy flavor production , Charm baryons , Canonical ensemble , Proton-proton
collisions

##  1 Introduction

Factorization provides an approach to compute the production of heavy-flavor
(HF) hadrons in high-energy hadronic collisions within the framework of
perturbative QCD. This approach involves the convolution of the parton
distribution functions (PDFs) of the colliding hadrons, the hard-scattering
cross sections at the partonic level and the fragmentation functions (FFs).
While the PDFs are considered to be universal and usually taken from global
fits of various sets of empirical data [1], and the partonic cross sections
are calculated as perturbative series in powers of the strong-coupling
constant, the FFs modelling the nonperturbative hadronization processes
represent a major source of theoretical uncertainty.

FFs parameterize the percentage of the quark energy transferred to the
produced hadron and the fraction (weight) of a heavy quark hadronizing into a
particular hadron species. Under the assumption that they are universal across
different collision systems, FFs are usually tuned against electron-positron
($e^{+}e^{-}$) collision data [2]. However, this assumption of universality
has been challenged by the measurement of production of heavy baryons, which
is in particular sensitive to the fragmentation mechanism of heavy quarks
(HQs) [3]. For example, in the case of bottom quarks, the
$b\rightarrow\Lambda_{b}^{0}$ fragmentation fraction was reported not to be
compatible with that deduced from LEP data [4]. In the case of charm quarks,
fragmentation calculations based on fit of the $e^{+}e^{-}$ knowledge in the
$k_{T}$-factorization approach [5] or the general-mass variable-flavor number
scheme [6] reported a substantial deficiency in the $\Lambda_{c}$ production
cross section as compared to the measurement by the ALICE collaboration at
mid-rapidity in $\sqrt{s}=7$ TeV $pp$ collisions [3]; and therefore the the
pertinent $\Lambda_{c}/D^{0}\simeq 0.54$ ratio, which is also regarded as a
sensitive measure of novel recombination hadronization mechanisms in heavy-ion
collisions [7, 8, 9], was much underestiamted [5, 6]. String fragmentation
invoking color reconnection mechanisms accountes for multi-parton interactions
in high-energy $pp$ collisions and turns out to enhance the baryon production
in both strange and charm sector [10].

Another way of determining the fragmentation fractions of a heavy quark
hadronizing into various kinds of HF hadrons is given by their thermal weights
at a hadronization “temperature” $T_{H}$ in the statistical hadronization
model (SHM) [11, 12], given that relative chemical equilibrium among them may
be reached in the quark-rich environment such as high-energy $pp$ collisions.
It was shown [13] that by incorporating a large set of hitherto unobserved
charm-baryon states as inspired by relativistic quark model (RQM) [14]
computations into the grand-canonical SHM [15, 16], the large
$\Lambda_{c}/D^{0}\simeq 0.54$ reported by ALICE [3] could be well accounted
for, simply augmented by the enhanced feedowns from those “missing” charm-
baryons in particular in the low momentum regime.

While the aforementioned $\Lambda_{c}/D^{0}\simeq 0.54$ was measured in
minimum bias $pp$ collisions, measurement of the dependence of HF production
in such collisions on the charged particle multiplicity ($dN_{\rm ch}/d\eta$)
may provide insight into the multi-parton interactions as well as the
interplay between the hard and soft mechanisms in particle production [17, 18,
19, 20, 21]. Indeed, study of $pp$ collisions at varying $dN_{\rm ch}/d\eta$
is of considerable interest as phenomena typically indicative of the formation
of a deconfined quark-gluon medium have been observed in high-multiplicity
$pp$ collisions, such as the long-range two-particle azimuthal correlations on
the near side [22] and the enhanced production of multi-strange hadrons [23],
opening up an angle for examining the final state effects known from heavy-ion
collisions. As far as the $\Lambda_{c}/D^{0}$ is concerned, its dependence on
$dN_{\rm ch}/d\eta$ [24] would reflect the modification of the fragmentation
mechanisms of charm quarks as the surrounding environment changes.

In the present work, we generalize the grand-canonical SHM investigation [13]
of charm-baryon production to the case of canonical SHM, and explore in how
far the observed multiplicity dependence of $\Lambda_{c}/D^{0}$ in high-energy
$pp$ collisions [24] can be explained by the additional canonical suppression
of the production of charm-baryons relative to charm-mesons. Particles
entering the evaluation of the canonical partition function are taken from the
full listings of PDG [25] (excluding bottom particles), with charm-baryons
still augmented by the additional states predicted by the relativistic quark
model (RQM) computations [13, 14]. The existence of those “missing” charm-
baryon states is supported by the direct measurement of charm-baryon
spectroscopy on lattice [26] as well on the analysis of the partial pressure
of open-charm states and charm-quark susceptibilities in thermal lattice QCD
[27]. We calculate the chemical factors of each charm-hadron as a measure of
canonical suppression relative to the grand-canonical (large volume) limit. We
demonstrate that, while strict charm number conservation as imposed by
canonical treatment of SHM induces common suppression on both charm-baryons
and -mesons, requirement of strict baryon number conservation causes further
suppression on charm-baryons, resulting in a decreasing $\Lambda_{c}/D^{0}$
toward smaller system size (or multiplicity $dN_{\rm ch}/d\eta$). Similarly,
strict strangeness conservation causes further suppression on charm-strange
mesons and leads to a multiplicity-dependent $D_{s}/D^{0}$ ratio. Using thus
computed thermal densities of each charm-hadron as pertinent weights, we
conduct a simulation of the charm quark fragmentation within the FONLL scheme
[28, 29] and compute the transverse momentum ($p_{T}$) dependent
$\Lambda_{c}/D^{0}$ and $D_{s}/D^{0}$ at different $dN_{\rm ch}/d\eta$.

##  2 Canonical ensemble SHM

###  2.1 Canonical partition function

When applied to heavy-ion collisions, SHM is usually implemented in a grand-
canonical ensemble (GCE) of an ideal hadron resonance gas (HRG), in which
quantum charges, such as electric charge $Q$, baryon number $B$, strangeness
$S$ and charm number $C$ are conserved on average and regulated by the
corresponding chemical potentials
$\vec{\mu}=(\mu_{Q},\mu_{B},\mu_{S},\mu_{C})$. For hadrons whose mass is much
larger than the typical hadronization temperature $T_{H}$ and thus the quantum
statistics effect can be safely neglected, the primary mean number of the
$j$-th hadron produced from this GCE-SHM is then given by

$\langle
N_{j}\rangle^{GCE}=\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}e^{\vec{\mu}\cdot\vec{q_{j}}/T_{H}},$
(1)

where $\vec{q}_{j}=(Q_{j},B_{j},S_{j},C_{j})$ is the corresponding quantum
charges and $z_{j}$ is the one-particle partition function of the $j$-th
hadron of mass $m_{j}$ and spin $J_{j}$

$\displaystyle
z_{j}=(2J_{j}+1)\frac{VT_{H}}{2\pi^{2}}m_{j}^{2}K_{2}(\frac{m_{j}}{T_{H}}),$
(2)

which is just the chemical-equilibrium multiplicity of the $j$-th hadron in
the grand-canonical limit under the Boltzmann approximation, with $K_{2}$
being the modified Bessel function of the second order. In Eq. (1),
$\gamma_{s}$ and $\gamma_{c}$ are fugacities to account for the deviation of
strange and charm production from chemical equilibrium, respectively, for
hadrons containing $N_{sj}$ ($N_{cj}$) strange (charm) valence quarks and
antiquarks. As the hadronization temperature $T_{H}$ is usually assumed to be
a universal parameter in the GCE implementation of SHM, the hadron yield
ratios would be the same in all collision centralities.

However, it has long been pointed out [30, 31] that requirement of exact
conservation of quantum charges as implemented in the canonical ensemble (CE)
treatment of SHM becomes increasingly important for sufficiently small
reaction volumes like those in elementary $pp$ and $e^{+}e^{-}$ collisions;
thereby a system size (or charged particle multiplicity $dN_{ch}/d\eta$)
dependence of hadron production would follow from the canonical corrections
[32]. For a collision system having conserved quantum charges
$\vec{Q}=(Q,B,S,C)$, the CE-SHM partition function at a given hadronization
temperature $T_{H}$ reads [33, 34, 35]

$\displaystyle Z(\vec{Q})=\frac{1}{(2\pi)^{4}}\int_{0}^{2\pi}d^{4}\phi
e^{i\vec{Q}\cdot\vec{\phi}}{\rm
exp}[\sum_{j}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}e^{-i\vec{q}_{j}\cdot\vec{\phi}}z_{j}],$
(3)

where the summation $\sum_{j}$ is over all hadrons and the volume in $z_{j}$
now should be understood as the correlation volume in which the quantum
charges are strictly conserved (it will be still denoted as $V$ throughout the
work, except in Sec. 4 it is replaced with $V_{c}$). Here, again the Boltzmann
approximation is assumed, which causes an error of few percent for pions (for
typical hadronization temperature $T_{H}=160-170$ MeV) when resonance decays
are taken into account, and much less for all other hadrons [34]. Note that,
unlike the case of GCE-SHM, the logarithm of the canonical partition function
${\rm log}Z(\vec{Q})$ does not scale linearly with the volume.

The primary mean number of the $j$-th hadron produced from this CE-SHM is then
given by [34, 35]

$\displaystyle\langle
N_{j}\rangle^{CE}=\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}\frac{Z(\vec{Q}-\vec{q}_{j})}{Z(\vec{Q})}.$
(4)

Comparing Eq. (1) and Eq. (4) reveals that the chemical potential term in the
GCE-SHM associated with conservation of quantum charges on average is now
replaced with the chemical factor $Z(\vec{Q}-\vec{q}_{j})/Z(\vec{Q})$ in the
CE-SHM arising from the requirement of exact conservation of quantum charges.
It will be shown (see Sec. 3) that the chemical factors for charged hadrons
(i.e., $\vec{q}_{j}\neq 0$) are always less than unity in a completely neutral
system (i.e., $\vec{Q}=0$ and thus $\vec{\mu}=0$) and reach their grand
canonical value of unity at asymptotically large volumes, therefore
characterizing the canonical suppression of the production of charged hadrons
at finite volume.

###  2.2 Evaluation of the canonical partition function

To evaluate the canonical partition function, one divides hadrons to be summed
over in the exponential of Eq.(3) into three categories [36]. The first
category consists of completely neutral mesons ($M^{0}$) with
$\vec{q}_{j}=(Q_{j},B_{j},S_{j},C_{j})=(0,0,0,0)$, e.g., $\pi^{0}$, $\eta$,
$\rho^{0}$, $a_{0}(980)^{0}$, $\phi$, $J/\psi$. Clearly their contribution to
the canonical partition function can be factorized out, yielding a
multiplicative term

$Z_{0}={\rm
exp}[\sum_{j=M^{0}}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}]=\prod_{j=M^{0}}{\rm
exp}(\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}).$ (5)

The second category refers to charged mesons, including positively charged
mesons ($M^{+}$) with at least one of $Q_{j}=+1$, $S_{j}=+1$ or $C_{j}=+1$
(but definitely $B_{j}=0$), e.g., $\pi^{+}$, $\rho^{+}$, $K^{+}$, $K^{0}$,
$K^{*+}$, $K^{*0}$, $D^{0}$, $D^{+}$, $D^{*0}$, $D^{*+}$, $D_{s}^{+}$,
$D_{s}^{*+}$, and their anti-mesons ($M^{-}$), e.g., $\pi^{-}$, $\rho^{-}$,
$K^{-}$, $\bar{K}^{0}$, $K^{*-}$, $\bar{K}^{*0}$, $\bar{D}^{0}$, $D^{-}$,
$\bar{D}^{*0}$, $D^{*-}$, $D_{s}^{-}$, $D_{s}^{*-}$. Noting that $M^{+}$’s and
$M^{-}$’s always have opposite charges but the same $N_{sj}$ and $N_{cj}$, the
contribution of the charged mesons to the exponential in Eq. (3)

$\displaystyle{\rm
exp}[\sum_{j=M^{+},M^{-}}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}e^{-i(Q_{j}\phi_{Q}+S_{j}\phi_{S}+C_{j}\phi_{C})}]$
$\displaystyle={\rm
exp}[\sum_{j=M^{+}}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}e^{-i(Q_{j}\phi_{Q}+S_{j}\phi_{S}+C_{j}\phi_{C})}$
$\displaystyle+\sum_{j=M^{+}}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}e^{i(Q_{j}\phi_{Q}+S_{j}\phi_{S}+C_{j}\phi_{C})}]$
$\displaystyle={\rm
exp}[2\sum_{j=M^{+}}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}{\rm
cos}(Q_{j}\phi_{Q}+S_{j}\phi_{S}+C_{j}\phi_{C})].$ (6)

The third category comprises baryons ($B$) with $B_{j}=+1$ and antibaryons
($\bar{B}$) with $B_{j}=-1$. Obviously $B$’s and $\bar{B}$’s carry opposite
values of all other charges ($Q_{j}$, $S_{j}$ and $C_{j}$), so that their
contribution to the exponential in Eq. (3)

$\displaystyle{\rm
exp}[\sum_{j=B,{\bar{B}}}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}e^{-i(Q_{j}\phi_{Q}+B_{j}\phi_{B}+S_{j}\phi_{S}+C_{j}\phi_{C})}]$
$\displaystyle={\rm
exp}[e^{-i\phi_{B}}\sum_{j=B}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}e^{-i(Q_{j}\phi_{Q}+S_{j}\phi_{S}+C_{j}\phi_{C})}$
$\displaystyle+e^{i\phi_{B}}\sum_{j=B}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}e^{i(Q_{j}\phi_{Q}+S_{j}\phi_{S}+C_{j}\phi_{C})}$
$\displaystyle={\rm exp}[e^{-i\phi_{B}}\omega+e^{i\phi_{B}}\omega^{*}]$
$\displaystyle={\rm exp}[2|\omega|{\rm cos}(\phi_{B}-{\rm arg}\omega)],$ (7)

where

$\displaystyle\omega$
$\displaystyle=\omega(\phi_{Q},\phi_{S},\phi_{C})=|\omega|e^{i{\rm
arg}\omega}$
$\displaystyle=\sum_{j=B}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}e^{-i(Q_{j}\phi_{Q}+S_{j}\phi_{S}+C_{j}\phi_{C})}.$
(8)

Now combining contributions from all these three categories gives rise to the
canonical partition function

$\displaystyle
Z(\vec{Q})=\frac{Z_{0}}{(2\pi)^{3}}\int_{0}^{2\pi}d^{3}\phi{\cos}(Q\phi_{Q}+S\phi_{S}+C\phi_{C}+B{\rm
arg}\omega)$ $\displaystyle\times I_{B}(2|\omega|)\times
e^{2\sum_{j=M^{+}}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}{\rm
cos}(Q_{j}\phi_{Q}+S_{j}\phi_{S}+C_{j}\phi_{C})},$ (9)

where $I_{B}(2|\omega|)=\frac{1}{2\pi}\int_{0}^{2\pi}d\phi_{B}e^{2|\omega|{\rm
cos}\phi_{B}}{\rm cos}(B\phi_{B})$ is the modified Bessel function of the
first kind and order $B$.

Figure 1: The behavior of the canonical partition function at $V=20~{}{\rm
fm^{3}}$, $T_{H}=170$ MeV and $\gamma_{s}=\gamma_{c}=1$, as a function of
integer electric charge ($Q$), baryon number ($B$), strangeness ($S$) and
charm number ($C$), respectively, with remaining quantum charges set to zero.

Taking the full particle listings in PDG [25] with charm-baryons augmented by
additional states from relativistic quark model (RQM) computations [13, 14],
the canonical partition function is evaluated and the behavior as a function
of increasing $Q,B,S,C$ at $V=20~{}{\rm fm^{3}}$, $T_{H}=170$ MeV and
$\gamma_{s}=\gamma_{c}=1$ is elucidated in Fig. 1. The immediate observation
is, while the dependence of $Z(Q,0,0,0)$ on the electric charge is rather
mild, $Z(0,B,0,0)$ and $Z(0,0,S,0)$ display a much stronger suppression as
baryon number or strangeness increases. The suppression of the canonical
partition function at finite quantum charge relative to $Z(0,0,0,0)$ can be
understood as the chemical factor of a corresponding charged hadron in a
completely neutral system. For instance, $Z(0,1,0,0)/Z(0,0,0,0)$ can be seen
as the corresponding chemical factor of an antibaryon. Canonical treatment of
SHM implies that, the creation of a baryon requires the accompanying creation
of an antibaryon in order to fulfill exact conservation of baryon number,
which is of course energy-expensive given that the mass of the lightest baryon
(i.e., the proton) $\sim 1$ GeV is already several times the typical
hadronization temperature $T_{H}$; in contrast, creation of an electrically-
charged hadron can be compensated by generating an oppositely-electrically-
charged pion to exactly conserve electric charge, which is less expensive in
terms of energy and explains the milder dependence of $Z(Q,0,0,0)$. This
argument generalizes to the case of $Z(0,0,0,C)$, where exact charm number
conservation requires the simultaneous production of a very massive charm and
anticharm pair, explaining the strongest decrease ($\sim$ two orders of
magnitude, shown in a logarithmic scale in the lowest panel) of the
$Z(0,0,0,C)$ upon each unit of increase of the charm number.

##  3 Canonical suppression of charm-hadron multiplicities

The ALICE measurement of $\Lambda_{c}/D^{0}$ in $\sqrt{s}=13$ TeV $pp$
collisions was performed in the mid-rapidity region, where the matter created
is nearly baryon-free [24]. Hence we impose vanishing quantum charges
$\vec{Q}=(Q,B,S,C)=(0,0,0,0)$ (i.e., a completely neutral system) in the CE-
SHM calculations in the following.

###  3.1 Selective canonical suppression

$hadron$ | CF | $\langle N_{j}\rangle^{CE}$ | $\langle N_{j}\rangle^{GCE}$
---|---|---|---
$\pi^{0}$ | 1.00000 | 1.13033 | 1.13033
$\pi^{+}$ | 0.89365 | 1.00146 | 1.12064
$K^{0}$ | 0.80905 | 0.30430 | 0.37612
$K^{+}$ | 0.80961 | 0.30917 | 0.38187
$p$ | 0.66607 | 0.076786 | 0.11528
$\Lambda$ | 0.66377 | 0.033375 | 0.050281
$D^{0}$ | 0.00520 | 0.0000031128 | 0.00059834
$D^{+}$ | 0.00518 | 0.0000030250 | 0.00058382
$D_{s}^{+}$ | 0.00485 | 0.0000016972 | 0.00035016
$\Lambda_{c}^{+}$ | 0.00458 | 0.0000006056 | 0.00013214
Table 1: Chemical factors (CF) and direct multiplicities of a selective set of
light and charm hadrons at $V=20~{}{\rm fm^{3}}$, $T_{H}=170$ MeV and
$\gamma_{s}=\gamma_{c}=1$ in CE-SHM and in GCE-SHM, respectively

In Table 1, we show the chemical factors and direct (i.e., without feeddowns)
multiplicities $\langle N_{j}\rangle^{CE}$ and $\langle N_{j}\rangle^{GCE}$ of
a selective set of light and charm hadrons calculated from CE-SHM and in GCE-
SHM, respectively at $V=20~{}{\rm fm^{3}}$, $T_{H}=170$ MeV and
$\gamma_{s}=\gamma_{c}=1$. Clearly the chemical factor of a given hadron is
obtained upon dividing $\langle N_{j}\rangle^{CE}$ by $\langle
N_{j}\rangle^{GCE}$ in such a completely neutral system, which displays in
general an increasing suppression for more massive charged hadrons. While the
stronger canonical suppression of $\pi^{+}$ relative to $\pi^{0}$ and of
$D^{+}$ relative to $D^{0}$ may be attributed to their carrying additional
electric charge, further suppression of $K^{+}$ relative to $\pi^{+}$ and
$D_{s}^{+}$ relative to $D^{+}$ is due to the requirement of strangeness
conservation.

It is interesting to note that, when dividing out the common effects from
canonical strangeness or charm suppression, the relative suppression (i.e.,
the ratio of their chemical factors) of $\Lambda$ with respect to $K^{0}$
(denoted as $\Lambda/K^{0}$) $\simeq 0.820$ and of the
$\Lambda_{c}^{+}/D^{+}\simeq 0.886$ turn out to be very similar and both
comparable to the $p/\pi^{+}\simeq 0.745$, delineating the pure effect of
canonical suppression owing to exact baryon-number conservation. This can be
more quantitatively understood from the selective baryon-number canonical
ensemble SHM, in which only baryon-number is subject to exact conservation in
a canonical treatment but other quantum charges are still treated grand-
canonically and thus conserved only on average [32]. This selective baryon-
number canonical ensemble SHM can be obtained by putting all $\phi$’s other
than $\phi_{B}$ to zero in the integrand of Eq. (3) and disregarding the
corresponding integrals, yielding

$\displaystyle Z(B)$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}d\phi_{B}e^{iB\phi_{B}}{\rm
exp}[\sum_{j}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}e^{-iB_{j}\phi_{B}}z_{j}]$
$\displaystyle=\prod_{j=M^{0},M^{+},M^{-}}{\rm
exp}(\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j})\times I_{B}(x_{B}),$ (10)

where the contribution from mesons ($B_{j}=0$) is factorized out and the
technique of separating the baryons ($B_{j}=+1$) and antibaryons ($B_{j}=-1$)
leading to Eq. (2.2) is used to arrive at the modified Bessel function of the
first kind, with
$x_{B}=\sum_{j=B,{\bar{B}}}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}$ being
the sum of multiplicities of all baryons and antibaryons in the grand-
canonical limit. Substituting this into Eq. (4), the primary mean number of
the $j$-th hadron from the baryon-number canonical ensemble SHM is

$\displaystyle\langle
N_{j}\rangle^{B-CE}=\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}\frac{I_{B-B_{j}}(x_{B})}{I_{B}(x_{B})}.$
(11)

Therefore, with $x_{B}=2.38$ computed from the hadron spectrum we have
employed at $V=20~{}{\rm fm^{3}}$, $T_{H}=170$ MeV and
$\gamma_{s}=\gamma_{c}=1$, the production of baryons and antibaryons in a
completely neutral system is suppressed by a factor of
$I_{1}(x_{B})/I_{0}(x_{B})\simeq 0.752$ as a result of the exact conservation
of baryon number, which roughly explains the relative suppression of
$p/\pi^{+}$, $\Lambda/K^{0}$ and $\Lambda_{c}^{+}/D^{+}$.

Similarly, when conserving only the charm number, the selective charm
canonical ensemble SHM gives

$\displaystyle\langle
N_{j}\rangle^{C-CE}=\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}\frac{I_{C-C_{j}}(x_{C})}{I_{C}(x_{C})},$
(12)

where $x_{C}=\sum_{j,C_{j}=\pm 1}\gamma_{s}^{N_{sj}}\gamma_{c}^{N_{cj}}z_{j}$
is the sum of multiplicities of all charm-hadrons in the grand-canonical
limit. With $x_{C}=0.01254$ in the present calculation, the pertinent chemical
factor $I_{1}(x_{C})/I_{0}(x_{C})\simeq 0.00627$ also roughly explains the
relative suppression of $D_{s}^{+}/K^{+}\simeq 0.00599$,
$\Lambda_{c}^{+}/p\simeq 0.00688$ (where the common effect from canonical
strangeness or baryon-number suppression is divided out) and
$D^{+}/\pi^{+}\simeq 0.00580$ (all computed from Table. 1).

###  3.2 System-size dependence of charm-hadron production

$CF$ | V=10 ${\rm fm}^{3}$ | 20 | 50 | 100 | 200
---|---|---|---|---|---
$D^{0}$ | 0.025877 | 0.066239 | 0.190294 | 0.373107 | 0.627886
$D^{+}$ | 0.025439 | 0.065891 | 0.190002 | 0.372841 | 0.627669
$D_{s}^{+}$ | 0.015805 | 0.053178 | 0.178586 | 0.362376 | 0.619125
$\Lambda_{c}^{+}$ | 0.016956 | 0.055485 | 0.182039 | 0.365923 | 0.622147
$\Xi_{c}^{+0}$ | 0.009884 | 0.042956 | 0.167943 | 0.352535 | 0.611073
$\Omega_{c}$ | 0.003495 | 0.022604 | 0.130312 | 0.312514 | 0.576383
$\Lambda_{c}^{+}/D^{0}$ | 0.655254 | 0.837649 | 0.956620 | 0.980745 | 0.990860
$D_{s}^{+}/D^{0}$ | 0.610774 | 0.802820 | 0.938474 | 0.971239 | 0.986047
Table 2: Chemical factors (CF) of ground-state charm-hadrons at $T_{H}=170$
MeV, $\gamma_{s}=0.6$, $\gamma_{c}=15$ and varying volumes. The ratios of the
CF’s of $\Lambda_{c}^{+}$ and $D^{0}$ and of $D_{s}^{+}$ and $D^{0}$ are also
summarized in the last two rows.

Having presented the analysis of the canonical suppression mechanisms, we are
now in a position to study the system-size dependence of the charm-hadron
yields in a completely neutral system via CE-SHM under more realistic
conditions. The strangeness saturation parameter is still fixed to be
$\gamma_{s}=0.6$ [13], as it has been shown to be insensitive to the system
size for charged particle multiplicity $dN_{\rm ch}/d\eta\leq 50$ [32] which
covers the range of our current interest. For the charm fugacity factor
$\gamma_{c}$, without further knowledge of the total charm cross section in
the present $\sqrt{s}=13$ TeV $pp$ collisions, we take a typical value of
$\sim 15$ determined from the charm canonical ensemble SHM in semi-central
heavy-ion collisions at the LHC energy [37]. Given the flavor hierarchy in the
operational hadronization temperature as suggested based on comparisons of
light- and strange-quark susceptibilities [38], and the agreement of the
susceptibilities computed from lattice with the same charm-baryon spectrum as
used here up to temperatures of $170$ MeV [27], we still take $T_{H}=170$ MeV
as our default value of hadronization temperature for the production of charm-
hadrons [13]. Uncertainty from using a possibly lower $T_{H}\sim 160$ MeV may
be (partially) bracketed by the varied branching ratios of the excited charm-
baryons decaying into $\Lambda_{c}^{+}$ [13].

In Table 2, chemical factors (CF) of ground-state charm-hadrons with
$T_{H}=170$ MeV, $\gamma_{s}=0.6$, $\gamma_{c}=15$ are displayed at varying
volumes. At small volumes, $D_{s}^{+}$ and $\Lambda_{c}^{+}$ have smaller
chemical factors than $D^{0}$ and $D^{+}$ due to additional canonical
strangeness and baryon-number suppression, respectively. The chemical factors
of charm-baryons $\Xi_{c}$ and $\Omega_{c}$ are further progressively smaller
as the strangeness content increases, reflecting the canonical suppression
from conserving all three (charm, strangeness and baryon number) quantum
charges. As volume increases, the canonical suppression effects brought about
by exact strangeness and baryon-number conservation are diminishing, and hence
at asymptotically large volume, the chemical factors of all charm-harons tend
to the same residual value that arises solely from their common canonical
charm suppression. This is also seen from the gradual increase of the relative
canonical suppression (i.e., the ratios of the chemical factors) of
$\Lambda_{c}^{+}/D^{0}$ and of $D_{s}^{+}/D^{0}$ toward unity from small to
large volumes (see the last two rows of Table. 2). While the absolute value of
the chemical factors does depend on $\gamma_{c}$, we have checked that the
ratio of the chemical factors such as $\Lambda_{c}^{+}/D^{0}$ is nearly
independent of $\gamma_{c}$, essentially guaranteeing that the fragmentation
fraction of each charm-hadron is robust against variation of the $\gamma_{c}$.
This can be understood from the fact that the $\Lambda_{c}^{+}/D^{0}$, with
the common canonical charm suppression divided out (canonical suppression due
to electric-charge conservation is negligible for very massive hadrons), can
be attributed to the canonical baryon-number suppression on $\Lambda_{c}^{+}$
as discussed in Sec. 3.1 and thus intrinsically determined by the pertinent
chemical factor $I_{1}(x_{B})/I_{0}(x_{B})$ in Eq. (11). But the charm-baryon
contribution (proportional to $\gamma_{c}$) to the total $x_{B}$ turns out to
be less than 2% even for $\gamma_{s}=15$, so that the variation of
$\gamma_{c}$ has little impact on the value of the relevant chemical factor.

$n_{j}~{}(\cdot 10^{-4}{\rm fm}^{-3}$) | V=10 ${\rm fm}^{3}$ | 20 | 50 | 100 | 200 | GCE
---|---|---|---|---|---|---
$D^{0}$ | 0.445553 | 1.148287 | 3.310131 | 6.495330 | 10.934662 | 17.420949
$D^{+}$ | 0.194705 | 0.503847 | 1.453016 | 2.851351 | 4.800262 | 7.647869
$D_{s}^{+}$ | 0.075040 | 0.252484 | 0.847910 | 1.720531 | 2.939551 | 4.747914
$\Lambda_{c}^{+}({\rm BR}50\%)$ | 0.126963 | 0.439135 | 1.497132 | 3.045487 | 5.207572 | 8.415360
$\Lambda_{c}^{+}({\rm BR}100\%)$ | 0.149573 | 0.519555 | 1.776775 | 3.617118 | 6.187127 | 10.001702
$\Xi_{c}^{+0}$ | 0.016539 | 0.071955 | 0.281624 | 0.591389 | 1.025276 | 1.678110
$\Omega_{c}$ | 0.000756 | 0.004889 | 0.028184 | 0.067592 | 0.124662 | 0.216283
$\Lambda_{c}^{+}/D^{0}({\rm BR}50\%)$ | 0.284956 | 0.382426 | 0.452288 | 0.468873 | 0.476244 | 0.483060
$\Lambda_{c}^{+}/D^{0}({\rm BR}100\%)$ | 0.335702 | 0.452461 | 0.536769 | 0.556880 | 0.565827 | 0.574119
$D_{s}^{+}/D^{0}$ | 0.168420 | 0.219879 | 0.256156 | 0.264887 | 0.268829 | 0.272540
Table 3: Total (i.e., with feeddowns) thermal densities of ground-state charm-
hadrons at $T_{H}=170$ MeV, $\gamma_{s}=0.6$, $\gamma_{c}=15$ and varying
volumes in the CE-SHM and also in the GCE-SHM (last column). The ratios of the
thermal densities of $\Lambda_{c}^{+}$ and $D^{0}$ and of $D_{s}^{+}$ and
$D^{0}$ are also summarized in the last three rows. For $\Lambda_{c}^{+}$, two
scenarios of the branching ratios (BR=$50\%$ vs. BR=$100\%$) of the RQM-
augmented excited $\Lambda_{c}$’s and $\Sigma_{c}$’s above $DN$ threshold
decaying to the ground state $\Lambda_{c}^{+}$ are compared.

We have computed the the primary production numbers of all charm-hadrons
within the spectrum employed (i.e., PDG charm-mesons, PDG charm-baryons, and
additional charm-baryons from RQM predictions) with Eq. (4), and then the
excited states are decayed into the ground state charm-hadrons with branching
ratios specified in PDG. For excited $\Lambda_{c}$’s and $\Sigma_{c}$’s above
the $DN$ threshold, we still assume two scenarios of BR=$100$% or BR=$50$% for
their branching ratios decaying to the ground state $\Lambda_{c}^{+}$ [13]. In
Table 3, the total (i.e., with feeddowns from higher states) thermal densities
(number divided by the corresponding volume) of ground-state charm-hadrons
computed from the CE-SHM with $T_{H}=170$ MeV, $\gamma_{s}=0.6$,
$\gamma_{c}=15$ are summarized at varying volumes, in comparison with the GCE-
SHM ones shown in the last column. An immediate observation is that the
thermal density of each charm-hadron is increasing from small to large volumes
toward their grand-canonical limit. The translated yield ratios of
$\Lambda_{c}^{+}/D^{0}$ and $D_{s}^{+}/D^{0}$ also demonstrate a marked
system-size dependence, representing a $\sim 40$% reduction from the largest
to the smallest volume as shown here, essentially due to the additional
canonical baryon-number and strangeness suppression on charm-baryons and
charm-strange mesons in small systems. One notes that, while their thermal
densities at $V=200~{}{\rm fm^{3}}$ are still far from the corresponding GCE-
SHM values, the yield ratios of both $\Lambda_{c}^{+}/D^{0}$ and
$D_{s}^{+}/D^{0}$ at this volume are already almost the same as the grand-
canonical values, simply because the effects of canonical baryon-number and
strangeness suppression already become vanishing and the residual canonical
charm suppression is common for all charm-hadrons.

##  4 Fragmenation and Decay Simulation

To calculate the charm-hadron $p_{T}$ spectra and then their ratios, we take
the charm quark $p_{t}$ spectrum generated by FONLL [28, 29] in $\sqrt{s}=13$
TeV $pp$ collisions, and then perform its fragmentation into all kinds of
charm-mesons and charm-baryons using the same procedure specified in [13]. The
fragmentation function [39] reads

$\displaystyle D_{c\rightarrow
H}(z)=N_{H}\frac{rz(1-z)^{2}}{[1-(1-r)z]^{6}}[6-18(1-2r)z$
$\displaystyle+(21-74r+68r^{2})z^{2}\qquad\quad$
$\displaystyle-2(1-r)(6-19r+18r^{2})z^{3}$
$\displaystyle\quad+3(1-r)^{2}(1-2r+2r^{2})z^{4}]\,$ (13)

where $z=p_{T}/p_{t}$ is the fraction of the hadron ($H$) momentum ($p_{T}$)
with respect to the quark momentum ($p_{t}$). For the parameter $r$ in Eq.
(4), we take the same value $r_{D^{0}}$=0.1 and $r_{\Lambda_{c}^{+}}$=0.16 as
used in [13], which were tuned to fit the slope of the measured $p_{T}$
spectra of the ground-state $D^{0}$ and $\Lambda_{c}^{+}$ in $\sqrt{s}=5.02$
TeV $pp$ collisions. The value of $r$ for other charm-mesons ($M$) follows
from the mass scaling:
$r_{M}/r_{D^{0}}=((m_{M}-m_{c})/m_{M})/((m_{D^{0}}-m_{c})/m_{D^{0}})$ [39],
where $m_{c}=1.5$ GeV is the charm-quark mass used in our calculations, and
likewise for charm baryons ($B$):
$r_{B}/r_{\Lambda_{c}^{+}}=((m_{B}-m_{c})/m_{B})/((m_{\Lambda_{c}^{+}}-m_{c})/m_{\Lambda_{c}^{+}})$.

The normalization coefficient, $N_{H}$, for a given hadron, is determined by
requiring its $p_{T}$ integrated yield to be equal to the fraction of its CE-
SHM density at a given volume relative to the total density of all charm-
hadrons; that is, fragmentation fractions are determined by the thermal
weights as computed from CE-SHM. To associate the total mean number of the
ground-state charm-hadrons $\langle N_{j}^{\rm tot}\rangle^{CE}$ computed in
the CE-SHM with the rapidity density $dN_{j}/dy$ as experimentally measured in
high-energy $pp$ collisions, one still need to establish the connection
between the correlation volume $V_{c}$ entering the CE-SHM calculation in
which quantum charges are strictly conserved and the $dV/dy$ corresponding to
one unit of rapidity. Causal argument suggests that $V_{c}$ may correspond to
a few units of rapidity [40]; canonical thermal fits also suggests that
$V_{c}=(1$-$3)dV/dy$ and further $dV/dy\simeq 2.4dN_{\rm ch}/d\eta~{}({\rm
fm^{3}})$ [32]. In the following, $V_{c}=2dV/dy$ is chosen to demonstrate our
final results, which has been found to yield best fit to the splitting of
$\Lambda_{c}/D^{0}$ for the two bins of $dN_{\rm ch}/d\eta$ under
consideration.

Each charm hadron produced from fragmentation is then decayed into ground-
state particles by assuming a constant matrix element, with the decay
kinematics solely determined by phase space and the branching ratios as
discussed in Sec. 3.2. This way, the final total $p_{T}$-spectra (their
respective absolute value associated with the total charm cross section is not
our primary concern in the present study) of the ground state charm-hadrons
are obtained and the $p_{T}$-dependent ratios of $\Lambda_{c}/D^{0}$ and
$D_{s}^{+}/D^{0}$ are generated, as demonstrated in Fig. 2 with
$V_{c}=2dV/dy$, with the uncertainty band in the $\Lambda_{c}/D^{0}$
encompassing the assumed range of the branching ratios (BR$=50$%-$100$%) of
the RQM-augmented excited $\Lambda_{c}$’s and $\Sigma_{c}$’s decaying into the
ground state $\Lambda_{c}^{+}$. We have also checked that varying
$V_{c}=(2$-$3)dV/dy$ but fixing BR=100% would generate almost the same central
value and similar uncertainty band for the $\Lambda_{c}/D^{0}$.

Figure 2: (Upper panel) The $\Lambda_{c}^{+}/D^{0}$ ratio computed from
fragmentation (and decay) simulations with CE-SHM weights determined at
$T_{H}=170$ MeV, $\gamma_{s}=0.6$, $\gamma_{c}=15$ for two bins of charged
particle multiplicity ($dN_{\rm ch}/d\eta=3.9$ and $44.0$) in $\sqrt{s}=13$
TeV $pp$ collisions, in comparison with ALICE measurement [24]. The
correlation volume is taken to be $V_{c}=2dV/dy$ and its further relation with
$dN_{\rm ch}/d\eta$ is quoted from [32]. The uncertainty band refers to the
variation of the branching ratios (BR=$50$%-$100$%) of the RQM-augmented
excited $\Lambda_{c}$’s and $\Sigma_{c}$’s decaying to the ground state
$\Lambda_{c}^{+}$. For $dN_{\rm ch}/d\eta=44.0$, result from using only PDG
listings of charm-baryons is also shown. (Lower panel) The same result for the
$D_{s}^{+}/D^{0}$ ratio.

From Fig. 2, the role of the RQM-augmented charm-baryons is reaffirmed that
generate a substantial enhancement of the $\Lambda_{c}/D^{0}$ relative to the
case with only PDG listings of charm-baryons through their feeddowns to the
ground state $\Lambda_{c}^{+}$ [13]. As the main result of the present study,
the splitting of the $\Lambda_{c}/D^{0}$ between two bins of $dN_{\rm
ch}/d\eta$ as indicated by the ALICE measurements is roughly described by our
CE-SHM calculations, due to the additional canonical baryon-number suppression
on charm-baryons that becomes more pronounced toward smaller system size. Yet
the computed curve shows a significant deficiency relative to ALICE data at
intermediate $p_{T}\simeq 2$-$5$ GeV in the $dN_{\rm ch}/d\eta=44.0$ bin,
which might be remedied by possible radial flow in such high multiplicity
events hardening the charm-baryon spectrum more than the charm-meson one.
Similarly, the calculated $D_{s}^{+}/D^{0}$ also exhibits a significant
splitting between the two considered bins of $dN_{\rm ch}/d\eta$, which is
attributed to the additional canonical strangeness suppression on charm-
strange mesons.

##  5 Summary

We have computed the charged-particle multiplicity dependent hadro-chemistry
of charm-hadrons (as represented by the $\Lambda_{c}/D^{0}$ and
$D_{s}^{+}/D^{0}$ ratios) in high-energy $pp$ collisions within the approach
of the canonical ensemble statistical hadronization model (CE-SHM). Taking
full particle listings of PDG complemented by additional charm-baryons from
relativistic quark model predictions, the canonical partition function has
been evaluated and the chemical factors of typical hadrons that measure the
pertinent canonical suppression computed. In particular, the mechanisms
underlying the relative suppression between typical pairs of hadrons have been
examined in connection with the selective canonical ensemble SHM that
implements strict conservation of a single quantum charge (charm or baryon
number) that distinguishes between the particles within the pair.

The chemical factors of various charm-hadrons have been shown to exhibit a
pronounced system-size dependence that translates to charged-particle
multiplicity dependence of charm-hadron ratios, highlighting the importance of
strict conservation of quantum charges in small systems. Taking thus computed
thermal densities of charm-hadrons as fragmentation weights, fragmentation and
decay simulations yield $p_{T}$-dependent $\Lambda_{c}/D^{0}$ and
$D_{s}/D^{0}$ ratios at different charged-particle multiplicities in fair
agreement with ALICE measurements. The pronounced splitting of
$\Lambda_{c}/D^{0}$ and $D_{s}/D^{0}$ between different multiplicity bins is
attributed to the stronger canonical suppression of $\Lambda_{c}$ and $D_{s}$
in consequence of strict baryon-number and strangeness conservation (besides
their common canonical charm suppression), respectively, toward smaller
system-size.

Acknowledgments: This work was supported by the NSFC under grant 12075122. We
are indebted to Ralf Rapp for useful communications.

## References

  * [1] R. Placakyte [H1and for the ZEUS Collaboration], arXiv:1111.5452 [hep-ph].
  * [2] M. Lisovyi, A. Verbytskyi and O. Zenaiev, Eur. Phys. J. C 76 (2016) no.7, 397.
  * [3] S. Acharya et al. [ALICE Collaboration], JHEP 1804 (2018) 108.
  * [4] Y. Amhis et al. [HFLAV Collaboration], Eur. Phys. J. C 77 (2017) no.12, 895.
  * [5] R. Maciula and A. Szczurek, Phys. Rev. D 98 (2018) no.1, 014016.
  * [6] B. A. Kniehl, G. Kramer, I. Schienbein and H. Spiesberger, Phys. Rev. D 101 (2020) no.11, 114021.
  * [7] S. Acharya et al. [ALICE Collaboration], Phys. Lett. B 793 (2019) 212.
  * [8] A. M. Sirunyan et al. [CMS Collaboration], Phys. Lett. B 803 (2020) 135328.
  * [9] M. He and R. Rapp, Phys. Rev. Lett. 124 (2020) no.4, 042301.
  * [10] C. Bierlich and J. R. Christiansen, Phys. Rev. D 92 (2015) no.9, 094010.
  * [11] P. Braun-Munzinger, K. Redlich and J. Stachel, “Particle production in heavy ion collisions,” In *Hwa, R.C. (ed.) et al.: Quark gluon plasma* 491-599, [nucl-th/0304013].
  * [12] A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, Nature 561 (2018) no.7723, 321.
  * [13] M. He and R. Rapp, Phys. Lett. B 795 (2019) 117.
  * [14] D. Ebert, R. N. Faustov and V. O. Galkin, Phys. Rev. D 84 (2011) 014025.
  * [15] I. Kuznetsova and J. Rafelski, Eur. Phys. J. C 51 (2007) 113.
  * [16] A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, Phys. Lett. B 659 (2008) 149.
  * [17] J. Adam et al. [ALICE Collaboration], JHEP 1509 (2015) 148.
  * [18] S. Acharya et al. [ALICE Collaboration], Phys. Lett. B 810 (2020) 135758.
  * [19] A. M. Sirunyan et al. [CMS Collaboration], arXiv:2007.04277 [hep-ex].
  * [20] R. Aaij et al. [LHCb Collaboration], arXiv:2009.06619 [hep-ex].
  * [21] S. Deb, R. Sahoo, D. Thakur, S. Tripathy and A. Khuntia, arXiv:2010.08984 [nucl-ex].
  * [22] V. Khachatryan et al. [CMS Collaboration], Phys. Lett. B 765 (2017) 193.
  * [23] J. Adam et al. [ALICE Collaboration], Nature Phys. 13 (2017) 535.
  * [24] C. Hill, for ALICE Collaboration, HP2020 proceedings
  * [25] M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98, no.3 (2018) 030001.
  * [26] P. Madanagopalan, R. G. Edwards, N. Mathur and M. J. Peardon, PoS LATTICE 2014 (2015) 084.
  * [27] A. Bazavov et al., Phys. Lett. B 737 (2014) 210.
  * [28] M. Cacciari, M. Greco and P. Nason, JHEP 9805 (1998) 007.
  * [29] M. Cacciari, S. Frixione, N. Houdeau, M. L. Mangano, P. Nason and G. Ridolfi, JHEP 1210 (2012) 137.
  * [30] J. Rafelski and M. Danos, Phys. Lett. 97B (1980) 279.
  * [31] R. Hagedorn and K. Redlich, Z. Phys. C 27 (1985) 541.
  * [32] V. Vovchenko, B. Donigus and H. Stoecker, Phys. Rev. C 100 (2019) no.5, 054906.
  * [33] F. Becattini, Z. Phys. C 69 (1996) no.3, 485.
  * [34] F. Becattini and U. W. Heinz, Z. Phys. C 76 (1997) 269.
  * [35] F. Becattini, arXiv:0901.3643 [hep-ph].
  * [36] A. Keranen and F. Becattini, Phys. Rev. C 65 (2002) 044901, Erratum: [Phys. Rev. C 68 (2003) 059901].
  * [37] A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, Nucl. Phys. A 789 (2007) 334.
  * [38] R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz and C. Ratti, Phys. Rev. Lett. 111 (2013) 202302.
  * [39] E. Braaten, K. m. Cheung, S. Fleming and T. C. Yuan, Phys. Rev. D 51, 4819 (1995).
  * [40] P. Castorina and H. Satz, Int. J. Mod. Phys. E 23 (2014) no.4, 1450019.


* Zhang et al. (2020) Fan Zhang, Conggai Li, Ying Zhang, Lu Qin, and Wenjie Zhang. 2020.  Finding Critical Users in Social Communities: The Collapsed Core and Truss Problems.  _IEEE Transactions on Knowledge and Data Engineering_ 32, 1 (2020), 78–91.
  * Zhang et al. (2017b) Fan Zhang, Wenjie Zhang, Ying Zhang, Lu Qin, and Xuemin Lin. 2017b.  OLAK: an efficient algorithm to prevent unraveling in social networks.  _Proceedings of the VLDB Endowment_ 10, 6 (2017), 649–660.
  * Zhang et al. (2017a) Fan Zhang, Ying Zhang, Lu Qin, Wenjie Zhang, and Xuemin Lin. 2017a.  Finding critical users for social network engagement: The collapsed k-core problem. In _AAAI_.
  * Zhang et al. (2018b) Fan Zhang, Ying Zhang, Lu Qin, Wenjie Zhang, and Xuemin Lin. 2018b.  Efficiently reinforcing social networks over user engagement and tie strength. In _ICDE_.
  * Zhang et al. (2019) Na Zhang, Xuefeng Guan, Jun Cao, Xinglei Wang, and Huayi Wu. 2019.  Wavelet-HST: A wavelet-based higher-order spatio-temporal framework for urban traffic speed prediction.  _IEEE Access_ 7 (2019), 118446–118458.
  * Zhang and Yu (2019) Yikai Zhang and Jeffrey Xu Yu. 2019.  Unboundedness and efficiency of truss maintenance in evolving graphs. In _SIGMOD_.
  * Zhao et al. (2020) Jun Zhao, Renjie Sun, Qiuyu Zhu, Xiaoyang Wang, and Chen Chen. 2020.  Community identification in signed networks: a k-truss based model. In _CIKM_.
  * Zhao et al. (2021) Kangfei Zhao, Zhiwei Zhang, Yu Rong, Jeffrey Xu Yu, and Junzhou Huang. 2021.  Finding critical users in social communities via graph convolutions.  _IEEE Transactions on Knowledge and Data Engineering_ 35, 1 (2021), 456–468.
  * Zheng et al. (2017) Zibin Zheng, Fanghua Ye, Rong-Hua Li, Guohui Ling, and Tan Jin. 2017.  Finding weighted k-truss communities in large networks.  _Information Sciences_ 417 (2017), 344–360.
  * Zhou et al. (2019) Yaoming Zhou, Junwei Wang, and Hai Yang. 2019.  Resilience of transportation systems: concepts and comprehensive review.  _IEEE Transactions on Intelligent Transportation Systems]_ 20, 12 (2019), 4262–4276.
  * Zhou et al. (2022) Zhongxin Zhou, Wenchao Zhang, Fan Zhang, Deming Chu, and Binghao Li. 2022.  VEK: a vertex-oriented approach for edge k-core problem.  _World Wide Web_ 25, 2 (2022), 723–740.
  * Zhu et al. (2022) Difeng Zhu, Guojiang Shen, Jingjing Chen, Wenfeng Zhou, and Xiangjie Kong. 2022.  A higher-order motif-based spatiotemporal graph imputation approach for transportation networks.  _Wireless Communications and Mobile Computing_ 2022 (2022), 1–16.
  * Zhu et al. (2018) Weijie Zhu, Chen Chen, Xiaoyang Wang, and Xuemin Lin. 2018\.  K-core minimization: An edge manipulation approach. In _CIKM_.
  * Zhu et al. (2019) Weijie Zhu, Mengqi Zhang, Chen Chen, Xiaoyang Wang, Fan Zhang, and Xuemin Lin. 2019\.  Pivotal Relationship Identification: The K-Truss Minimization Problem.. In _IJCAI_.